kogad blog

進化計算とかプログラミングとか，勉強したことを書きます．

Julia の多重ディスパッチを使って木構造で表現される数式を評価する

Julia プログラミング 遺伝的プログラミング

はじめに

以前，遺伝的プログラミングを実装したとき，木構造で表現された数式を評価する処理を，再帰を使って書いたことがありました．

よく考えると，多重ディスパッチも使うとスッキリ書けそうだなと思ったので，試しにやってみました．

実装

演算子，変数，定数となる頂点の型（構造体）を定義して，それぞれの型ごとに関数を定義してやるだけです．

  • 演算子の頂点であれば，子の頂点に演算子を適用する
  • 変数の頂点であれば，その頂点が持つ変数に対応する数値を取り出す
  • 定数の頂点であれば，その定数をそのまま返す

変数と，変数に代入する値はディクショナリで管理しています．

gist711c8bf24be47dfcfbceddd469e5ff71

遺伝的アルゴリズムでクラスタリング（GA-clustering）

遺伝的アルゴリズム プログラミング

はじめに

この記事では遺伝的アルゴリズム（Genetic Algorithm; GA）を用いてクラスタリングを行う GA-clustering についてまとめます．

以下の論文を参考にして書いています．この記事になにかおかしい点があれば，この論文を参照してください．

U.Maulik, S.Bandyopadhyay, Genetic algorithm-based clustering technique, Pattern Recog. 33 (2000) 1455–1465

CiteSeerXで読めます (http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.19.1878&rep=rep1&type=pdf)．

また，関連手法としてもっとうまい手法もあります（が，とりあえずはこれを紹介させて欲しい）．

GA-clustering とは

GA-clusteringとは，遺伝的アルゴリズムと K-means 法を組み合わせたクラスタリング手法です．

K-means 法の，局所最適解にハマってしまうことがあるという欠点を，遺伝的アルゴリズムで補おうというアイディアです．

K-means 法

Wikipedia 参照

ja.wikipedia.org

GA-clustering の流れ

おおまかな流れは基本的な遺伝的アルゴリズムと同じで，初期集団を生成してから評価→選択→交叉→突然変異を終了条件を満たすまで繰り返します．

各個体は，すべてのクラスタの重心を並べた配列とします．

以下，各操作の中身を見ていきます．

初期集団の生成

クラスタリングしたいデータを K 個のクラスタにランダムに割り振って，それぞれのクラスタの重心を計算し，個体を生成します．

これを集団の個体数 N だけ繰り返し，初期集団とします．

個体の評価

適応度を計算する前に，K-means 法を1ステップ分だけ実行して重心を更新しておきます．

クラスタリングの評価基準として，各クラスタ内の各点と重心の距離の和を定義します．

各クラスタを C_1, C_2, \ldots, C_K，その重心を \boldsymbol{z}_1, \boldsymbol{z}_2, \ldots, \boldsymbol{z}_K とすると，

$$ \mathscr{M}(C_1, C_2, \ldots, C_K) = \sum_{i = 1}^{K} \sum_{\boldsymbol{x}_i \in C_i} || \boldsymbol{x}_i - \boldsymbol{z}_i || $$

\boldsymbol{z}_i はクラスタ C_i の重心です．

この \mathscr{M} が小さいほど良いクラスタが得られたと考えます．

したがって，適応度は\displaystyle\frac{1}{\mathscr{M}(C_1, C_2, \ldots, C_K)} とします．

選択

ルーレット選択です．

交叉

交叉率 p_c で二点交叉を適用します．

突然変異

突然変異率は p_m とします．

突然変異を，v を遺伝子の位置（座標），\delta を [0, 1] の一様乱数として，

$$ \begin{align} v \pm 2 \delta v, ~( v \neq 0 ) \\ v \pm 2 \delta, ~( v = 0 )~ \end{align} $$

と定義します．\pm は等確率でどちらかが選ばれます．

実装

ここに置いときます．

github.com

わりと汚いですけど，勘弁してくださいな．

（ちなみに，plotのところはmarker_zを使ってあげるとスッキリ書けます．でも，いろいろ試したけど色が気に食わなかったので......）

試しにクラスタリング

適当にデータを作ってクラスタリングしてみました．

個体数200，世代数200，交叉率0.9，突然変異率0.001です．

左が元のデータ，真ん中が GA-clustering によるクラスタリング，右が K-means 法によるクラスタリングです．

f:id:kogad:20190504211410p:plain

K-means 法より GA-clustering のほうが良さげな感じ出てますね．

