Gibbs sampling について - 情報理論関連をぐだぐだと

ボルツマンマシンの状態数

前回はなしたボルツマンマシンだけれど、実際に学習を実行しようとなると、あるパラメータにおける期待値を計算する必要がある。期待値は例えば、

${ \displaystyle E[X_i X_j] = \sum_{x\in\mathcal{X}}p(x)x_ix_j }$

というようなもので、ボルツマンマシンは $\{ 0,1\}^n$ 上の確率分布だったから、 $n$ が大きくなると足す必要のある数は爆発的に大きくなる。要するに、実際の計算機ではやっていられなくなる。そこでちょっとした工夫が必要になる。

全結合ボルツマンマシンの例

今回は次のボルツマンマシンを使う。 $n = 4$

${ \displaystyle b = ( 1, -3.5, -1, 2) }$

$w =$


0	-2	-1	3
-2	0	-1	2
-1	-1	0	2
3	2	2	0

このボルツマンマシンの分布は次のようになる。

状態	確率（％）
0000	0.06
0001	0.16
0010	0.00
0011	0.03
0100	0.02
0101	0.02
0110	0.00
0111	0.00
1000	0.42
1001	23.02
1010	0.09
1011	37.95
1100	1.15
1101	23.02
1110	0.09
1111	13.96

ちなみにこの分布を求めて、そこからサンプリングして経験分布を得るC++14のコードは次のようになる*1。

#include <iostream>
#include <bitset>
#include <cmath>
#include <iomanip>
#include <random>
using namespace std;

int main() {
    random_device rd;
    mt19937 mt(rd());
    uniform_real_distribution<double> score(0.0,1.0);
    
    double energy[16] = {0};
    double e_ene[16] = {0};
    double prob[16];
    double i_prob[16] = {0};
    double t_e_ene = 0;
    double b[4] = {1,-3.5,-1,2};
    double w[4][4] = {
        { 0, 2,-1, 3},
        { 2, 0,-1, 2},
        {-1,-1, 0, 2},
        { 3, 2, 2, 0}
    };
    
    for(int i = 0; i < 16; ++i){
        bitset<4> sample(i);
        for(int j = 0; j < 4; ++j){
            energy[i] -= sample[j] * b[j];
            for(int k = j + 1; k < 4; ++k){
                energy[i] -= sample[j] * sample[k] *w[j][k];
            }
        }
        e_ene[i] = exp( -energy[i] );
        t_e_ene += e_ene[i];
    }
    for(int i = 0; i < 16; ++i){
        prob[i] = e_ene[i] / t_e_ene;
        if( i ) {i_prob[i] = i_prob[i-1] +  prob[i];}
        else{i_prob[i] =  prob[i];}
    }
    
    int n_count[16] = {0};
    int sample_num = 1000000;
    for(int j = 0; j < sample_num; ++j){
        double sc = score(mt);
        for(int i = 0; i < 16; ++i){
            if(i_prob[i] > sc){
                n_count[i]++ ;
                break;
            }
        }
    }
    
    for(int i = 0; i < 16; ++i){
        cout << static_cast< bitset<4> >(i) << fixed << setprecision(2)<<" " << (double) 100 * n_count[i] / sample_num <<" " << prob[i] * 100 <<  endl;
    }
    return 0;
}

Gibbs sampling

このサンプリング法では、すべての状態についての確率を求めていた*2。しかし、 $n$ が増えれば増えるほど、このやり方は無理になっていく。

ここで、ボルツマンマシンの構造からギブスサンプリングを用いることを考える。

ギブスサンプリングとは、ある状態のあるユニットに注目し、他のユニットの値は定まっているものとして、その注目しているユニットの値を更新していく方法で、十分長い間この更新を繰り返すと所望の分布から生成されたサンプルが得られると理論的には保証されている方法だ。

あるユニット $x_i$ に注目したときに、他のユニットの値が定まっているときのユニット $x_i$ の分布は、

${ \displaystyle {\rm Pr}(X_i = 0|b,w,\{x_j\}_{j \neq i}) = \frac{1}{ 1 + {\rm exp}(b_i + \sum_{j \neq i}w_{i,j}x_j)} }$

で定まる。

この分布を用いてサンプルを行い、もとの分布と比較するコードを次に紹介する。

#include <iostream>
#include <bitset>
#include <cmath>
#include <iomanip>
#include <random>
using namespace std;

int main() {
    random_device rd;
    mt19937 mt(rd());
    uniform_real_distribution<double> score(0.0,1.0);
    
    double energy[16] = {0};
    double e_ene[16] = {0};
    double prob[16];
    double i_prob[16] = {0};
    double t_e_ene = 0;
    double b[4] = {1,-3.5,-1,2};
    double w[4][4] = {
        { 0, 2,-1, 3},
        { 2, 0,-1, 2},
        {-1,-1, 0, 2},
        { 3, 2, 2, 0}
    };
    double p_ene[4];
    bitset<4> g_sample (0);

    int g_n_count[16] = {0};
    int g_sample_num = 1000000;
    int g_s_t = 0;
    for(int i = 0; i <  g_sample_num; ++i){
        for(int j = 0; j < 4; ++j){
            p_ene[j] = b[j];
            for(int k = 0; k < 4; ++k){
                if(k != j)    p_ene[j] += w[k][j] * g_sample[k];
            }
        }
        int i_site = i % 4;
        double p = 1.0 / ( 1.0 + exp(p_ene[i_site]));
        double sc = score(mt);
        g_sample[i_site] = (sc > p);
        if((i > g_sample_num / 10) && (i % 51 == 0) ){
            int nnn = g_sample.to_ulong();
            g_n_count[nnn]++;
            g_s_t++;
        }
    }
    
    for(int i = 0; i < 16; ++i){
        bitset<4> sample(i);
        for(int j = 0; j < 4; ++j){
            energy[i] -= sample[j] * b[j];
            for(int k = j + 1; k < 4; ++k){
                energy[i] -= sample[j] * sample[k] *w[j][k];
            }
        }
        e_ene[i] = exp( -energy[i] );
        t_e_ene += e_ene[i];
    }
    for(int i = 0; i < 16; ++i){
        prob[i] = e_ene[i] / t_e_ene;
        if( i ) {i_prob[i] = i_prob[i-1] +  prob[i];}
        else{i_prob[i] =  prob[i];}
    }
    for(int i = 0; i < 16; ++i){
        cout << static_cast< bitset<4> >(i) << fixed << setprecision(2)<<" " << (double) 100 * g_n_count[i] / g_s_t <<" " << prob[i] * 100 <<  endl;
    }
    return 0;
}

このコードでは状態を0000から始めて、更新回数の10分の1をburn-in期間において51回の更新毎にサンプリングしてある。結果の一例を次に書く*3。

状態	Gibbs samplingによる経験分布の確率（％）	計算による真の確率（％）
0000	0.06	0.06
0001	0.17	0.16
0010	0.00	0.00
0011	0.02	0.03
0100	0.03	0.02
0101	0.02	0.02
0110	0.00	0.00
0111	0.00	0.00
1000	0.37	0.42
1001	23.06	23.02
1010	0.14	0.09
1011	37.85	37.95
1100	1.14	1.15
1101	23.43	23.02
1110	0.07	0.09
1111	13.65	13.96