Science

制限ボルツマンマシン(RBM)の実装と理論【Pythonで実装】

制限ボルツマンマシン(RBM)の実装と理論【Pythonで実装】

 

本記事では、制限ボルツマンマシン(RBM)の理論とPythonによる実装を初学者の方でも理解できるように解説しました。

Pythonによる実装では、RBMの実装でよく使用される手法であるパーシステント・コンストラスティブ・ダイバージェンス法(PCD法)と呼ばれる手法を使用して実装しました。

後半部分では、PythonによるRBMの具体的なコードを載せているので、自分で適当に動かして遊んでみてください。

 

制限ボルツマンマシン(RBM)の定義と学習法

本章では、制限ボルツマンマシンの定義と学習法を説明していきます。

  1. 制限ボルツマンマシンの定義
  2. 制限ボルツマンマシンの学習法
  3. 制限ボルツマンマシンの学習方程式の導出

 

制限ボルツマンマシン(RBM)の定義

 

制限ボルツマンマシン(Restricted Boltzmann Machine ; RBM)は、主にデータの背後の隠れた分布を近似する学習モデルです。

簡単にデータの背後にある確率分布と言っていますが、必ずしも各データが背後にある確率分布に従っている保証はないです。

しかし、生成モデルと呼ばれるタスクでは基本的にデータの背後にはデータの生成過程を支配する背後の確率分布が存在していることを仮定します。

 

そこで、制限ボルツマンマシンは、以下のパラメータパラメータ\(\mathbf{\theta} = (\mathbf{w} \in \mathbb{R}^{N_{v} \times N_{h}}, \mathbf{b} \in \mathbb{R}^{N_{v}}, \mathbf{c} \in \mathbb{R}^{N_{h}} )\)によって特徴付けられる確率分布を定義します。

$$P(\mathbf{v}, \mathbf{h}\mid \mathbf{\theta}) = \frac{1}{Z} e^{- E(\mathbf{v}, \mathbf{h} ; \mathbf{\theta})} $$

ここで、関数\(E(\mathbf{v}, \mathbf{h} ; \mathbf{\theta})\)は一般にエネルギー関数と呼ばれ、RBMの場合以下のように定義されます。

$$E(\mathbf{v}, \mathbf{h} ; \mathbf{\theta}) =~ – \sum_{i=1}^{N_{v}} \sum_{j=1}^{N_{h}} w_{ij} v_{i} h_{j} ~- \sum_{i=1}^{N_{v}} b_{i} v_{i} ~- \sum_{j=1}^{N_{h}} c_{j} h_{j}$$

具体的にネットワークの結合をグラフィカルモデルで可視化したものをいかに示します。

RBMのネットワーク図

ここで、一般に変数\(v_{i}\), \(h_{j}\)はそれぞれ可視変数、隠れ変数と呼ばれ\(v_{i} \in \{0, 1\}\), \(h_{j} \in \{0, 1\}\)の値をとります。

また、\(Z\)は一般に分配関数と呼ばれ以下のように定義されます。

$$Z = \sum_{\mathbf{v}, \mathbf{h}} e^{- E(\mathbf{v}, \mathbf{h} ; \mathbf{\theta})} $$

ここで、\(\sum_{\mathbf{v}, \mathbf{h}} \cdots \)は可視変数と隠れ変数に関して可能な状態全てに関する和を取ることを意味しています。

制限ボルツマンマシンの確率分布をグラフィカルモデルを利用して表したものが以下のようになります。

そして、制限ボルツマンマシンは隠れ変数に関して周辺化した以下の確率分布を学習モデルとしてデータの背後にある確率分布を近似的に実現します。

$$P(\mathbf{v} \mid \mathbf{\theta}) = \sum_{\mathbf{h}} P(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta})$$

 

RBMの条件付き独立性

 

RBMのようなネットワークを考える最大のメリットは、条件付き独立性です。

具体的に可視変数に関する条件付き確率を計算すると以下のように表せます。

$$P(\mathbf{v} \mid \mathbf{h}, \mathbf{\theta}) = \prod_{i} \frac{e^{(b_{i} + \sum_{j} w_{ij} h_{j})v_{i}}}{1 + e^{(b_{i} + \sum_{j} w_{ij} h_{j})}} $$

