分布推定アルゴリズム
分布推定アルゴリズム。遺伝的アルゴリズムを改良した物です。個体の集合を交叉・突然変異させるのではなく、個体の生成確率を進化させます。最適化問題のアルゴリズムです。以下、自分へのメモです。わかったことが増えたら追記するかも。
ビットストリング
計算量に関しては、ビット数をn、反復数をTとしています。
- Population-Based Incremental Learning (PBIL)
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.61.8554
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.44.5424
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.43.1108
- Population-based incremental learning - Wikipedia
- http://mikilab.doshisha.ac.jp/dia/research/person/msano/reports/41/report41.html
- 1994年。これが全ての始まり。シンプルGAより成績良し。アルゴリズムも単純。Negative Learning Rate について、遺伝的アルゴリズム - その理論と先端的手法でふれられていないです。計算量は。
- 強化学習の追跡手法(強化学習のhttp://www.cs.ualberta.ca/~sutton/book/ebook/node23.html)と似た手法です。追跡手法に、Negative 学習率と突然変異が追加になっています。
- 探索に関する基本戦略は、
- 最初は広く浅く、徐々に深掘り。
- 時にはランダム探索。
- 負例からも学習。
- Hill Climbing with Learning (HCwL)
- 単変量周辺分布アルゴリズム (Univariate Marginal Distribution Algorithm; UMDA)
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.7.7030
- 1996年。PBILをよりシンプルにしたもの。計算量は。
- Mutual Information Maximimization for Input Clustering (MIMIC)
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.89.4429
- 1997年。2変数の条件付き確率を元に、個体生成確率を決めます。どのビットから順番に決めるかは、条件付きエントロピーの総和が最小になる順列を貪欲法で決めます。計算量は。
- Compact Genetic Algorithm (cGA)
- 二変量周辺分布アルゴリズム (Bivariate Marginal Distribution Algorithm; BMDA)
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.55.1151
- 1999年。カイ二乗を使い、変数間の関連を求め、カイ二乗による最大全域木を作り、依存関係の木を作り、依存関係にある物は条件付き確率を使います。カイ二乗ではなく、MIMICのようにエントロピー()を使った方が良かったのでは?エントロピー(相互情報量)は足し算できますが、カイ二乗はしてはいけないはず。計算量は最小全域木にプリム法を使えばですが、論文の方法だと。
- Edmonds' algorithm - Wikipediaによると、エドモンドス法を使えば、有向グラフの最小全域木も無向グラフの時と同じ計算量だそうで、カイ二乗や相互情報量ではなく、条件付きエントロピーを使い、有向グラフとして問題を解いた方が結果が良くなると思います。BOAはベイジアンネットワークを使っていますが、計算量が多すぎると思います。
- ベイジアン最適化アルゴリズム (Bayesian Optimization Algorithm; BOA)
- http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.52.2148
- ベイジアン最適化アルゴリズムに関する調査
- http://mikilab.doshisha.ac.jp/dia/research/person/suyara/report2002/8_30/index.html
- 実装は、http://www.illigal.uiuc.edu/web/source-code/1999/03/24/bayesian-optimization-algorithm-boa-code-in-c/。
- 遺伝的アルゴリズム - その理論と先端的手法では、分子の構造エネルギーの最小値を計算しています。
- 1999年。PBILと並んで論文の被引用数が多い。
- 単純な確率ベクトルではなく、ベイジアンネットワークを元に個体を生成します。Learning Bayesian networks: The combination of knowledge and statistical dataによる、ベイジアンネットワークの学習方法で、個体からベイジアンネットワークを学習します。ベイジアンネットワークの評価には、ベイジアン・ディリクレメトリック(Bayesian Dirichlet metric)を提案。高速化のため、貪欲法を使用しています。
- Estimation of Bayesian Networks Algorithm (EBNA)
- Extended Compact Genetic Algorithm (ECGA)
- Hierarchical Bayesian Optimization Algorithm (hBOA)
- Incremental Bayesian Optimization Algorithm (iBOA)
実数値ベクトル
- Stochastic hill climbing with learning by vectors of normal distributions (SHCLVND)
- http://eprints.kfupm.edu.sa/66958/
- http://mikilab.doshisha.ac.jp/dia/research/person/msano/reports/44/report44.html
- 1997年。HCwL をビットストリングから実数値ベクトルにする。値は、正規分布で、標準偏差を世代ごとに減少させ、収束させる。HCwLより成績悪い。
- Real-coded PBIL
- Probabilistic Incremental Program Evolution (PIPE)
- Estimation of Gaussian Networks Algorithm (EGNA)
書籍
- 遺伝的アルゴリズム - その理論と先端的手法
- Advances in Evolutionary Algorithms: Theory, Design and Practice (Studies in Computational Intelligence)
- A Survey of Optimization by Building and Using Probabilistic Models
- Scalable Optimization via Probabilistic Modeling: From Algorithms to Applications (Studies in Computational Intelligence)
PBILを実装した感想
PBILを実装してみた感想。主要部分は50行程度。論文通り、シンプルGAよりもいい感じ。人工知能・機械学習で定番のノイズが乗った関数の最小値にはいい感じ。
試してみた関数。どれも、特に問題なく最小値を出します。計算時間は一瞬(1秒未満)。個体数100、世代数1000、Negative Learning Rate は 0.075。それ以外はPBILの論文通り。実数は普通に2進数で表現。グレイコードは使いませんでした。(その後、グレイコード版を作ってみました。計算時間は、ほぼ同じです。結果は、どっこいどっこい。)
- x^2 + y^2 + 2。xとyの定義域は-1 〜 1、解像度は2^24。
- x^3 - x。xの定義域は-1.1 〜 1.1と-1.2 〜 1.2、解像度は2^24。
- 0.993851231 + exp(-0.01 * x^2) * sin(10 * x) * cos(8 * x)。xの定義域は-3 〜 3、解像度は2^24。HCwLに載っていた関数。
- 0.993851231 + exp(-0.001 * z^2) * sin(10 * z) * cos(8 * z)。z = x + 9。xの定義域は-10 〜 10、解像度は2^24。SHCLVNDに載っていた関数。これのグラフは以下の通り。
主要部分のソースコードはこんな感じ。
/** Population-Based Incremental Learning (PBIL) */ public class PBIL { private final int N, ITER_COUNT, totalBits; private final double[][] domains; private final double learnRate, negLearnRate, mutProb, mutShift; private final Random rand = new Random(); private final DoubleCostFunction costFunc; private final double[] bestCosts; private boolean[] bestGene; private double bestCost; public PBIL(int N, int ITER_COUNT, double[][] domains, DoubleCostFunction costFunc) { this(N, ITER_COUNT, domains, costFunc, 0.1, 0.075, 0.02, 0.05); } public PBIL(int N, int ITER_COUNT, double[][] domains, DoubleCostFunction costFunc, double learnRate, double negLearnRate, double mutProb, double mutShift) { this.N = N; this.ITER_COUNT = ITER_COUNT; this.domains = domains; this.costFunc = costFunc; this.learnRate = learnRate; this.negLearnRate = negLearnRate; this.mutProb = mutProb; this.mutShift = mutShift; this.totalBits = getTotalBits(domains); this.bestCosts = new double[ITER_COUNT]; } public void optimize() { bestCost = POSITIVE_INFINITY; final double[] probVec = new double[totalBits]; Arrays.fill(probVec, 0.5); for (int i = 0; i < ITER_COUNT; i++) { // N個の個体生成 final boolean[][] genes = new boolean[N][totalBits]; for (boolean[] gene : genes) { for (int k = 0; k < gene.length; k++) { if (rand.nextDouble() < probVec[k]) gene[k] = true; } } // コスト計算 TODO 並列化 final double[] costs = new double[N]; for (int j = 0; j < N; j++) { costs[j] = costFunc.cost(toRealVec(genes[j], domains)); } // コスト最小・最大の個体を探す boolean[] minGene = null, maxGene = null; double minCost = POSITIVE_INFINITY, maxCost = NEGATIVE_INFINITY; for (int j = 0; j < N; j++) { double cost = costs[j]; if (minCost > cost) { minCost = cost; minGene = genes[j]; } if (maxCost < cost) { maxCost = cost; maxGene = genes[j]; } } // 全体の最小と比較 if (bestCost > minCost) { bestCost = minCost; bestGene = minGene; } bestCosts[i] = bestCost; // 最小・最大コストで確率ベクトルを更新 for (int j = 0; j < totalBits; j++) { if (minGene[j] == maxGene[j]) { probVec[j] = probVec[j] * (1d - learnRate) + (minGene[j] ? 1d : 0d) * learnRate; } else { final double learnRate2 = learnRate + negLearnRate; probVec[j] = probVec[j] * (1d - learnRate2) + (minGene[j] ? 1d : 0d) * learnRate2; } } // 突然変異 for (int j = 0; j < totalBits; j++) { if (rand.nextDouble() < mutProb) { probVec[j] = probVec[j] * (1d - mutShift) + (rand.nextBoolean() ? 1d : 0d) * mutShift; } } } } public boolean[] getBestBits() { return bestGene; } public double[] getBestRealVec() { return toRealVec(bestGene, domains); } public double getBestCost() { return bestCost; } public double[] getBestCosts() { return bestCosts; } }
UMDAのアルゴリズム
- を長さnの確率ベクトルとし、に初期化。
- 以下、所定の回数ループを繰り返す。
- の確率ベクトルに基づき、長さnのビット列をN個ランダムに生成します。
- 評価関数にて、ビット列を評価します。評価値は実数。
- ビット列を評価値でソートする。全ループを通じて最良だったビット列を保存します。
- 上位M個のビット列の各ビットごとの1の割合をとする。MはN/2など。
- にて更新する。は学習率、0.1など。
- 2に戻る。
評価関数の引数を実数ベクトルにしたい場合は、定義域を決め、ビット列を実数に変換すればよいです。
UMDAのソースコード
/** * Univariate Marginal Distribution Algorithm。 * 単変量周辺分布アルゴリズム。 */ public class UMDA { private final int M, N, LOOP_COUNT, totalBits; private final double learningRate; private final double[][] domains; private final DoubleCostFunction costFunc; private final Random rand = new Random(); private boolean[] bestBits; private double bestCost; public UMDA(int N, int LOOP_COUNT, double[][] domains, DoubleCostFunction costFunc) { this(N, LOOP_COUNT, domains, costFunc, 0.1, N / 2); } public UMDA(int N, int LOOP_COUNT, double[][] domains, DoubleCostFunction costFunc, double learningRate, int M) { this.N = N; this.LOOP_COUNT = LOOP_COUNT; this.domains = domains; this.costFunc = costFunc; this.learningRate = learningRate; this.M = M; this.totalBits = getTotalBits(domains); } public void optimize() { bestCost = Double.MAX_VALUE; final double[] uniFreq = new double[totalBits]; Arrays.fill(uniFreq, 0.5); Gene[] genes = new Gene[N]; for (int iter = 0; iter < LOOP_COUNT; iter++) { // 新しい個体作成とコスト計算 TODO 並列化 for (int j = 0; j < N; j++) { genes[j] = createNewGene(uniFreq); genes[j].cost = costFunc.cost(toRealVec(genes[j].bits, domains)); } // コストでソート Arrays.sort(genes); // 全体の最小と比較 if (bestCost > genes[0].cost) { bestCost = genes[0].cost; bestBits = genes[0].bits; } // 単変量頻度の更新 double[] nextUniFreq = calcUniFreq(genes); for (int k = 0; k < totalBits; k++) { uniFreq[k] += learningRate * (nextUniFreq[k] - uniFreq[k]); } } } /** 単変量頻度 */ private double[] calcUniFreq(Gene[] genes) { double[] uniFreq = new double[totalBits]; for (int k = 0; k < totalBits; k++) { int varCount = 0; for (int j = 0; j < M; j++) { if (genes[j].bits[k]) varCount++; } uniFreq[k] = (double) varCount / (double) M; } return uniFreq; } /** 新しい個体作成 */ private Gene createNewGene(double[] uniFreq) { boolean[] bits = new boolean[totalBits]; for (int k = 0; k < totalBits; k++) { if (rand.nextDouble() < uniFreq[k]) bits[k] = true; } return new Gene(bits); } /** 遺伝子=ビット列+コスト */ private static class Gene implements Comparable { final boolean[] bits; double cost; private Gene(boolean[] bits) { this.bits = bits; } public int compareTo(Object o) { Gene g1 = (Gene) o; return (int) signum(cost - g1.cost); } } public boolean[] getBestBits() { return bestBits; } public double getBestCost() { return bestCost; } }
MIMIC について
UMDAまであった、学習率がなくなっています。その代わり、一度に全ての個体を入れ替えるのではなく、一部だけ(ごく少数だけ)を入れ替えるようになっています。これで、学習率と同じ効果を生んでいます。例えば、1度に下位10%を交換し、上位90%から確率分布を生成します。なので、UMDAと同一の評価回数にするには、反復回数をそれに合わせて、この場合は10倍に増やす必要があります。
一番大事な部分は論文の4ページ目の(2)のこの式です。この確率に基づいて個体を生成します。
そして、上記の順列は、下記条件付きエントロピーの総和を最小化することで決めます。ただし、計算量がになるように、後ろの項から貪欲法で決めます。貪欲法ではなく、ちゃんと計算すると、有向グラフの全点対間最短路となり、のはず。
このエントロピーの定義は、Numerical Recipes (ニューメリカルレシピ・イン・シー 日本語版―C言語による数値計算のレシピ)日本語版第1版のp.468に書いてあります。
MIMIC のソースコード
/** Mutual Information Maximizing Input Clustering (MIMIC) */ public class MIMIC { private final int N, LOOP_COUNT, totalBits, truncateSize; private final double[][] domains; private final DoubleCostFunction costFunc; private final Random rand = new Random(); private boolean[] bestBits; private double bestCost; public MIMIC(int N, int LOOP_COUNT, double[][] domains, DoubleCostFunction costFunc) { this(N, LOOP_COUNT, domains, costFunc, 0.1); } public MIMIC(int N, int LOOP_COUNT, double[][] domains, DoubleCostFunction costFunc, double learnRate) { this.N = N; this.LOOP_COUNT = LOOP_COUNT; this.domains = domains; this.costFunc = costFunc; this.truncateSize = (int)(N * learnRate); this.totalBits = getTotalBits(domains); } public void optimize() { bestCost = Double.MAX_VALUE; // N個の個体をランダムに作成 Gene[] genes = createRandomGenes(); for (int iter = 0; iter < LOOP_COUNT; iter++) { // コスト計算 TODO 並列化 for (int j = (iter == 0 ? 0 : N - truncateSize); j < N; j++) { genes[j].cost = costFunc.cost(toRealVec(genes[j].bits, domains)); } // コストでソート Arrays.sort(genes); // 全体の最小と比較 if (bestCost > genes[0].cost) { bestCost = genes[0].cost; bestBits = genes[0].bits; } // 単変量頻度 double[] uniFreq = calcUniFreq(genes); // 二変量頻度 double[][][] biFreq = calcBiFreq(genes); // ビットの順列を決める int[] order = calcOrder(uniFreq, biFreq); // 新しい個体作成 for (int j = N - truncateSize; j < N; j++) { genes[j] = createNewGene(order, uniFreq, biFreq); } } } /** N個の個体をランダムに作成 */ private Gene[] createRandomGenes() { Gene[] genes = new Gene[N]; for (int j = 0; j < N; j++) { boolean[] bits = new boolean[totalBits]; for (int k = 0; k < bits.length; k++) { bits[k] = rand.nextBoolean(); } genes[j] = new Gene(bits); } return genes; } /** 単変量頻度 */ private double[] calcUniFreq(Gene[] genes) { double[] uniFreq = new double[totalBits]; for (int k = 0; k < totalBits; k++) { int varCount = 0; for (int j = 0; j < N - truncateSize; j++) { if (genes[j].bits[k]) varCount++; } uniFreq[k] = (double) varCount / (double) (N - truncateSize); } return uniFreq; } /** 二変量頻度 */ private double[][][] calcBiFreq(Gene[] genes) { double[][][] biFreq = new double[totalBits][totalBits][4]; for (int m = 0; m < totalBits; m++) { for (int n = 0; n < totalBits; n++) { int[] varCount = new int[4]; for (int j = 0; j < N - truncateSize; j++) { int x = genes[j].bits[m] ? 2 : 0; int y = genes[j].bits[n] ? 1 : 0; varCount[x + y]++; } for (int j = 0; j < 4; j++) { biFreq[m][n][j] = (double) varCount[j] / (double) (N - truncateSize); } } } return biFreq; } /** ビットの決定順を決める */ private int[] calcOrder(double[] uniFreq, double[][][] biFreq) { int[] order = new int[totalBits]; boolean[] used = new boolean[totalBits]; // 最初の1ビット目を決める double minEntropy = Double.POSITIVE_INFINITY; for (int k = 0; k < totalBits; k++) { double p = uniFreq[k]; double entropy = -p * safeLog(p) - (1 - p) * safeLog(1 - p); if (entropy < minEntropy) { minEntropy = entropy; order[0] = k; } } int curBit = order[0]; used[curBit] = true; // 2ビット目以降を決める for (int m = 1; m < totalBits; m++) { double p1 = uniFreq[curBit]; minEntropy = Double.POSITIVE_INFINITY; for (int n = 0; n < totalBits; n++) { if (used[n]) continue; double entropy = 0; for (int i = 0; i < 4; i++) { double p2 = biFreq[curBit][n][i]; entropy += -p2 * (safeLog(p2) - safeLog(i >= 2 ? p1 : 1 - p1)); } if (entropy < minEntropy) { minEntropy = entropy; order[m] = n; } } curBit = order[m]; used[curBit] = true; } return order; } /** 新しい個体作成 */ private Gene createNewGene(int[] order, double[] uniFreq, double[][][] biFreq) { boolean[] bits = new boolean[totalBits]; // 最初のビットの中身を決定 int curBit = order[0]; if (rand.nextDouble() < uniFreq[curBit]) bits[curBit] = true; // 2ビット目以降 for (int k = 1; k < totalBits; k++) { double p = biFreq[curBit][order[k]][(bits[curBit] ? 2 : 0) + 1] / (bits[curBit] ? uniFreq[curBit] : 1 - uniFreq[curBit]); if (rand.nextDouble() < p) bits[order[k]] = true; curBit = order[k]; } return new Gene(bits); } /** 遺伝子=ビット列+コスト */ private static class Gene implements Comparable { final boolean[] bits; double cost; private Gene(boolean[] bits) { this.bits = bits; } public int compareTo(Object o) { Gene g1 = (Gene) o; return (int) signum(cost - g1.cost); } } private static double safeLog(double v) { return v <= 0 ? 0 : log(v); } public boolean[] getBestBits() { return bestBits; } public double getBestCost() { return bestCost; } }
比較
PBIL, UMDA, MIMIC の比較。個人的には色々やったのですが、どのケースでもPBILがベストでした。のフーリエ変換を紹介します。xの定義域は、-1〜1で、sinとcosを10個または40個でフーリエ変換します。
sinとcosが合計10個の時の、フーリエ変換した結果の誤差。縦軸が誤差で、横軸が反復回数です。なるべく速く0に近づいた方がよいです。MIMICはめちゃくちゃです。UMDAは学習率0.1です。
フーリエ変換した結果のグラフです。PBILとUMDAは同じ結果に収束しているので、グラフが重なっています。赤い線がフーリエ変換すべき対象のグラフです。
sinとcosが合計40個の時の、フーリエ変換した結果の誤差。計算時間が長すぎるのでMIMICは割愛しました。UMDAは学習率0.1, 0.2, 0.4の3つを試してみました。学習率を上げると速く収束しますが、おかしな結果に収束しやすくなります。これらの中で、PBILが最も結果が良かったです。
フーリエ変換した結果のグラフです。グラフがグチャグチャになるので、UMDAは最も良い結果に収束した学習率0.1を載せます。PBILの方が圧倒的によいです。
Multiple-Restart Stochastic Hillclimbing
上記の実験は、Multiple-Restart Stochastic Hillclimbing (MRSH) でもやってみました(上記のグラフには載っていません)。しかし、圧倒的に PBIL のほうがいいです。MRSH は深さ優先探索に近い手法ですが、PBIL のような幅優先探索の方が、よい結果にたどり着きます。MRSH のアルゴリズムの詳細は、http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.43.1108 をご覧ください。また、このアルゴリズムは、遺伝的アルゴリズム (バイオインフォマティクス・シリーズ) の p.132 にて、ロイヤルロード関数で、遺伝的アルゴリズムの10倍高速に探索する、ランダム突然変異山登り法および、それの改良したものです。その、ロイヤルロード関数でも、PBIL の方が、MRSH よりも探索が速いです。