(cache) コンストラスティブ・ダイバージェンス法を用いた制限ボルツマンマシン(RBM)の実装

この記事はDeep Learning Advent Calendar 2015の1日目です。

早いもので2015年も最後の月になってしまいました。2015年は「ディープラーニング」という単語が（残念ながらバズワードとして）それなりに広まった一年だったと思います。Deep Learning Advent Calendarがそんな2015年の締めくくりにふさわしいものになりますように。

この記事ではDeep Learningの一種である深層ボルツマンマシン(DBM)の基礎となる制限ボルツマンマシン(RBM)について理論から説明していき、C++による実装を示します。実装を行う関係上どうしても数式を多用しますがなるべく式変形を丁寧に示すつもりです。数式は基本的に

深層学習: Deep Learning

作者: 麻生英樹,安田宗樹,前田新一,岡野原大輔,岡谷貴之,久保陽太郎,ボレガラダヌシカ,人工知能学会,神嶌敏弘
出版社/メーカー: 近代科学社
発売日: 2015/11/05
メディア: 単行本
この商品を含むブログ (1件) を見る

の第2章・第3章に従い、式変形を補足するような形で進めていこうと思います。それなりに長くなりますがお付き合いください。

制限ボルツマンマシンのグラフィカルモデル

機械学習ではデータを観測することでデータの分布に関する何らかのモデルを作り、画像が何を写しているのか当てるなどといった問題を解くことが目標になります。制限ボルツマンマシンのような生成モデルでは、観測されるデータが何らかの確率分布から生成されていると考えます。観測データの確率分布を近似するのがRBMの目標になります。

RBMを説明する前に、より一般的なマルコフ確率場と呼ばれる確率モデルを説明します。マルコフ確率場を図示したグラフは無向グラフになります。ここでノードは0または1の値しか取らないことにします。番目のノードを X_i とし、 X_i が取る値 $x_1, x_2 \cdots$ を集めたベクトルを $\mathbf{X}$ とします。
f:id:kivantium:20151129221405p:plain:w300
このときエネルギー関数と呼ばれる関数を
$\Phi(\mathbf{X}) = -\sum_{i\in \Omega} \phi_i(x_i) - \sum_{\{i, j\}\in \varepsilon}\psi_{ij}(x_i, x_j)$
と定義します。 $\Omega$ はノードの集合、 $\epsilon$ はリンクの集合を表します。また、 $\phi_i$ はノードごとの偏りを表す関数、 $\psi_{ij}$ はノードごとの関連性を表す何らかの関数です。

このエネルギー関数を用いて $\mathbf{X}$ の確率分布を
$p(\mathbf{X}) = \frac{1}{Z}\exp(-\Phi(\mathbf{X}))$
のように定義します。
ただし、は確率の和が1になるようにするための規格化定数で、
$Z = \sum_{\mathbf{X}}\exp(-\Phi(\mathbf{X}))$
と表されます。
このような分布をボルツマン分布と呼びます。ボルツマン分布はもともと統計物理で使われる確率分布ですが、その理論が機械学習で援用されているようです。

制限ボルツマンマシンはマルコフ確率場のうちで、グラフの結合や関数の形に制限を加えたものになります。
f:id:kivantium:20151129223341p:plain:w500
制限ボルツマンマシンのグラフは上の図のように可視層と隠れ層の二つの層に分かれています。その名の通り可視層の値が観測され、隠れ層の値は観測されません。隠れ層を用意するのはモデルの複雑さを上げてよりよい近似を行えるようにするためです。ここでポイントになるのは、可視層のノードと隠れ層のノードを結ぶリンクのみが存在し、同じ層のノードを結ぶリンクが存在しないことです。これが後で説明する条件付き独立が成立するための重要な条件となります。

制限ボルツマンマシンのエネルギー関数の形は、
$\Phi(\mathbf{V}, \mathbf{H} | \mathbf{\theta}) = -\sum_{i=1}^n b_i(v_i) - \sum_{j=1}^m c_j h_j - \sum_{i=1}^n\sum_{j=1}^m w_{ij}v_i h_j$
となります。 $b_i, c_j, w_{ij} \in \theta$ はそれぞれパラメータです。確率分布の式や、 v_i, h_j が0または1しか取らないことはマルコフ確率場と同様です。