隠れ変数に関しても以下のように表すことができます。

$$P(\mathbf{h} \mid \mathbf{v}, \mathbf{\theta}) = \prod_{j} \frac{e^{(c_{j} + \sum_{i} w_{ij} v_{i})h_{j}}}{1 + e^{(c_{j} + \sum_{i} w_{ij} v_{i})}} $$

つまり、RBMはグラフの構造上、一方の層を固定したときに他方の層の変数が因数分解の形(独立)で表すことができます。

これが、RBMの条件付き独立性と呼ばれる性質で、学習を効率化させる性質の一つです。

 

RBMの学習方程式

 

データ\(\mathcal{D} = \{\mathbf{v}^{1},\ldots, \mathbf{v}^{P}\}\)が与えられたときに、RBMのパラメータ\(\mathbf{\theta}\)の学習は、一般に以下の対数尤度最大化のよってパラメータを決定します。

$$\mathcal{L}_{\mathcal{D}}(\mathbf{\theta}) \equiv \frac{1}{P} \sum_{\mu=1}^{P} \ln p(\mathbf{v}^{\mu} \mid \mathbf{\theta})$$

上式を最大化するようなパラメータを見つけるためにパラメータに関して微分して\(0\)となる条件を求めると以下のような条件式を導くことができます。

$$\mathbb{E}_{\text{data}} \left[ \frac{\partial E(\mathbf{v}^{\mu}, \mathbf{h} ; \mathbf{\theta})}{\partial \mathbf{\theta}} \right] = \mathbb{E}_{\text{RBM}} \left[ \frac{\partial E(\mathbf{v}^{\mu}, \mathbf{h} ; \mathbf{\theta})}{\partial \mathbf{\theta}} \right]$$

ここで、\(\mathbb{E}_{\text{data}}[\cdots]\)は以下のような平均操作を表します。

$$\mathbb{E}_{\text{data}}[\cdots] = \frac{1}{P} \sum_{\mu} \sum_{\mathbf{h}} p(\mathbf{h} \mid \mathbf{v}^{\mu}, \mathbf{\theta}) \cdots$$

また、\(\mathbb{E}_{\text{RBM}}[\cdots]\)は以下のような平均操作を表します。

$$\mathbb{E}_{\text{RBM}}[\cdots] =  \sum_{\mathbf{v}, \mathbf{h}} p(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta}) \cdots$$

実際に各パラメータに関して計算すると以下のような連立方程式を得ることができます。

\begin{align}&\frac{\partial \mathcal{L}_{\mathcal{D}}(\mathbf{\theta})}{\partial w_{ij}} : \mathbb{E}_{\text{data}}[v_{i} h_{j}] = \mathbb{E}_{\text{RBM}}[v_{i} h_{j}] \\ &\frac{\partial \mathcal{L}_{\mathcal{D}}(\mathbf{\theta})}{\partial b_{i}} : \mathbb{E}_{\text{data}}[v_{i}] = \mathbb{E}_{\text{RBM}}[v_{i}] \\ &\frac{\partial \mathcal{L}_{\mathcal{D}}(\mathbf{\theta})}{\partial c_{j}} : \mathbb{E}_{\text{data}}[h_{j}] = \mathbb{E}_{\text{RBM}}[ h_{j}]\end{align}

この連立方程式を一般に『RBMの学習方程式』と呼びます。

 

RBMの学習法

 

RBMの学習方程式を得られましたが、一般に学習方程式を解析的に解くのは困難なため、数値的に以下の勾配上昇方によってパラメータを逐次更新し学習を実行します。

RBMの学習更新式(勾配上昇法)

\begin{align}&w_{ij}^{(t+1)} \leftarrow w_{ij}^{(t)} + \eta \big( \mathbb{E}_{\text{data}}[v_{i} h_{j}] – \mathbb{E}_{\text{RBM}}[v_{i} h_{j}] \big) \\ &b_{i}^{(t+1)} \leftarrow b_{i}^{(t)} + \eta \big(\mathbb{E}_{\text{data}}[v_{i}] – \mathbb{E}_{\text{RBM}}[v_{i}] \big) \\ &c_{j}^{(t+1)} \leftarrow c_{j}^{(t)} + \eta \big( \mathbb{E}_{\text{data}}[h_{j}] – \mathbb{E}_{\text{RBM}}[ h_{j}] \big) \end{align}