何回か実行してみるとK-means 法でもいい感じの結果は出ますが，GAの方はほぼ毎回いい感じになるはずです．

すべてのビット列を出力する

プログラミング Julia

ある長さのすべてのビット列のパターンが欲しい場面があったんだけど，意外とハマってしまった． 正直とても簡単だけど，怪電波を受信してしまって訳わからんことをしてしまったので，自戒の念を込めて記事にする．

長さ L のすべてのビット列が欲しければ，1 ~ L の整数をビット列に変換して端から L ビットだけとってくればよい．

Julia であれば，

L = 4
for i = 0:2^L-1
    str = bitstring(i)
    println(str[length(str)-L+1:length(str)])
end

とすれば，

0000
0001
0010
0011
0100
0101
0110
0111
1000
1001
1010
1011
1100
1101
1110
1111

のような感じになる．

bitstring()の返り値の型は String なので，整数の配列とかで欲しければ適当に parse() してやればよい．

もっと簡単な方法があるような気もする．

追記 (2021/02/07) digits(n; base, pad) がの方が簡単でよいです．

Julia で機械学習 〜 Flux.jl で XOR を学習させてみる

機械学習 Julia プログラミング

はじめに

　この記事では， Julia の機械学習用フレームワークの Flux.jl を使って，簡単なニューラルネットワークに XOR を学習させてみます．

　Flux.jl の詳細についてはドキュメントを参照してください．また，この記事はドキュメントを参考にしています．

fluxml.ai

コード

　とりあえずこれがコードです．

using Flux
using Flux: @epochs, mse, train!

m = Chain(Dense(2,4, tanh),
          Dense(4,1,σ))

loss(x,y) = mse(m(x), y)
ps = params(m)

xs = [[0,0], [0,1], [1,0], [1,1]]
ys = [0,1,1,0]
data = zip(xs, ys)

η = 0.1
opt = Descent(η)

@epochs 10000 train!(loss, ps, data, opt)

# 結果
for i = 1:4
    println(m(xs[i]))
end

解説

m = Chain(Dense(2,4, tanh),
          Dense(4,1,σ))

　この部分でニューラルネットワークのモデルを定義しています．

Dense() は y = Wx + b の入出力を持つ層を意味します．引数は順番に，入力数・出力数・活性化関数です．

　 Chain() は引数で渡した層・関数を繋げてくれる関数です． Flux.jl で定義されている層の他に，自作の関数を用いることもできます．

σ はシグモイド関数です．Flux.jl で定義されています．

　　

ps = params(m)

　params()はモデルを引数とし，そのモデルのパラメータをまとめたものを返してくれます．

　　

xs = [[0,0], [0,1], [1,0], [1,1]]
ys = [0,1,1,0]
data = zip(xs, ys)

　データセットです．今回はXORの入出力です．

　　

η = 0.1
opt = Descent(η)

　ここではパラメータの最適化手法を定義しています．

Descent() は勾配降下法で，η は学習率です．

ちなみに，Flux.jl では勾配降下法以外にも，ADAM等が用意されています．

　　

@epochs 10000 train!(loss, ps, data, opt)

　この部分で学習を行います． @epochs マクロで指定したエポック数だけ train!() を実行します．

引数で指定した損失関数・パラメータ・データセット・最適化手法に対し学習を行います．

結果

　結果はこんな感じになりました．

Float32[0.00527357] (tracked)
Float32[0.986211] (tracked)
Float32[0.986653] (tracked)
Float32[0.0163712] (tracked)

　 この出力に対応する入力は，上から順に[0, 0], [0, 1], [1, 0], [1, 1] なので，ちゃんと学習できているみたいです．

(tracked) は Flux.jl によるものです．詳しくはドキュメントを．

　

　

スキーマ・ビルディングブロック・リンケージ（遺伝的アルゴリズム）

遺伝的アルゴリズム

はじめに

　この記事では，遺伝的アルゴリズムの中で重要な概念であるスキーマ・ビルディングブロック・リンケージについてまとめます．これらは，より性能の良い遺伝的アルゴリズムを考える上で非常に大切です．

目次

スキーマ

　スキーマは，遺伝子を表現する文字列と文字 * からなる個体と同じ長さの文字列で，以下で説明するように複数の個体を表します．

　例えば遺伝子がビット列で表現される場合，スキーマは \{ 0, 1, * \} の3つの文字からなる文字列です． * は 0, 1の両方の文字を表します．このとき，次のスキーマ

H = *0**1

は，以下の個体すべてを表します．

\{ 00001, 00011, 00101, 00111, 10001, 10011, 10101, 10111 \}