制限ボルツマンマシンの学習でやるべきことは、この形で表せる確率分布の範囲内で最も観測データの分布に近くなるようなパラメータを探すことです。以下その手法について説明していきます。

制限ボルツマンマシンの条件付き独立性

制限ボルツマンマシンを表す無向グラフでは同じ層のノードをつなぐリンクがないと説明しましたが、この制限から生まれる非常に重要な性質が条件付き独立性です。これは、片方の層の値を固定するともう片方の層の確率確率変数が互いに独立になることを表します。式で書くと、
${\displaystyle P(\mathbf{h}|\mathbf{v}, \mathbf{\theta}) = \prod_{j=1}^m P(h_j | \mathbf{v}, \mathbf{\theta}) } \\ {\displaystyle P(\mathbf{v}|\mathbf{h}, \mathbf{\theta}) = \prod_{i=1}^n P(v_i | \mathbf{h}, \mathbf{\theta}) }$
となります。この性質のおかげでサンプリングが簡単に行えるようになります。

グラフィカルモデルを知っている方であれば、片方の層を固定したときに条件付き独立性が成り立つことは固定されていない層同士を結ぶ経路が遮断されていることから直感的に分かると思います。詳細はPRML8.3.1などを参照してください。

では条件付き独立性がどのような式になるかを見ていきます。
$\begin{align} P(\mathbf{h}|\mathbf{v}, \mathbf{\theta}) &= \frac{P(\mathbf{v}, \mathbf{h} | \mathbf{\theta})}{P(\mathbf{v} | \mathbf{\theta})} \\ &= \frac{P(\mathbf{v}, \mathbf{h} | \mathbf{\theta})}{\sum_{\mathbf{h}}P(\mathbf{v}, \mathbf{h} | \mathbf{\theta})} \\ \\ &= \frac{\frac{1}{Z(\mathbf{\theta})} \exp(\sum_{i=1}^n b_i v_i + \sum_{j=1}^m c_j h_j + \sum_{i=1}^n\sum_{j=1}^m w_{ij}v_j h_j)}{\sum_{\mathbf{h}} \frac{1}{Z(\mathbf{\theta})}\exp(\sum_{i=1}^n b_i v_i + \sum_{j=1}^m c_j h_j + \sum_{i=1}^n\sum_{j=1}^m w_{ij}v_j h_j)} \\ \\ &= \frac{\exp(\sum_{i=1}^n b_i v_i) \exp(\sum_{j=1}^m h_j (c_j + \sum_{i=1}^n w_{ij} v_i))}{\sum_{\mathbf{h}} \exp(\sum_{i=1}^n b_i v_i) \exp(\sum_{j=1}^m h_j (c_j + \sum_{i=1}^n w_{ij} v_i))} \end{align}$
ここで表記を簡単にするために
$\lambda_j = c_j + \sum_{i=1}^n w_{ij} v_i$
とおくと上の式はさらに
$\begin{align} \frac{\exp(\sum_{j=1}^m h_j \lambda_j)}{\sum_{\mathbf{h}}\exp(\sum_{j=1}^m h_j \lambda_j)} &= \frac{\prod_{j=1}^m \exp(h_j \lambda_j)}{\sum_{\mathbf{h}} \prod_{j=1}^m \exp(h_j \lambda_j)} \\ &= \frac{\prod_{j=1}^m \exp(h_j \lambda_j)}{\prod_{j=1}^m 1+\exp(\lambda_j)} \\ &= \prod_{j=1}^m \frac{\exp(h_j \lambda_j)}{1+\exp(\lambda_j)} \\ &= \prod_{j=1}^m P(h_j | \mathbf{v}, \mathbf{\theta}) \end{align}$
と変形されます。これで条件付き独立性が成立することと、各確率変数が従う具体的な確率分布が求まりました。
ここで分かりにくいのが1行目から2行目の変形ですが、これは h_j が0または1の値しか取らないことに注目すると、
$\prod_{j=1}^m 1+\exp(\lambda_j) = (\exp(0\cdot\lambda_1) + \exp(1\cdot\lambda_1))(\exp(0\cdot\lambda_2) + \exp(1\cdot\lambda_2)) \cdots$
を展開していった各項が $\mathbf{h}$ の取りうる全てのパターンに対応していることから等しいと分かります。

