Hatena::ブログ(Diary)

yukobaのブログ このページをアンテナに追加 RSSフィード Twitter

2009年04月26日 分布推定アルゴリズム

分布推定アルゴリズム

分布推定アルゴリズム遺伝的アルゴリズムを改良した物です。個体の集合を交叉・突然変異させるのではなく、個体の生成確率進化させます。最適化問題アルゴリズムです。以下、自分へのメモです。わかったことが増えたら追記するかも。

ビットストリング

計算量に関しては、ビット数をn、反復数をTとしています。

実数値ベクトル

書籍

ソースコード

以下順に、以下の3つのソースコードがあります。(最初の頃のソースコードがバグってました。すいません。)

  1. PBIL
  2. UMDA
  3. MIMIC

この3つの中では、PBILが最も良いと思います。比較は一番下です。

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に載っていた関数。これのグラフは以下の通り。

f:id:yukoba:20090427021625p:image

主要部分のソースコードはこんな感じ。

/** 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のアルゴリズム

  1. Pを長さnの確率ベクトルとし、P_i=0.5初期化
  2. 以下、所定の回数ループを繰り返す。
  3. Pの確率ベクトルに基づき、長さnのビット列をN個ランダムに生成します。
  4. 評価関数f()にて、ビット列を評価します。評価値は実数。
  5. ビット列を評価値でソートする。全ループを通じて最良だったビット列を保存します。
  6. 上位M個のビット列の各ビットごとの1の割合をP’とする。MはN/2など。
  7. P_i = P_i + ¥alpha (P’_i - P_i) にて更新する。¥alphaは学習率、0.1など。
  8. 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)のこの式です。この確率に基づいて個体を生成します。

p_{¥pi}(X) = p(X_{i_1}|X_{i_2})p(X_{i_2}|X_{i_3}) ¥dots p(X_{i_{n-1}}|X_{i_n}) p(X_{i_n})

そして、上記の順列は、下記条件付きエントロピーの総和を最小化することで決めます。ただし、計算量がO(n^2)になるように、後ろの項から貪欲法で決めます。貪欲法ではなく、ちゃんと計算すると、有向グラフの全点対間最短路となり、O(n^3)のはず。

h(X_{i_1}|X_{i_2}) + h(X_{i_2}|X_{i_3}) + ¥dots + h(X_{i_{n-1}}|X_{i_n}) + h(X_{i_n})

このエントロピーの定義は、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(x - 0.32)(x + 0.8)(x - 0.66)(x + 0.21)フーリエ変換を紹介します。xの定義域は、-1〜1で、sinとcosを10個または40個でフーリエ変換します。

sinとcosが合計10個の時の、フーリエ変換した結果の誤差。縦軸が誤差で、横軸が反復回数です。なるべく速く0に近づいた方がよいです。MIMICはめちゃくちゃです。UMDAは学習率0.1です。

f:id:yukoba:20090613003024p:image

フーリエ変換した結果のグラフです。PBILとUMDAは同じ結果に収束しているので、グラフが重なっています。赤い線がフーリエ変換すべき対象のグラフです。

f:id:yukoba:20090613003025p:image

sinとcosが合計40個の時の、フーリエ変換した結果の誤差。計算時間が長すぎるのでMIMICは割愛しました。UMDAは学習率0.1, 0.2, 0.4の3つを試してみました。学習率を上げると速く収束しますが、おかしな結果に収束しやすくなります。これらの中で、PBILが最も結果が良かったです。

f:id:yukoba:20090613003028p:image

フーリエ変換した結果のグラフです。グラフがグチャグチャになるので、UMDAは最も良い結果に収束した学習率0.1を載せます。PBILの方が圧倒的によいです。

f:id:yukoba:20090613003027p:image

Multiple-Restart Stochastic Hillclimbing

上記の実験は、Multiple-Restart Stochastic Hillclimbing (MRSH) でもやってみました(上記のグラフには載っていません)。しかし、圧倒的に PBIL のほうがいいです。MRSH は深さ優先探索に近い手法ですが、PBIL のような幅優先探索の方が、よい結果にたどり着きます。MRSH のアルゴリズムの詳細は、CiteSeerX — An Empirical Comparison of Seven Iterative and Evolutionary Function Optimization Heuristics をご覧ください。また、このアルゴリズムは、遺伝的アルゴリズム (バイオインフォマティクス・シリーズ) の p.132 にて、ロイヤルロード関数で、遺伝的アルゴリズムの10倍高速に探索する、ランダム突然変異山登り法および、それの改良したものです。その、ロイヤルロード関数でも、PBIL の方が、MRSH よりも探索が速いです。

トラックバック - http://d.hatena.ne.jp/yukoba/20090426