ここで、\(\eta \in \mathbb{R}\)は一般に学習率と呼ばれ小さな値をとります。

この更新式を収束するまで更新し続けます。

また、RBMではデータに関する平均の項を『ポジティブフェーズ』と呼び、RBMの確率分布に関する項を『ネガティブフェーズ』と呼びます。

ここまでの結果をまとめると学習方程式を解析的に解く問題が、各更新でネガティブフェーズとポジティブフェーズを評価する問題になりました。

ポジティブフェーズとネガティブフェーズの評価は簡単?

 

パラメータを学習する問題は、各更新でネガティブフェーズとポジティブフェーズを評価する問題に変わりました。

では、ネガティブフェイズとポジティブフェーズの評価は簡単でしょうか?

まずは、positive phaseの評価に関して考察していきます。

 

Positive phaseの評価について

 

positive phaseは以下の期待値を評価することに相当することを思い出しましょう。

$$\mathbb{E}_{\text{data}}[\cdots] = \frac{1}{P} \sum_{\mu} \sum_{\mathbf{h}} p(\mathbf{h} \mid \mathbf{v}^{\mu}, \mathbf{\theta}) \cdots$$

まず、データに関する標本平均は簡単に取ることができます。

また、隠れ変数に関する条件付き確率に関する平均もRBMの条件付き独立性のため簡単に評価することができます。

あまりにデータ数が多い場合は、mini-batch法などを使用すれば基本的に問題はありません。

 

Negative phaseの評価について

 

negative phaseの以下の期待値を評価することを思い出しましょう。

$$\mathbb{E}_{\text{RBM}}[\cdots] =  \sum_{\mathbf{v}, \mathbf{h}} p(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta}) \cdots$$

つまり、\(p(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta})\)の計算を行う必要があります。

しかし、\(p(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta})\)を評価するためには以下の分配関数を評価する必要があります。

$$Z = \sum_{\mathbf{v}, \mathbf{h}} e^{- E(\mathbf{v}, \mathbf{h} ; \mathbf{\theta})} $$

繰り返しになりますが、\(\sum_{\mathbf{v}, \mathbf{h}} \cdots\)は\(\mathbf{v}\), \(\mathbf{h}\)の可能な状態に関する和を取ることを意味しています。

仮に可視変数と隠れ変数が\(20\)個あったら、状態数は約\(1.0 \times 10^{6}\)となってしまいます。

このような状態和の計算は一般に困難で『計算量爆発』と呼ばれたりします。

そのため、Negative Phaseに関しては近似的学習法により評価することで学習を行います。

次章では、Negative Phaseの近似的学習法に関して説明していきます。

 

制限ボルツマンマシンの近似的学習法

前章では、制限ボルツマンマシンの学習は、Negative Phaseの計算が計算量爆発になり困難であることを説明しました。

本章では、Negative Phaseの計算を近似的に行う手法を紹介していきます。

しかし、Negative Phaseの効率的かつ精度の高い評価方法は現在も進行中の研究のため、現時点で最も使われている近似的な評価方法を紹介していきます。

具体的には以下の流れで説明していきます。

  1. マルコフ連鎖モンテカルロ法による期待値の近似評価
  2. ギブスサンプリングによる評価と問題点
  3. コンストラスティブダイバージェンス法(CD法)
  4. パーシステントコンストラスティブダイバージェンス法(PCD法)

 

マルコフ連鎖モンテカルロ法による期待値評価

 

ネガティブフェイズの期待値評価を厳密に行うことは困難なので、マルコフ連鎖モンテカルロ法を使用し、期待値を近似的に評価することを目指します。

『マルコフ連鎖モンテカルロ法??』という方は下記を参考にしてください。

 

モンテカルロ法の基本をわかりやすく説明してみた|マルコフ連鎖モンテカルロ法まで本記事では、初心者でも理解できるようにモンテカルロ法の基本的な概念からマルコフ連鎖モンテカルロ法までを段階的に説明しました。さらに、適宜Pythonによるシミュレーションを行いながら視覚的に理解できるように努めました。...