　ここで，スキーマの定義長と次数を定義します．

  • 定義長: スキーマの * 以外の文字のうち，左端にある文字と右端にある文字の距離
  • 次数: スキーマの * 以外の文字の数

上で示したスキーマ H = *0**1 の定義長は3，次数は2となります．

ビルディングブロック

　スキーマのうち，次数が小さい・定義長が短い・平均適応度が高いものをビルディングブロックと呼びます．平均適応度は，スキーマが表す個体の適応度の平均です．

　遺伝的アルゴリズムでは，ビルディングブロックをうまく処理することでより良い解を求めることができます．

　ビルディングブロックの例を１つ示します． 個体を20 bit のビット列 s で表し，その適応度を4 bit のトラップ関数を5つ並べた評価関数

\begin{equation} f(s) = \sum_{i = 1}^{4} {\rm{trap}}_5(u_i) \end{equation}

\begin{equation} \mathrm{trap}_5(u) = \begin{cases} 3 - u & ( u = 0, 1, 2, 3) \\ 4 & ( u = 4) \end{cases} \end{equation}

で計算するものとします．u_i は文字列 s を4等分した文字列中の1の数です．

　このとき，5つのトラップ関数のそれぞれの最適解のみを固定したスキーマ

\begin{align}　 H_1 &= 1111 ~ **** ~ **** ~ **** ~ **** \\ H_2 &= **** ~ 1111 ~ **** ~ **** ~ **** \\ H_3 &= **** ~ **** ~ 1111 ~ **** ~ **** \\ H_4 &= **** ~ **** ~ **** ~ 1111 ~ **** \\ H_5 &= **** ~ **** ~ **** ~ ****\ ~ 1111 \end{align}

はビルディングブロックと考えられます．これらのビルディングブロックを組み合わせることで，最適解

\begin{equation} s^{*} = 1111~1111~1111~1111~1111 \end{equation}

を得ることができます．

　遺伝的アルゴリズムでは，最適解につながるビルディングブロックに適合する個体を増やしていき，それらを交叉によって組み合わせることで良い解を求めることができます．逆に，交叉によってビルディングブロックが破壊されてしまうと解の探索がうまく進みません．

　ビルディングブロックが破壊されてしまう簡単な例として，上記のビルディングブロック H_2 を含む個体の一点交叉を考えます．

\begin{align} \require{color} 1001 ~ \textcolor{blue}{11} | \textcolor{blue}{11} ~ 0011 ~ 1010 ~ 0101 \\ 1101 ~ 00 | 00 ~ 1100 ~ 1001 ~ 1011 \end{align}

個体中に示した | で交叉させると，

\begin{align} 1011 ~ \textcolor{red}{1100} ~ 1100 ~ 1001 ~ 1011 \\ 1101 ~ 0011 ~ 0011 ~ 1010 ~ 0101 \end{align}

となり，ビルディングブロック H_2 が破壊されてしまいます（ちなみに，上記の2個体の場合，ビルディングブロックを作り出すことができる交叉点も存在します）．

リンケージ

　リンケージとは，複数の遺伝子座が集まってビルディングブロックを構成する傾向のことです．リンケージが分かれば，どの遺伝子座がビルディングブロックを構成しているかがわかります．ビルディングブロックを構成している遺伝子座がわかれば，ビルディングブロックを壊さないような交叉手法を適用することができます．

　 　例として，ビルディングブロックの例で扱った評価関数の部分関数（トラップ関数）それぞれに寄与する遺伝子座の集合をリンケージと考えることができます．遺伝子座の番号を左から1, ２, ......, 20とすると，{1, 2, 3, 4}, {5, 6, 7, 8}, {9, 10, 11, 12}, {13, 14, 15, 16}, {17, 18, 19, 20} がリンケージとなります．

　十分な数の初期個体群をランダムに生成することにより集団内にビルディングブロックを確保し，ビルディングブロックを壊さないように交叉を行うことで効率の良い解の探索ができます．初期集団にビルディングブロックが含まれていない場合でも，突然変異を適用することでビルディングブロックの供給が望めます．

　リンケージの長さが k ビットなら，O(k 2^k ) 程度の初期集団を確保すれば十分と考えられます．

参考

棟朝雅晴『遺伝的アルゴリズムーその理論と先端的手法ー』[POD版] , 森北出版

CGA（Compact GA）〜遺伝子の分布を推定する遺伝的アルゴリズム

遺伝的アルゴリズム プログラミング Julia

目次

はじめに