同様の計算により、 $\lambda_i = b_i + \sum_{j=1}^m w_{ij} h_j$ とすれば、
$P(\mathbf{v}|\mathbf{h}, \mathbf{\theta}) = \prod_{i=1}^n \frac{\exp(v_i \lambda_i)}{1+\exp(\lambda_i)}$
であることが分かります。

こうしてRBMの条件付き独立性が分かったので、確率分布に従ってサンプリングを行うことができるようになりました。
以下では具体的にサンプリングの方法を説明していきます。

隠れ層のサンプリング

数式をたくさん書いて疲れたのでC++のコードを示しながら説明したいと思います。

まず、RBMを表すクラスを作ります。

class RBM {
    const int n_visible; // 可視変数の数
    const int n_hidden;  // 隠れ変数の数
    std::vector<float> b, c; // パラメータb, c
    std::vector<std::vector<float>> w; // パラメータw

    public:
    RBM(int visible_size, int hidden_size);
};

RBM::RBM(int visible_size, int hidden_size)
 : n_visible(visible_size), n_hidden(hidden_size) {
    // bを0で初期化(range-based for文)
    b.resize(n_visible);
    for(auto &ref : b) ref = 0.0;

    // cを0で初期化
    c.resize(n_hidden);
    for(auto &ref : c) ref = 0.0;

    // 重みを平均0, 標準偏差0.01の正規分布に従う乱数で初期化
    std::random_device rd;
    std::mt19937 mt(rd());
    std::normal_distribution<float> rand(0, 0.01);

    w.resize(n_visible);
    for(auto &ref : w) {
        ref.resize(n_hidden);
        for(auto &ref2 : ref) {
            ref2 = rand(mt);
        }
    }
}

重みの初期化についてはHintonのA Practical Guide to Training Restricted Boltzmann Machinesを参考にしました。

次に、可視層が与えられたときに隠れ層が1を取る確率を計算する関数を定義します。これは上の式変形からも分かるようにそれぞれの要素について計算することが出来るので次のようになります。

std::vector<float> RBM::get_probability_hidden(const std::vector<int> &visible) {
    // 入力サイズの確認
    assert(visible.size()==static_cast<unsigned>(n_visible));

    // 確率を入れるためのvectorの準備
    std::vector<float> prob;
    prob.resize(n_hidden);

    // 確率の計算
    for(int j=0; j<n_hidden; ++j) {
        float sum = c[j];
        for(int i=0; i<n_visible; ++i) {
            sum += w[i][j] * visible[i];
        }
        prob[j] = exp(sum)/(1.0+exp(sum));
    }
    return prob;
}

文字はなるべく上の式変形と揃えたので何をやっているか確認してみてください。

この関数を使って得た確率を格納したvectorを次の関数に与えて実際に隠れ層が取る値のサンプルを求めます。

std::vector<int> RBM::get_sample(const std::vector<float> &prob) {
    // サンプルを格納するためのvector
    std::vector<int> sample(prob.size());
    // サンプリングのための乱数の用意
    std::random_device rd;
    std::mt19937 mt(rd());
    std::uniform_real_distribution<float> rand(0, 1.0);

    // サンプルの生成
    for(int i=0; i<prob.size(); ++i) {
        if(rand(mt) <= prob[i]) {
            sample[i] = 1;
        } else {
            sample[i] = 0;
        }
    }
    return sample;
}

この実装だと返り値のvectorが値渡しになるためコピーのコストがかかって遅くなっている点など、まだ高速化の余地がありますが実装が楽になることを優先しました。

このようにして、可視層が与えられた時に隠れ層の確率分布を求めた後で、その確率分布に従うサンプルを得ることができます。
隠れ層が与えられた時に可視層のサンプルを得る方法も同様です。