ギブスサンプリングによる評価と問題点

 

制限ボルツマンマシンは、マルコフ連鎖モンテカルロ法の中でもギブスサンプリングという手法を利用することが多いです。

具体的には、以下のようなアルゴリズムで\(P(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta})\)からのサンプル列を取得します。

制限ボルツマンマシンのギブスサンプリング

  1. \(\mathbf{v}\)をランダムに初期化し、\(\mathbf{v}^{0}\)と定義する
  2. \(P(\mathbf{h} \mid \mathbf{v}^{0}) \)から、\(\mathbf{h}^{0}\)をサンプルする
  3. \(P(\mathbf{v} \mid \mathbf{h}^{1}) \)から、\(\mathbf{v}^{1}\)をサンプルする
  4. ①〜②を繰り返す

*条件付き独立の性質から各サンプルの要素\(v_{i}^{t}\)を一つ一つサンプルするのではなく、並列的に全ての要素\(\mathbf{v}^{t}\)をサンプルすることができます(\(\mathbf{h}^{t}\)も同様)

イメージ的には、以下のようになります。

 

ギブスサンプリングのイメージギブスサンプリングのイメージ

 

ここで、マルコフ連鎖の性質から、最初の段階は定常分布\(P(\mathbf{v}, \mathbf{h} \mid \mathbf{\theta})\)からサンプルではないので、定常状態になっていない初期のサンプルは標本平均を計算するためのサンプル列の候補から外します。

定常状態に到達するまでのおおよその時間を『バーンイン(burn-in)』と言います。

また、一つのサンプル\(\mathbf{v}^{0}\)をサンプルした後に、\(\mathbf{v}^{1}\)を標本平均を計算するためのサンプル列の候補として選ぶのもNGです。

理由は、直後のサンプル\(\mathbf{v}^{1}\)は\(\mathbf{v}^{0}\)と統計的に相関しているからです。

そのため、サンプルを一つ得たら、十分間隔を開けて新たなサンプルをサンプル列に加える必要があります。

それを避けるための方法として、ランダムな多数の初期値から独立に走らせた複数のマルコフ連鎖から、独立に一つずつサンプルを採用する方法を使用することができます。

しかし、このようなギブスサンプリングを勾配更新毎に実行するのでは計算コストは減ったものの実用的なレベルではありません。

そのためさらに大胆な近似手法であるコンストラスティブ・ダイバージェンス法(Constrastive Divergence ; CD法)を使用します。

 

コンストラスティブ・ダイバージェンス法(Constrastive Divergence ; CD法)

 

コンストラスティブ・ダイバージェンス法(Constrastive Divergence ; CD法)では、ネガティブファイズの評価を以下のK回の連鎖から得られるサンプルの値で近似するという考え方です。

 

CD法の概念図k-CD法のイメージ

 

具体的な勾配更新式は以下のようになります。

k-CD法による勾配更新①

\begin{align}&w_{ij}^{(t+1)} \leftarrow w_{ij}^{(t)} + \eta \big( \mathbb{E}_{\text{data}}[v_{i} h_{j}] – \mathbb{E}_{P_{k}}[v_{i} h_{j}] \big) \\ &b_{i}^{(t+1)} \leftarrow b_{i}^{(t)} + \eta \big(\mathbb{E}_{\text{data}}[v_{i}] – \mathbb{E}_{P_{k}}[v_{i}] \big) \\ &c_{j}^{(t+1)} \leftarrow c_{j}^{(t)} + \eta \big( \mathbb{E}_{\text{data}}[h_{j}] – \mathbb{E}_{P_{k}}[ h_{j}] \big) \end{align}

ここで、\(P_{k}\)は$k$回の連鎖から得られたサンプル列による平均を意味します。

具体的には、以下を意味します。

$$\mathbb{E}_{P_{k}}[\cdots] = \frac{1}{k} \sum_{t=1}^{K} \cdots $$

一般的に、k回の連鎖から得られるサンプルを利用するCD法をk-CD法と言います。

また、CD法では、隠れ変数のサンプル値\(\mathbf{h}^{t}\)を使用するのではなく、条件付き確率\(P(\mathbf{h} \mid \mathbf{v}^{t}, \mathbf{\theta})\)を使用します。