　CGA（Compact GA）とは，確率ベクトルから個体を生成・評価し確率ベクトルを修正していくことで最適化を行う遺伝的アルゴリズム（GA）の1つです．

　CGA の特徴は，同時に存在する個体は2つのみであるという点です．そのため，他の GA と比べてメモリ消費量が非常に少なくなります．

　この記事では CGA のアルゴリズムを紹介し，実際に簡単な最適化問題に適用してみたいと思います．

アルゴリズム

　CGAでは交叉や突然変異といった遺伝的操作は行わず，確率ベクトルを元に2つの個体を生成し，それらの適応度を確率ベクトルにフィードバックしていくことで最適化を進めていきます．以下ににその手順を示します．

  1. 確率ベクトル p に対し， i 番目の要素 p[i] = 0.5 （ i = 1, 2, ..., L ) として初期化する
  2. 確率 p[i] より，2つの固体 a, b を生成する
  3. 個体 a, b を評価し，評価値の高い方を winner， 低い方を loserとする
  4. winner, loser の i 番目の遺伝子 winner[i], loser[i] について， winner[i] ≠ loser[i] の場合，次のようにして確率ベクトルを更新する
    1. winner[i] = 1 の場合， p[i] = p[i] + 1/n
    2. winner[i] = 0 の場合， p[i] = p[i] - 1/n
  5. p の全ての要素に対し， p[i] = 1 または p[i] = 0 となれば終了．そうでなければ2へ戻る

　L は個体の長さ，n は GA における集団のサイズを想定した定数です．CGA は集団サイズ n のGAをシミュレートしていることになります．

実装

　ここでは100ビットの OneMax 問題を解きます．100ビットなので個体の長さ L = 100 とし，集団サイズは n = 100 としました．

　Julia (1.1.0) で実装しました．上記のアルゴリズムを素直に書いただけのプログラムです．

mutable struct Individual
    chrom::Array{Int64,1}
    fitness::Float64
end

function fitnessFunc(chrom::Array{Int64,1})
    sum(chrom)
end

function isConverged(prob::Array{Float64,1})
    for i = 1:length(prob)
        if 0 < prob[i] < 1 
            return false
        end
    end
    return true
end

const CHROM_LENGTH = 100
const POPULATION_SIZE = 100

function main()
    prob = Array{Float64, 1}(undef, CHROM_LENGTH)

    fill!(prob, 0.5)

    ind_a = Individual(Array{Int64,1}(undef, CHROM_LENGTH), 0.0)
    ind_b = Individual(Array{Int64,1}(undef, CHROM_LENGTH), 0.0)

    gen = 0
    while !isConverged(prob)
        gen += 1
        println(gen)

        for i = 1:CHROM_LENGTH
            if prob[i] > rand()
                ind_a.chrom[i] = 1
            else
                ind_a.chrom[i] = 0
            end

            if prob[i] > rand()
                ind_b.chrom[i] = 1
            else
                ind_b.chrom[i] = 0
            end
        end

        ind_a.fitness = fitnessFunc(ind_a.chrom)
        ind_b.fitness = fitnessFunc(ind_b.chrom)

        if ind_a.fitness >= ind_b.fitness
            winner = ind_a
            loser = ind_b
        else
            winner = ind_b
            loser = ind_a
        end

        for i = 1:CHROM_LENGTH
            if winner.chrom[i] != loser.chrom[i]
                if winner.chrom[i] == 1
                    prob[i] += 1/POPULATION_SIZE
                else
                    prob[i] -= 1/POPULATION_SIZE
                end
            end
        end
    end

    println("gen: $(gen), prob: $(round.(prob))")
end

main()

　以下が実行結果です．

gen: 1560, prob: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]

　1560世代で全てのビットが1になりました．数回実行してみましたが，だいたい1500世代前後で収束します．

　集団サイズをより小さくすると，ある程度のサイズまではより早く収束しました（だいたい40くらいまで）． 集団サイズ大きくすると世代数も大きくなりますが，一応収束しました．（集団サイズ1000000で試してみると19425050世代で収束しました）

　評価関数を書き換えればいろんな問題に適用できるので，パラメータと合わせていろいろ変えて試してみると面白いと思います．

参考

[1] 棟朝雅晴『遺伝的アルゴリズムーその理論と先端的手法ー』[POD版] , 森北出版

[2] G.R. Harik , F.G. Lobo , D.E. Goldberg. The compact genetic algorithm. IEEE Transactions on Evolutionary Computation ( Volume: 3 , Issue: 4 , Nov 1999 ) p.287-297. 1999