RBMの学習

RBMの学習の目的は、パラメータを調整してRBMが表す確率分布が入力データを生成している確率分布の近似になるようにすることです。
近似のために入力データを生成する真の分布 $q(\mathbf{v})$ とRBMの表す確率分布 $p(\mathbf{v}|\mathbf{\theta})$ のカルバック・ライブラー・ダイバージェンス(KLダイバージェンス）を小さくすることを考えます。

KLダイバージェンスの表式は次の通りです。
$\begin{align} D_{KL}(q(\mathbf{v})||p(\mathbf{v}|\mathbf{\theta})) &= \sum_{\mathbf{v}} q(\mathbf{v}) \log\frac{q(\mathbf{v})}{p(\mathbf{v}|\mathbf{\theta})} \\ &= \mathbf{E}_{q(\mathbf{v})}[\log q(\mathbf{v})] - \mathbf{E}_{q(\mathbf{v})}[\log p(\mathbf{v}|\mathbf{\theta})] \end{align}$

KLダイバージェンスを小さくするために $\mathbf{\theta}$ で偏微分して勾配 $\Delta \theta$ を求めます。第一項が $\mathbf{\theta}$ に依存しないことから、勾配は
$- \frac{\partial D_{KL}}{\partial \mathbf{\theta}} = \frac{\partial \mathbf{E}_{q(\mathbf{v})}[\log p(\mathbf{v}|\mathbf{\theta})]}{\partial \mathbf{\theta}}$
と表されます。

具体的に確率分布の式を代入すると、
$\begin{align} \frac{\partial \mathbf{E}_{q(\mathbf{v})}[\log p(\mathbf{v}|\mathbf{\theta})]}{\partial \mathbf{\theta}} &= \frac{\partial \mathbf{E}_{q(\mathbf{v})}[\log \sum_{\mathbf{h}}p(\mathbf{v}, \mathbf{h}|\mathbf{\theta})]}{\partial \mathbf{\theta}} \\ &= \frac{\partial \mathbf{E}_{q(\mathbf{v})}[\log \sum_{\mathbf{h}}\exp(-\Phi(\mathbf{v}, \mathbf{h}|\mathbf{\theta}))]}{\partial \mathbf{\theta}} - \frac{\partial Z(\mathbf{\theta})}{\partial \theta} \\ &= -\mathbf{E}_{q(\mathbf{v})}\biggl[\sum_{\mathbf{h}} \frac{\exp(-\Phi(\mathbf{v}, \mathbf{h} | \mathbf{\theta}))}{\sum_{\mathbf{h'}}\exp(-\Phi(\mathbf{v}, \mathbf{h'}| \mathbf{\theta}))} \frac{\partial \Phi(\mathbf{v}, \mathbf{h}| \mathbf{\theta})}{\partial \theta}\biggr] + \sum_{\mathbf{v}, \mathbf{h}} \frac{\exp(-\Phi(\mathbf{v}, \mathbf{h}| \mathbf{\theta}))} {\sum_{\mathbf{v'}, \mathbf{h'}} \exp(-\Phi(\mathbf{v'}, \mathbf{h'}| \mathbf{\theta}))} \frac{\partial \Phi(\mathbf{v}, \mathbf{h} | \mathbf{\theta})}{\partial \mathbf{\theta}} \\ &= -\mathbf{E}_{p(\mathbf{h}|\mathbf{v}, \mathbf{\theta})q(\mathbf{v})} \biggl[ \frac{\partial \Phi(\mathbf{v}, \mathbf{h} | \mathbf{\theta})}{\partial \mathbf{\theta}}\biggr] + \mathbf{E}_{p(\mathbf{v}, \mathbf{h} | \mathbf{\theta})} \biggl[ \frac{\partial \Phi(\mathbf{v}, \mathbf{h} | \mathbf{\theta})} {\partial \mathbf{\theta}} \biggr] \end{align}$
のように計算されます。