あらに、得られたサンプル列の標本平均ではなく、一つのサンプル\(\mathbf{v}^{k}\)を使用してネガティブフェイズを評価します。

気持ちとしては、勾配上昇法で更新を繰り返すうちに多数のサンプルの効果が間接的に取り込めるという所に起因しています。

ここまでをまとめると、実用的なk-CD法による勾配上昇式は以下のようになります。

k-CD法による勾配更新②

\begin{align}&w_{ij}^{(t+1)} \leftarrow w_{ij}^{(t)} + \eta \big( v_{i}^{\mu} P(h_{j}=1 \mid \mathbf{v}^{\mu}, \mathbf{\theta}) – v_{i}^{k} P(h_{j}=1 \mid \mathbf{v}^{k}, \mathbf{\theta}) \big) \\ &b_{i}^{(t+1)} \leftarrow b_{i}^{(t)} + + \eta \big( v_{i}^{\mu} – v_{i}^{k} \big) \\ &c_{j}^{(t+1)} \leftarrow c_{j}^{(t)} + \eta \big(P(h_{j}=1 \mid \mathbf{v}^{\mu}, \mathbf{\theta}) – P(h_{j}=1 \mid \mathbf{v}^{k}, \mathbf{\theta}) \big) \end{align}

パラメータを更新した後は、訓練サンプルから新たなサンプル\(\mathbf{v}\)を選び更新を繰り返します。

 

ミニバッチ勾配上昇法を利用する場合は、ミニバッチに含まれるデータを初期値とした複数のマルコフ連鎖を構成し、得られたサンプルの標本平均をとり勾配更新を行います。

なんとk-CD法は1-CD法(連鎖を一回)でも十分な機能を持ち実用的なことが経験的に知られています。

このk-CD法は、勾配更新のたびに\(k\)回の連鎖を走らせるだけなので、ギブスサンプリングに比べて計算量を大幅に減少できます。

 

パーシステント・コンストラスティブ・ダイバージェンス法(PCD法)

 

パーシステント・コンストラスティブダイバージェンス法(Persistent Constrastive Divergence)法は、CD法が定常状態からのサンプリングができないという問題を部分的に解決した方法でCD法以上の効率・制度を経験的に確保することができると言われています。

具体的には以下のような更新を行います。

k-PCD法のイメージk-PCD法のイメージ

 

k-PCD法では前のパラメータ更新で使用したサンプルを初期値として利用し、再び勾配更新を行う方法です。

現在、RBMを実装する際はPCD法が使用されることが多いです。

 

制限ボルツマンマシンの誤差関数

 

制限ボルツマンマシンの学習が順調に進んでいるかどうかを確認するためには、対数尤度または、尤度をチェックする必要があります。

しかし、対数尤度・尤度を計算するためには、分配関数を計算する必要があります…

そのため、直接計算することは計算量困難です。

そこで、近似的に対数尤度を評価する方法が開発されています。

擬似対数尤度(Pseudo-likelihood)

 

擬似対数尤度は、全ての可視ユニットが独立であると仮定して計算した対数尤度で、対数尤度近似する手法です。

具体的には、以下のように近似します。

疑似対数尤度

$$\text{PLL}({\bf v}) \equiv \log \prod_{i=1}^{m} p(v_{i} | {\bf v}_{-i}) = \sum_{i=1}^{m} \log p(v_{i} | {\bf v}_{-i})$$

 

ここで、\( {\bf v}_{-i} \)は、\(v_{i}\)以外のユニットの集合です。

この擬似対数尤度は、訓練データ数が大きければ、漸近的に対数尤度に一致することが示されます。

制限ボルツマンマシンの場合、\(p(v_{i}|{\bf v}_{-i})\)は、以下のように表せます。

制限ボルツマンマシンの疑似対数尤度

$$p(v_{i} | {\bf v}_{-i} ) = \sigma \big(F(v_{i}=0, {\bf v}_{-i}) ~ – F(v_{i}=1, {\bf v}_{-i}) \big) $$

 

これでも計算量を大幅に減らすことができますが、入力データの次元が大きい時は依然として計算は大変です。

そのため、擬似対数尤度(PLL)を確率的に近似したものがよく使用されます。

$$\hat{\text{PLL}}({\bf v}) = D \cdot \log \sigma \big(F(\hat{{\bf v}}_{\tau} – F({\bf v})  \big)  $$

  • \( \tau \)は、各データの要素に関する一様分布を持つ確率変数
  • \(\hat{{\bf v}}_{\tau} \)は、\(\tau\)で指定されるユニットを反転したもの
    (e.g. 0 → 1, 1 → 0)

Note : 

$$\mathbb{E}_{\tau} [\hat{\text{PLL}}({\bf v}) ] = \text{PLL}({\bf v})  $$

こうすることで、計算量を大きく削減することができました。

それ以外の評価指標は、

  • Reconstruction error
  • Annealed Importance Sampling
  • Validationデータとtrainingデータの間のaverage自由エネルギー差

 

PytorchによるRBMの実装

 

RBMの学習では行列計算を多く含むためGPUを簡単に利用できるPytorchがおすすめです。

GPUを利用したPytorchによるRBMの実装例を共有していきます。

必要なライブラリをインポート

 

まずは、下記のコードで必要なライブラリをインポートしてください。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import seaborn as sns
# 時間を測る
import time
# 実行時間のプログレスバーを表示
from tqdm.notebook import tqdm
import torch
import torch.nn.functional as F
from torchvision import datasets, transforms
from torch.utils.data import TensorDataset, DataLoader
from torch import nn, optim

 

『time』はエポックあたりの計算時間を測るためにインポートしました。

また、『tqdm』はプログレスバーを表示するためにインポートしました。

 

RBMの実装

 

Pytorchを使用したRBMのコードは具体例は以下のようになります。

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import time
from tqdm.notebook import tqdm
import torch
import torch.nn.functional as F
from torchvision import datasets, transforms
from torch.utils.data import TensorDataset, DataLoader
from torch import nn, optim
import seaborn as sns

class RBM(nn.Module):
    """RBM with PyTorch
    Attributes:
        b : 可視ユニットのバイアス項
        c : 隠れユニットのバイアス項
        w : 重みパラメータ
        n_vis : 可視ユニットの数 
        n_hid : 隠れユニットの数
        epoch : エポック数
        learning_rate : 学習率
        batch_size : ミニバッチサイズ
        initial_std : 初期状態の分散
        seed : 乱数のシード
        device : デバイスを指定
    """
    def __init__(self, n_vis=784, n_hid=128, epoch=10, learning_rate=0.1, 
                 batch_size=100, initial_std=0.01, seed=0, device='cpu'):
        super(RBM, self).__init__()
        self.n_hid = n_hid
        self.n_vis = n_vis
        self.device = device
        self.b = torch.zeros(1, n_vis,  device=device)
        self.c = torch.zeros(1, n_hid,  device=device)
        self.w = torch.empty((n_hid, n_vis), device=device).normal_(mean=0, std=initial_std)
        self.k = k
        self.epoch = epoch
        self.learning_rate = learning_rate
        self.batch_size = batch_size

    def visible_to_hidden(self, v):
        """
        可視ユニットから隠れユニットをサンプル
        """
        p = torch.sigmoid(F.linear(v, self.w, self.c))
        return p.bernoulli()

    def hidden_to_visible(self, h):
        """
        隠れユニットから可視ユニットをサンプル
        """
        p = torch.sigmoid(F.linear(h, self.w.t(), self.b))
        return p.bernoulli()
    
    def com_hiddens(self, v):
        """
        P(h=1|v)を計算
        """
        return torch.sigmoid(F.linear(v, self.w, self.c))

    def sample_v(self, v, gib_num=1):
        """データを生成
        Attributes:
        v : 訓練データ
        gib_num : マルコフ連鎖の更新回数
        """
        v = v.view(-1, self.n_vis)
        v = v.to(self.device)
        h = self.visible_to_hidden(v)
        for _ in range(gib_num):
            v_gibb = self.hidden_to_visible(h)
            h = self.visible_to_hidden(v_gibb)
        return v, v_gibb

    def free_energy(self, v):
        """自由エネルギーを計算
        """
        v_term = torch.matmul(v, self.b.t())
        w_x_h = torch.matmul(v, self.w.t()) + self.c
        h_term = torch.sum(F.softplus(w_x_h), dim=1)
        return -h_term - v_term

    def loss(self, v):
        """疑似対数尤度を計算
        """
        flip = torch.randint(0, v.size()[1], (1,))
        v_fliped = v.clone()
        v_fliped[:, flip] = 1 - v_fliped[:, flip]
        free_energy = self.free_energy(v)
        free_energy_fliped = self.free_energy(v_fliped)
        return  v.size()[1]*F.softplus(free_energy_fliped - free_energy)

    def batch_fit(self, v_pos):
        """ミニバッチあたりの学習更新
        """
        # positive part
        ph_pos = self.com_hiddens(v_pos)
        # negative part
        v_neg = self.hidden_to_visible(self.h_samples) 
        ph_neg = self.com_hiddens(v_neg)
        lr = (self.learning_rate) / v_pos.size()[0]
        # 更新
        update = torch.matmul(ph_pos.t(), v_pos) - torch.matmul(ph_neg.t(), v_neg)
        self.w += lr * update
        self.b += lr * torch.sum(v_pos - v_neg, dim=0)
        self.c += lr * torch.sum(ph_pos - ph_neg, dim=0)
        # PCDのために隠れユニットの値を保持
        self.h_samples = ph_neg.bernoulli()

    def fit(self, train_loader, train_tensor):
        '''学習
        Attributes:
        train_loader : dataloader
        train_tensor : 訓練データ
        '''
        self.losses = []
        self.h_samples = torch.zeros(self.batch_size, self.n_hid, device=device)
        for epoch in tqdm(range(self.epoch)):
            running_loss = 0.0
            start = time.time()
            for id, (data, target) in enumerate(train_loader):
                data = data.view(-1, self.n_vis)
                data = data.to(self.device)
                target = target.to(self.device)
                self.batch_fit(data)
            # 疑似対数尤度を計算
            for _, (data, target) in enumerate(train_loader):
                data = data.view(-1, self.n_vis)
                data = data.to(self.device)
                target = target.to(self.device)
                running_loss += self.loss(data).mean().item()/train_tensor.size()[0]
            self.losses.append(running_loss)
            end = time.time()
            print(f'Epoch {epoch+1} pseudo-likelihood = {running_loss:.4f}, {end-start:.2f}s')

 

コードの概要を以下で簡単に説明します。

  1. 1-PCDを使用した勾配更新
  2. ミニバッチ勾配降下法を使用
  3. 疑似対数尤度を評価基準に採用
  4. 訓練データとDataLoaderを入力として学習を実行

次は、『MNIST』という手書き文字データを使った学習例を紹介します。

 

MNISTを使用した数値実験

 

MNISTとは以下のような画像サンプルです。

MNISTの具体例MNISTデータセット

 

この手書き文字42000枚を訓練データとしてRBMの学習を行います。

まずは、訓練データからDataLoaderを作成してください。

作成後以下のコードを実行することで学習を行うことができます。

# Deviceの指定
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
# RBMインスタンスの作成
model = RBM(initial_std=0.01, learning_rate=0.01, epoch=100, batch_size=100, device=device)
# 学習の実行
losses = model.fit(loader_train, X_train)

 

すると学習が実行され終了すると以下のように表示されます。

学習の出力

 

学習が完了したら下記のコードで訓練データを実際に生成することができます。

images = next(iter(loader_train))[0]
# データを生成
v, v_gibb = model.sample_v(images, gib_num=100)
# 訓練データ
v = v.detach().cpu().numpy()
# 生成されたデータ
v_gibb = v_gibb.detach().cpu().numpy()

 

上記の『v_gibb』にはRBMによって生成された画像が代入されているので、\(28 \times 28\)のサイズに変更し可視化してみると以下のようになります。

RBMに生成されたサンプル

手書き数字に類似した文字が生成されているのがわかります。

また、疑似対数尤度がエポック毎にどのように変化したのかを以下のコードで確認します。

epoch_num = torch.arange(1, model.epoch+1)
losses = torch.Tensor(model.losses)
fig, ax = plt.subplots()
ax.plot(epoch_num, losses, marker='o')
ax.set_xlabel('epoch')
ax.set_ylabel('PLL')
ax.set_xlabel('epoch', fontsize=30)
ax.set_ylabel('PLL', fontsize=30)
plt.tick_params(axis='y', which='major', labelsize=25)
plt.tick_params(axis='x', which='major', labelsize=25)

 

<出力>

疑似対数尤度を確認

まだ、疑似対数尤度は増加途中のように見えますが、着実に増加していることが確認できます。

参考資料

 

参考資料を紹介していきます。

参考文献

 

参考文献を紹介します。

ボルツマンマシン (シリーズ 情報科学における確率モデル 2)

 

制限ボルツマンマシンだけでなく、深層ボルツマンマシンやディープビリーフネット等の応用的なモデルの解説が詳しく載っています。

 

機械学習スタートアップシリーズ これならわかる深層学習入門

 

制限ボルツマンマシン以外にも深層学習に関する基本的な知識を概観できます。

 

深層学習 (機械学習プロフェッショナルシリーズ)

 

後半部分に制限ボルツマンマシンと深層ボルツマンマシンの解説が載っています。

 

参考論文

 

参考にした論文を紹介します。

 

 

まとめ

 

制限ボルツマンマシンの理論とPythonによる実装を説明しました。

筆者が勘違いしている可能性は0ではないので、本文は間違いを含むかもしれません。

もし、間違いや理解の違いを見つけた方は、本記事のコメント欄または、努力のガリレオのTwitterにDMしていただけると助かります。

 

ABOUT ME
努力のガリレオ
【運営者】 : 東大で理論物理を研究中(経歴)東京大学, TOEIC950点, NASA留学, カナダ滞在経験有り, 最優秀塾講師賞, オンライン英会話講師試験合格, ブログと独自コンテンツで収益6桁達成 【編集者】: イングリッシュアドバイザーとして勤務中(経歴)中学校教諭一種免許取得[英語],カナダ留学経験あり, TOEIC650点
【必見】オンライン英会話を利用しないのは時代遅れ?

 

 

英語学習でオンライン英会話を利用しないのは、時代遅れです…

 

『スピーキングができないから、オンライン英会話はまだまだ後回し…』と思っている方は大間違いです。

 

最近のオンライン英会話は、教材が豊富で、もはや参考書・問題集は不要になりつつあります…

また、多くのオンライン英会話が最新メソッドを導入して、スピーキングのみならず、リーディング・リスニング・ライティングの能力も上がる教材を導入しています。

さらに、効率的な記憶ツールや学習カウンセリングのサービスが提供されることもあります!!

 

特に初心者の方は、下記の記事で紹介しているオンライン英会話をチェックしてみてください。

 

【東大生が厳選】初心者にオススメなオンライン英会話10選実は、最適なオンライン英会話はユーザーによって異なります。本記事では、複数のオンライン英会話を経験した私が、他社と比較しつつ特徴を明確にして本当にオススメできるオンライン英会話を10個厳選しました。...

 

中には、2週間の無料体験期間を提供しているオンライン英会話もあります。

英語学習が滞っている方や英語学習を習慣化できていない方は無料体験期間で良いので挑戦してみてください。

何事も行動することが大切です。

 

また、多くのオンライン英会話がTOEICやTOEFL等の外部英語試験の対策レッスンを提供してきています。

中には、追加料金なしで、本番さながらの対策や面接対策を行うことができるオンライン英会話もあります。

 

特にTOEIC対策に強いオンライン英会話は下記を参考にしてください。

 

【徹底比較】TOEIC対策に強いオンライン英会話5選|東大生が厳選!オンライン英会話はTOEICのスコアを効率的に上げる最強のツールです。本記事では、TOEICにオンライン英会話が効果的な理由とTOEIC対策に強いオンライン英会話を厳選しました。詳しい内容は記事を参考にしてください。...

 

TOEFL対策に強いオンライン英会話に関しては下記を参考にしてください。

 

【徹底比較】TOEFL対策に強いオンライン英会話5選|東大生が厳選!TOEFL対策にオンライン英会話は有効です。事実、私もTOEFL対策のためにオンライン英会話を使用してTOEFL ibtの4技能、TOEFL itpのスコアを爆上げすることができました。この記事では、私が思うTOEFL対策に最適なオンライン英会話を紹介しました。...