ここで第一項は $q(\mathbf{v})$ をサンプル分布に置き換えて $p(\mathbf{h}|\mathbf{v}, \mathbf{\theta})$ の $\mathbf{v}$ をサンプルから計算することで近似可能ですが、第二項の期待値を計算するためには $2^{m+n}$ 回の足し算が必要な規格化定数を求める必要があり、実用上計算することが出来ません。

この計算量爆発の問題を解決するために考えられたのが、コンストラスティブ・ダイバージェンス法（CD法）です。

コンストラスティブ・ダイバージェンス法（CD法）

RBMが表す確率分布を真の分布に近づけるためにパラメータの勾配を求めようとしてきましたが、期待値計算で膨大な計算が必要になってしまうことが分かりました。
この期待値を大胆に近似して大幅に計算量を減らす手法がコンストラスティブ・ダイバージェンス法（CD法）です。

CD法ではサンプリングをk回行い、その期待値で先ほど求めた勾配の第2項を近似します。多く使われるのはサンプリングを1回のみ行う手法です。このようなざっくりした近似であってもうまくいくことが経験上知られています。（コンストラスティブ・ダイバージェンスという別の損失関数を計算することで正当化がされていますが、最終的に「経験的に小さいことが分かっている項」を無視して近似しているので経験則と言ってもそんなに間違っていないかと思います。）

具体的なアルゴリズムはBengioのLearning Deep Architectures for AIの60ページに分かりやすく記述されているのでそちらを参照してください。

このアルゴリズムに従ってRBMの訓練を行うコードは次のようになります。

void RBM::train(const std::vector<int> &visible0, const float rate) {
    // サンプルの計算
    auto prob0 = get_probability_hidden(visible0);
    auto hidden0 = get_sample(prob0);
    auto visible1 = get_sample(get_probability_visible(hidden0));
    auto prob1 = get_probability_hidden(visible1);

    // wの更新
    for(int i=0; i<n_visible; ++i) {
        for(int j=0; j<n_hidden; ++j) {
            w[i][j] += rate * (visible0[i]*hidden0[j] - visible1[i]*prob1[j]);
        }
    }
    // bの更新
    for(int i=0; i<n_visible; ++i) {
        b[i] += rate * (visible0[i] - visible1[i]);
    }
    // cの更新
    for(int j=0; j<n_hidden; ++j) {
        c[j] += rate * (hidden0[j] - prob1[j]);
    }
}

以上で理論の説明はおしまいです。（長く苦しい戦いだった……）

ここまでで作成したRBMクラスの全コードを以下に示します。RBM.hという名前で保存してください。

#include <vector>
#include <random>
#include <cassert>
#include <cmath>
class RBM {
    const int n_visible; // 可視変数の数
    const int n_hidden;  // 隠れ変数の数
    std::vector<float> b, c; // パラメータb, c
    std::vector<std::vector<float>> w; // パラメータw
    std::vector<float> get_probability_hidden(const std::vector<int> &visible);
    std::vector<float> get_probability_visible(const std::vector<int> &hidden);
    std::vector<int> get_sample(const std::vector<float> &prob);

    public:
    RBM(int visible_size, int hidden_size);
    void train(const std::vector<int> &visible0, const float rate=0.01);
    std::vector<int> get_reconstruction(const std::vector<int> &visible);
};

RBM::RBM(int visible_size, int hidden_size) : n_visible(visible_size), n_hidden(hidden_size) {
    // bを0で初期化(range-based for文)
    b.resize(n_visible);
    for(auto &ref : b) ref = 0.0;

    // cを0で初期化
    c.resize(n_hidden);
    for(auto &ref : c) ref = 0.0;

    // 重みを平均0, 標準偏差0.01の正規分布に従う乱数で初期化
    std::random_device rd;
    std::mt19937 mt(rd());
    std::normal_distribution<float> rand(0, 0.01);

    w.resize(n_visible);
    for(auto &ref : w) {
        ref.resize(n_hidden);
        for(auto &ref2 : ref) {
            ref2 = rand(mt);
        }
    }
}

std::vector<float> RBM::get_probability_hidden(const std::vector<int> &visible) {
    // check input
    assert(visible.size()==static_cast<unsigned>(n_visible));

    // prepare variables
    std::vector<float> prob;
    prob.resize(n_hidden);

    // caluculate probability
    for(int j=0; j<n_hidden; ++j) {
        float sum = c[j];
        for(int i=0; i<n_visible; ++i) {
            sum += w[i][j] * visible[i];
        }
        prob[j] = exp(sum)/(1.0+exp(sum));
    }
    return prob;
}
std::vector<float> RBM::get_probability_visible(const std::vector<int> &hidden) {
    // check input
    assert(hidden.size()==static_cast<unsigned>(n_hidden));

    // prepare variables
    std::vector<float> prob;
    prob.resize(n_visible);

    // caluculate probability
    for(int i=0; i<n_visible; ++i) {
        float sum = b[i];
        for(int j=0; j<n_hidden; ++j) {
            sum += w[i][j] * hidden[j];
        }
        prob[i] = sum;
    }
    return prob;
}
std::vector<int> RBM::get_sample(const std::vector<float> &prob) {
    // prepare variables
    std::vector<int> sample(prob.size());
    std::random_device rd;
    std::mt19937 mt(rd());
    std::uniform_real_distribution<float> rand(0, 1.0);

    // create sample
    for(int i=0; i<prob.size(); ++i) {
        if(rand(mt) <= prob[i]) {
            sample[i] = 1;
        } else {
            sample[i] = 0;
        }
    }
    return sample;
}

void RBM::train(const std::vector<int> &visible0, const float rate) {
    // サンプルの計算
    auto prob0 = get_probability_hidden(visible0);
    auto hidden0 = get_sample(prob0);
    auto visible1 = get_sample(get_probability_visible(hidden0));
    auto prob1 = get_probability_hidden(visible1);

    // wの更新
    for(int i=0; i<n_visible; ++i) {
        for(int j=0; j<n_hidden; ++j) {
            w[i][j] += rate * (visible0[i]*hidden0[j] - visible1[i]*prob1[j]);
        }
    }
    // bの更新
    for(int i=0; i<n_visible; ++i) {
        b[i] += rate * (visible0[i] - visible1[i]);
    }
    // cの更新
    for(int j=0; j<n_hidden; ++j) {
        c[j] += rate * (hidden0[j] - prob1[j]);
    }
}

std::vector<int> RBM::get_reconstruction(const std::vector<int> &visible) {
    auto hidden = get_sample(get_probability_hidden(visible));
    return get_sample(get_probability_visible(hidden));
}

RBMの使い道

こうして実装することができたRBMですが、どのような使い道があるのか見ておきます。

生成モデルとして使う: RBMは入力の分布を近似していると考えられるので、サンプリングしてデータを作るなどの使い方ができます
教師あり学習に使う: 観測データにそのデータが属するラベルを追加して学習することで、未知のデータに対して各クラスに分類される確率を求めることができるようになる
自己符号化器として使う: RBMの隠れ層は可視変数の特徴を抽出したものと考えられるので次元圧縮の手法としても使えます。
深層ボルツマンマシンの事前学習に使う: 深層ボルツマンマシンをRBMがたくさん重なったものとみて層ごとに学習する

最後の項目について補足しておきます。深層ボルツマンマシンというRBMの隠れ層をさらに積んだようなDeep Learningのモデルがあります。これをRBMがたくさん重なったものとみなして順番に学習していくことで学習を効率的に進めることができます。RBM自体はDeep Learningの手法ではないのですが、Deep Learningを支える技術の一つなので少し無理やりですがDeep Learning Advent Calendarの記事とさせてもらいました。深層ボルツマンマシンについてはそのうち説明できたらいいな……

RBMを用いたデータ再構成

可視層にデータを与えて隠れ層のサンプリングを行った後、隠れ層のデータから再び可視層のデータを復元できることができます。うまく学習ができていれば隠れ層が可視層の特徴をうまく抽出してくれるので復元したデータは最初に与えたデータと近いものになるはずです。そこで先ほど作ったRBMクラスを使って入力が再構成できるか実験してみます。

#include <iostream>
#include "RBM.h"
using namespace std;

int main() {
    // 可視層20, 隠れ層15のRBMを作成
    RBM myRBM(20, 15);

    // 入力データ
    vector<int> v0{1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    vector<int> v1{0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
    vector<int> v2{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0};
    vector<int> v3{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1};

    // 学習率を変えて学習するのが効果的
    for(int i=0; i<100; ++i) {
        myRBM.train(v0, 0.1);
        myRBM.train(v1, 0.1);
        myRBM.train(v2, 0.1);
        myRBM.train(v3, 0.1);
    }
    for(int i=0; i<1000; ++i) {
        myRBM.train(v0, 0.001);
        myRBM.train(v1, 0.001);
        myRBM.train(v2, 0.001);
        myRBM.train(v3, 0.001);
    }

    // それぞれのデータについて5回サンプリング
    for(int i=0; i<5; ++i) {
        auto sample = myRBM.get_reconstruction(v0);
        for(auto &ref : sample) cout << ref << " ";
        cout << endl;
    }
    for(int i=0; i<5; ++i) {
        auto sample = myRBM.get_reconstruction(v1);
        for(auto &ref : sample) cout << ref << " ";
        cout << endl;
    }
    for(int i=0; i<5; ++i) {
        auto sample = myRBM.get_reconstruction(v2);
        for(auto &ref : sample) cout << ref << " ";
        cout << endl;
    }
    for(int i=0; i<5; ++i) {
        auto sample = myRBM.get_reconstruction(v3);
        for(auto &ref : sample) cout << ref << " ";
        cout << endl;
    }
}

出力は

1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 0 0 0 0 0 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1 
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 1 1 1 1

となり、ちゃんと入力を復元できていることが分かりました。
学習回数をうまく設定しないと復元がうまく行かないのでそれなりにパラメータに依存しているようです。

これだけだと面白くないので画像を入力して復元できるかどうかも実験したいと思います。

#include <opencv2/opencv.hpp>
#include <iostream>
#include <vector>
#include "RBM.h"
 
using namespace cv;
using namespace std;
 
int main(int argc, char **argv) {
    // 画像の読み込み
    const int size = 50;
    Mat src_img = imread(argv[1], 0);
    // 画像の縮小
    Mat bin_img;
    resize(src_img, bin_img, Size(size, size));
    // 画像の二値化
    threshold(bin_img, bin_img, 0, 255, THRESH_BINARY|THRESH_OTSU);
    // 元画像の保存
    imwrite("original.png", bin_img);

    // データフォーマットの変換(Mat->vector)
    vector<int> image(size*size);
    for(int y=0; y<bin_img.rows; y++){
        for(int x=0; x<bin_img.cols; x++){
            if ((int) bin_img.data[y*bin_img.cols + x] == 0) image[y*bin_img.cols + x] = 0;
            else image[y*bin_img.cols + x] = 1;
        }
    }

    // RBMの訓練
    RBM myRBM(size*size, size*size/2);
    // 学習率を変えて2段階でやるとうまくいった(経験則)
    for(int i=0; i<100; ++i) {
        myRBM.train(image, 0.1);
    }
    for(int i=0; i<1000; ++i) {
        myRBM.train(image, 0.001);
    }

    // 復元画像の生成
    auto sample = myRBM.get_reconstruction(image);
    for(int y=0; y<bin_img.rows; y++){
        for(int x=0; x<bin_img.cols; x++){
            if (sample[y*bin_img.cols + x] == 0) bin_img.data[y*bin_img.cols + x] = 0;
            else bin_img.data[y*bin_img.cols + x] = 255;
        }
    }

    // 復元画像の保存
    imwrite("reconstruct.png", bin_img);
}

入力にはこのAdvent Calendarの主催者の@olanleedさんのアイコンを使いました。二値化して50x50にした後で復元でいるかどうかを試します。
元画像とフォーマット変換後の画像は以下の通りです。
f:id:kivantium:20151130234259p:plain:w300 f:id:kivantium:20151130234312p:plain:w300