🥃

あちらのお客様から制限ボルツマンマシンです。

2025/01/24に公開
2

初めに

2024年のノーベル物理学賞は人工知能の分野で受賞されました。「え!?物理学賞なのに人工知能なの!?」と思うかもしれませんが、少なくとも人工知能の根幹であるニューラルネットワークと統計力学は、同様の確率構造で議論され、情報統計力学とういう分野があるくらい、人工知能と物理学は密接につながっています。よって統計力学を勉強してきた筆者からすれば、物理学賞で人工知能が受賞されても不思議ではないと考えています。

今回はこうしたニューラルネットワークにちなんで特に統計力学と接点のあるボルツマンマシンについて述べていこうかなと思います。

ボルツマンマシンとは?

まず何をやりたいか説明します。これは昨今の生成AIのベースの考え方と同じです。とある情報(画像や音声など)がランダムに得られるとき、その情報を生成する元の確率分布q(\mathbf{x})は何かを考えます。具体的には確率分布q(\mathbf{x})から得られたデータ列\mathbf{x}^{(1)},\cdots,\mathbf{x}^{(N)}から元の確率分布q(\mathbf{x})を再現することを考えます。

ボルツマンマシンではこの確率分布を無向グラフをベースに考えていきます。

エポックごとのRMSE

ここで無向グラフ上のエネルギーを次式で定義します。

\Phi(\mathbf{x})=-\sum_{i,j=1}^{n}w_{ij}x_{i}x_{j}-\sum_{i=1}^{n}b_{i}x_{i}

ここでi,jは無向グラフ上のノードの位置を表すもので、図ではi,j\in\{1,2,3,4,5\}です。x_{i}x_{j}も)はベクトル\mathbf{x}\in\mathbb{R}^{n}の要素で各ノードが取る値です[1]。具体的には0または1の2つの値を取ります。nはノードの数をあら合わしており、図に照らし合わせるならばn=5です。b_{i}\in\mathbb{R}はとあるノードが持つエネルギーの大きさを表すパラメータ、w_{ij}\in\mathbb{R}は2つのノード間のエネルギーの大きさを表すパラメータです。ノード間に向きはないのでw_{ij}=w_{ji}です。

ここでとある行列\mathbf{W}\in\mathbb{R}^{n\times n}とベクトル\mathbf{b}\in\mathbb{R}^{n}について、\mathbf{W}=(w_{ij}),\mathbf{b}=(b_{i})と置けば、エネルギーは

\Phi(\mathbf{x})=-\mathbf{x}^{\top}\mathbf{W}\mathbf{x}-\mathbf{b}^{\top}\mathbf{x}

と書くことができます[2]。ここで無向グラフ内のノードが\mathbf{x}という状態を取る確率を

p(\mathbf{x})=\frac{1}{Z}\exp[-\Phi(\mathbf{x})]

と置きます。これはボルツマン分布(またはギブス分布)と呼ばれます。Zは規格化の定数で

Z=\sum_{\mathbf{x}} \exp[-\Phi(\mathbf{x})] =\sum_{x_{1}=0}^{1}\cdots\sum_{x_{n}=0}^{1}\exp[-\Phi(\mathbf{x})]

です。無向グラフにおいてとある状態\mathbf{x}を取る状態が高いというのは、ボルツマン分布の式よりとあるエネルギーを取る値が小さいことが意味します。もし画像・音声といった情報について情報\mathbf{x}周りの情報が良く得られるならば、\mathbf{x}を得る確率は高いはずです。ボルツマンマシンではこの\mathbf{x}がよく得られるようにパラメータであるw_{ij},b_{i}を学習させていくのが目的です。

制限ボルツマンマシンとは?

制限ボルツマンマシンでは無向グラフを可視変数を持つノードと隠れ変数を持つノードに分け、各ノードを可視層と隠れ層にグループ分けします。

エポックごとのRMSE

可視層のノードは実際に外からのデータの受け渡しをこなし、隠れ層はデータの受け渡しはせず、外から見えない層となっております。制限ボルツマンマシンでは可視層内のノードの結合はありません。隠れ層も同様です。結合は可視層と隠れ層間のみ存在します。ここでこの無向グラフの状態\mathbf{x}

\mathbf{x}=\begin{pmatrix} \mathbf{v}\\ \mathbf{h} \end{pmatrix}

可視層の状態(情報)\mathbf{v}\in\mathbb{R}^{n}と隠れ層の状態\mathbf{h}\in\mathbb{R}^{m}に分割することができます。この無向グラフのエネルギーは可視層と隠れ層間の結合がないことを考慮して、

\Phi(\mathbf{v},\mathbf{h};\bm{\theta})=-\mathbf{v}^{\top}\mathbf{W}\mathbf{h}-\mathbf{b}^{\top}\mathbf{v}-\mathbf{c}^{\top}\mathbf{h}

と書けます。ここで\mathbf{W}\in\mathbb{R}^{m\times n},\mathbf{b}\in\mathbb{R}^{n},\mathbf{c}\in\mathbb{R}^{m}というパラメータです。これらのパラメータを\bm{\theta}=(\mathbf{W},\mathbf{b},\mathbf{c})でまとめています。ここで無向グラフが\mathbf{v},\mathbf{h}を取る確率は、条件付確率の形式で、

p(\mathbf{v},\mathbf{h}|\bm{\theta})=\frac{1}{Z(\bm{\theta})}\exp[-\Phi(\mathbf{v},\mathbf{h};\bm{\theta})]

と書けます。ただ実際には我々が得る情報というのは可視層のみの情報\mathbf{v}のみで、隠れ層の事情は知りません。これを確率分布として表現すると隠れ層に関する周辺確率で、

p(\mathbf{v}|\bm{\theta})= \sum_{\mathbf{h}}p(\mathbf{v},\mathbf{h}|\bm{\theta}) =\sum_{\mathbf{h}}\frac{1}{Z(\bm{\theta})}\exp[-\Phi(\mathbf{v},\mathbf{h};\bm{\theta})]

と書けます。制限ボルツマンマシンではとある情報を生み出す真の確率q(\mathbf{v})とこのp(\mathbf{v}|\bm{\theta})がどれだけ近いか調べることになります。

制限ボルツマンマシンの学習

制限ボルツマンマシンで達成すべき事項は次のKLダイバージェンスを最小とするp(\mathbf{v}|\bm{\theta})を得ることです。

D_{KL}(q||p)=\sum_{\mathbf{v}}q(\mathbf{v})\log\frac{q(\mathbf{v})}{p(\mathbf{v}|\bm{\theta})}

これを目的関数として扱い、勾配降下法的にパラメータを更新していきます。具体的には次のようにパラメータを更新していきます。

\bm{\theta}_{t+1}\leftarrow \bm{\theta}_{t}-\eta\frac{\partial D_{KL}(q||p)}{\partial\bm{\theta}}

tは更新回数を表します。\etaは学習率です。ここでKLダイバージェンスは、

D_{KL}(q||p)= \mathbb{E}_{q(\mathbf{v})}\left[\log q(\mathbf{v})\right] -\mathbb{E}_{q(\mathbf{v})}\left[\log p(\mathbf{v}|\bm{\theta})\right]

であり右辺第1項は\bm{\theta}に依存しないので、

\Delta\bm{\theta}=\frac{\partial D_{KL}(q||p)}{\partial\bm{\theta}} =-\frac{\partial \mathbb{E}_{q(\mathbf{v})}\left[\log p(\mathbf{v}|\bm{\theta})\right]}{\partial\bm{\theta}} =-\mathbb{E}_{q(\mathbf{v})}\left[ \frac{\partial\log p(\mathbf{v}|\bm{\theta})}{\partial\bm{\theta}} \right]

となります。この\Delta\bm{\theta}はがんばって計算すると、

\Delta\bm{\theta}= \mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\right] - \mathbb{E}_{q(\mathbf{v})}\left[ \mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]\right]

となります。

詳細な導出

まず\frac{\partial\log p(\mathbf{v}|\bm{\theta})}{\partial\bm{\theta}}について、\frac{\partial p(\mathbf{v},\mathbf{h}|\bm{\theta})}{\partial\bm{\theta}}

\frac{\partial p(\mathbf{v},\mathbf{h}|\bm{\theta})}{\partial\bm{\theta}} =-p(\mathbf{v},\mathbf{h}|\bm{\theta})\frac{1}{Z(\bm{\theta})}\frac{\partial Z(\bm{\theta})}{\partial\bm{\theta}} -p(\mathbf{v},\mathbf{h}|\bm{\theta})\frac{\partial\Phi(\bm{\theta})}{\partial\bm{\theta}}

となるので、

\frac{\partial\log p(\mathbf{v}|\bm{\theta})}{\partial\bm{\theta}} =\frac{1}{p(\mathbf{v}|\bm{\theta})}\frac{\partial p(\mathbf{v}|\bm{\theta})}{\partial\bm{\theta}} =\frac{1}{p(\mathbf{v}|\bm{\theta})}\sum_{\mathbf{h}}\frac{\partial p(\mathbf{v},\mathbf{h}|\bm{\theta})}{\partial\bm{\theta}}

の関係を用いると、

\frac{\partial\log p(\mathbf{v}|\bm{\theta})}{\partial\bm{\theta}} =-\frac{1}{Z(\bm{\theta})}\frac{\partial Z(\bm{\theta})}{\partial\bm{\theta}} -\mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]

となります。式変形の途中でp(\mathbf{h}|\mathbf{v},\bm{\theta})=p(\mathbf{v},\mathbf{h}|\bm{\theta})/p(\mathbf{v}|\bm{\theta})を用いてます。Z(\bm{\theta})を含む項は\mathbf{v}に依存しないので、\Delta\theta

\Delta\bm{\theta}=-\frac{1}{Z(\bm{\theta})}\frac{\partial Z(\bm{\theta})}{\partial\bm{\theta}}- \mathbb{E}_{q(\mathbf{v})}\left[ \mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]\right]

であり、Z(\bm{\theta})を含む項は

-\frac{1}{Z(\bm{\theta})}\frac{\partial Z(\bm{\theta})}{\partial\bm{\theta}} =\sum_{\mathbf{v},\mathbf{h}}p(\mathbf{v},\mathbf{h}|\bm{\theta})\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} =\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\right]

なので、

\Delta\bm{\theta}= \mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\right] - \mathbb{E}_{q(\mathbf{v})}\left[ \mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]\right]

と書けます。

ボルツマンマシンの学習時の問題点

がんばって計算した\Delta\bm{\theta}を考察しましょう。

\Delta\bm{\theta}= \mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\right] - \mathbb{E}_{q(\mathbf{v})}\left[ \mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]\right]

まず第2項目ですが、q(\mathbf{v})を実際の分布ではなく次の経験分布\hat{q}(\mathbf{v})に置き換えることで計算可能です。

\hat{q}(\mathbf{v})=\frac{1}{N}\sum_{\mu=1}^{N}\delta(\mathbf{v} ,\mathbf{v}^{(\mu)})

ここで\delta\mathbf{v}=\mathbf{v}^{(\mu)}1を取りそれ以外は0を取る関数です。\mathbf{v}^{(\mu)}は観測によって得られたデータです。Nは全データ個数です。経験分布は\hat{q}(\mathbf{v})は得られたデータのみを選択する分布なので真の分布q(\mathbf{v})をなぞった確率分布といえます。ならば第2項目は、

\mathbb{E}_{q(\mathbf{v})}\left[ \mathbb{E}_{p(\mathbf{h}|\mathbf{v},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]\right] \approx \frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[ \frac{\partial\Phi(\mathbf{v}^{(\mu)},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \right]

と書くことができます。まだp(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})で期待値を取っている要素が気になりますが、後述しますが、この要素も簡潔に変形できるので、第2項目は容易に計算可能です。

問題は第1項目です。第1項目はp(\mathbf{v},\mathbf{h}|\bm{\theta})で期待値を取っています。実際に期待値を取るという計算を書き下すと、

\begin{align*} \mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\right] &=\sum_{\mathbf{v},\mathbf{h}}p(\mathbf{v},\mathbf{h}|\bm{\theta})\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}}\\ &=\sum_{v_{1}=0}^{1}\cdots\sum_{v_{n}=0}^{1}\sum_{h_{1}=0}^{1}\cdots\sum_{v_{m}=0}^{1}p(\mathbf{v},\mathbf{h}|\bm{\theta})\frac{\partial\Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial\bm{\theta}} \end{align*}

となります。ここでv_{i},h_{i}\mathbf{v},\mathbf{h}の要素です。この式を眺めれば、2^{n+m}回の和の計算をしないといけないことがわかります。つまりノードの数が増えれば増えるほど計算量が指数関数的に増えていくのです。制限ボルツマンマシンの計算にはこのような組合せ爆発の問題があります。

制限ボルツマンマシンの更新方法

前節で置き去りにしたp(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})で期待値をとる箇所を深掘りします。まず\bm{\theta}の要素である\mathbf{W},\mathbf{b},\mathbf{c}について各要素をw_{ij},b_{i},c_{i}と置くと、エネルギーは

\Phi(\mathbf{v},\mathbf{h};\bm{\theta})= -\sum_{i,j}w_{ij}v_{i}h_{j} -\sum_{i}b_{i}v_{i} -\sum_{j}c_{j}h_{j}

と書くことができます。なので

\frac{\partial \Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial w_{ij}}=v_{i}h_{j},\; \frac{\partial \Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial b_{i}}=v_{i} ,\; \frac{\partial \Phi(\mathbf{v},\mathbf{h};\bm{\theta})}{\partial c_{j}}=h_{j}

であり、各要素の更新式は

\begin{align*} \Delta w_{ij}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]\\ \Delta b_{i}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\\ \Delta c_{j}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]\\ \end{align*}

と書けます。

補足(1)

各パラメータの更新式は純粋に代入すれば、

\begin{align*} \Delta w_{ij}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[v_{i}^{(\mu)}h_{j}\right]\\ \Delta b_{i}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[v_{i}^{(\mu)}\right]\\ \Delta c_{j}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[h_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}^{(\mu)}\right]\\ \end{align*}

と書けますが、p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})での期待値を取る個所は\mathbf{h}が対象となるのでv_{i}^{(\mu)}は外に出せます。つまり、

\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[v_{i}^{(\mu)}h_{j}\right]=v_{i}^{(\mu)}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right],\; \mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[v_{i}^{(\mu)}\right]=v_{i}^{(\mu)}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[1\right]=v_{i}^{(\mu)}

であるので、

\begin{align*} \Delta w_{ij}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]\\ \Delta b_{i}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\\ \Delta c_{j}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[h_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]\\ \end{align*}

となるわけです。

\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]の詳細な形を求めましょう。まず確率分布p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})についてですが\mathbf{h}の要素h_{i}について独立です。これは隠れ層同士の結合がないため独立と考えることができます。これは制限ボルツマンマシンの特徴です。よって、

p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta}) =\prod_{i} p(h_{i}|\mathbf{v}^{(\mu)},\bm{\theta})

と書くことができます。したがって、

\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right] =\sum_{\mathbf{h}}p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})h_{j} =\sum_{h_{1}=0}^{1}\cdots\sum_{h_{m}=0}^{1}\prod_{i} p(h_{i}|\mathbf{v}^{(\mu)},\bm{\theta})h_{j}

と書くことができます。h_{j}=0,1p(0|\mathbf{v}^{(\mu)})+p(1|\mathbf{v}^{(\mu)})=1を忘れずに整理すると、

\mathbb{E}_{p(\mathbf{h}|\mathbf{v}^{(\mu)},\bm{\theta})}\left[h_{j}\right]=p(h_{j}=1|\mathbf{v}^{(\mu)},\bm{\theta})

となります。さらに気になるのはp(h_{j}=1|\mathbf{v}^{(\mu)},\bm{\theta})の具体的な形です。結論から言うと、

\sigma(x)=\frac{1}{1+\exp(-x)}

のシグモイド関数を使えば、

p(h_{j}=1|\mathbf{v}^{(\mu)},\bm{\theta})=\sigma(\lambda_{j}^{(\mu)}),\; \left(\lambda_{j}^{(\mu)}=\sum_{i}w_{ij}v_{i}^{(\mu)}+c_{j}\right)

と書くことができます。結果

\begin{align*} \Delta w_{ij}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\sigma(\lambda_{j}^{(\mu)})\\ \Delta b_{i}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[v_{i}\right]-\frac{1}{N}\sum_{\mu=1}^{N}v_{i}^{(\mu)}\\ \Delta c_{j}&=\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[h_{j}\right]-\frac{1}{N}\sum_{\mu=1}^{N}\sigma(\lambda_{j}^{(\mu)})\\ \end{align*}

という更新式になります。こうしてパラメータ更新の第1項目はデータ分だけ総和計算すればよいことがわかりました。

補足(2)

p(h_{i}=1|\mathbf{v}^{(\mu)},\bm{\theta})=\sigma(\lambda_{i})の導出をします。p(\mathbf{h}|\mathbf{v},\bm{\theta})=p(\mathbf{v},\mathbf{h}|\bm{\theta})/p(\mathbf{v}|\bm{\theta})の関係を用います。p(\mathbf{v}|\bm{\theta})について

\begin{align*} p(\mathbf{v}|\bm{\theta})&=\sum_{\mathbf{h}}\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i,j}w_{ij}v_{i}h_{j} +\sum_{i}b_{i}v_{i} +\sum_{j}c_{j}h_{j}\right)\\ &=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i}b_{i}v_{i}\right)\sum_{\mathbf{h}}\exp\left(\sum_{i,j}w_{ij}v_{i}h_{j} +\sum_{j}c_{j}h_{j}\right),\\ &=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i}b_{i}v_{i}\right)\sum_{\mathbf{h}}\exp\left[\sum_{j}\left(\sum_{i}w_{ij}v_{i} +c_{j}\right)h_{j}\right]\\ \end{align*}

であり、\lambda_{j}=\sum_{i}w_{ij}v_{i}+c_{j}と置けば、

p(\mathbf{v}|\bm{\theta})=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i}b_{i}v_{i}\right)\sum_{\mathbf{h}}\exp\left(\sum_{j}\lambda_{j}h_{j}\right)

と整理できます。ここで

\sum_{\mathbf{h}}\exp\left(\sum_{j}\lambda_{j}h_{j}\right) =\prod_{j}(1+\exp\lambda_{j})

となるので、

p(\mathbf{v}|\bm{\theta})=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i}b_{i}v_{i}\right)\prod_{j}(1+\exp\lambda_{j})

と書けます。同様の式変形で、

\begin{align*} p(\mathbf{v},\mathbf{h}|\bm{\theta})&=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i,j}w_{ij}v_{i}h_{j} +\sum_{i}b_{i}v_{i} +\sum_{j}c_{j}h_{j}\right)\\ &=\frac{1}{Z(\bm{\theta})}\exp\left(\sum_{i}b_{i}v_{i}\right)\exp\left(\sum_{j}\lambda_{j}h_{j}\right) \end{align*}

であるので、p(\mathbf{h}|\mathbf{v},\bm{\theta})は、

p(\mathbf{h}|\mathbf{v},\bm{\theta})=\frac{\exp\left(\sum_{j}\lambda_{j}h_{j}\right)}{\prod_{j}(1+\exp\lambda_{j})} =\prod_{j}\frac{\exp\left(\lambda_{j}h_{j}\right)}{1+\exp\lambda_{j}}

よって、p(h_{j}|\mathbf{v},\bm{\theta})は、

p(h_{j}|\mathbf{v},\bm{\theta})=\frac{\exp\left(\lambda_{j}h_{j}\right)}{1+\exp\lambda_{j}}

であり、p(h_{j}=1|\mathbf{v},\bm{\theta})は、

p(h_{j}=1|\mathbf{v},\bm{\theta})=\frac{1}{1+\exp(-\lambda_{j})}=\sigma(\lambda_{j})

となります。

組合せ爆発の対処法

\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[f(\mathbf{v},\mathbf{h})\right]の形の期待値計算について組合せ爆発が起こることは説明しました。こうした問題に対して期待値計算をモンテカルロ積分で近似したいと思います。こんな感じです。

\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[f(\mathbf{v},\mathbf{h})\right]\approx\frac{1}{K}\sum_{\nu=1}^{K}f(\mathbf{v}_{(\nu)},\mathbf{h}_{(\nu)})

ここで\mathbf{v}_{(\nu)},\mathbf{h}_{(\nu)}は何らかの方法でサンプリングした可視層と隠れ層の状態です[3]。ここで問題になるのがこの\mathbf{v}_{(\nu)},\mathbf{h}_{(\nu)}をどうやってサンプリングするかということですが、これはギブスサンプリングの考え方を採用します。

まず制限ボルツマンマシンでは可視層と隠れ層内にグループ分けできることを伝えました。特に制限ボルツマンマシンではp(\mathbf{v}|\mathbf{h},\bm{\theta})p(\mathbf{h}|\mathbf{v},\bm{\theta})の2つの条件付確率を計算でき、双方の条件付確率から\mathbf{v},\mathbf{h}を交互にサンプリングできます。

エポックごとのRMSE

また可視層内のノード同士、または隠れ層内のノード同士に結合はないので、各条件付き確率は、前節の補足(2)で示した計算方法に従えば、

\begin{align*} p(\mathbf{h}|\mathbf{v},\bm{\theta}) &=\prod_{j}p(h_{j}|\mathbf{v},\bm{\theta})= \prod_{j}\frac{\exp\left(\lambda_{j}h_{j}\right)}{1+\exp\lambda_{j}},\;\left(\lambda_{j}=\sum_{i}w_{ij}v_{i}+c_{j}\right)\\ p(\mathbf{v}|\mathbf{h},\bm{\theta}) &=\prod_{i}p(v_{i}|\mathbf{h},\bm{\theta})=\prod_{i}\frac{\exp\left(\lambda_{i}v_{i}\right)}{1+\exp\lambda_{i}},\; \left(\lambda_{i}=\sum_{j}w_{ij}h_{j}+b_{i}\right) \end{align*}

と書くことができます。よって\mathbf{v},\mathbf{h}内の要素v_{i},h_{j}は独立にサンプリングできます。ただこうしたギブスサンプリングであると、計算のたびK個のサンプリングを実施するため計算時間を要する問題があります。そこで提案されているのがコントラスティブ・ダイバージェンス法(CD法)というものです。

CD法ではギブスサンプリングに置ける初期データ\mathbf{v}_{(0)}を実際に観測したデータ\mathbf{v}^{(\mu)}としk回のギブスサンプリングを実施するという方法で、特にkは小さな整数でk=1がよく選ばれます。つまりk\ll Kであるので、複数回のサンプリングの必要性はなくなります。

エポックごとのRMSE

ここで\mu番目の観測したデータについて1回のギブスサンプリングで得た可視層と隠れ層の状態を\mathbf{v}^{(\mu)}_{(1)},\mathbf{h}^{(\mu)}_{(1)}としてモンテカルロ積分を次のように置きます。

\mathbb{E}_{p(\mathbf{v},\mathbf{h}|\bm{\theta})}\left[f(\mathbf{v},\mathbf{h})\right]\approx\frac{1}{N}\sum_{\mu=1}^{N}f(\mathbf{v}^{(\mu)}_{(1)},\mathbf{h}^{(\mu)}_{(1)})

こうすることで効率よく計算することができます。

PyTorchで実装する

実際に実装する前に各パラメータの更新式を明確にします。これまでの記述から次のように書けます。

\begin{align*} \Delta w_{ij}&=\frac{1}{N}\sum_{\mu=1}^{N}\left(v_{i(1)}^{(\mu)}h_{j(1)}^{(\mu)}-v_{i(0)}^{(\mu)}p_{j(0)}^{(\mu)}\right)\\ \Delta b_{i}&=\frac{1}{N}\sum_{\mu=1}^{N}\left(v_{i(1)}^{(\mu)}-v_{i(0)}^{(\mu)}\right)\\ \Delta c_{j}&=\frac{1}{N}\sum_{\mu=1}^{N}\left(h_{j(1)}^{(\mu)}-p_{j(0)}^{(\mu)}\right)\\ \end{align*}

だいぶすっきりした形になりましたね。ここでv_{i(\nu)}^{(\mu)}\mu番目のデータに関する\nu回サンプリングした後のベクトル\mathbf{v}_{(\nu)}^{(\mu)}i番目の要素です。h_{j(\nu)}^{(\mu)}も同様です。またp_{j(\nu)}^{(\mu)}

p_{j(\nu)}^{(\mu)}=\sigma(\lambda_{j(\nu)}^{(\mu)}),\; \lambda_{j(\nu)}^{(\mu)}=\sum_{i}w_{ij}v_{i(\nu)}^{(\mu)}+c_{j}

という確率です。さらにこれらを\Delta\mathbf{W}=(\Delta w_{ij}),\Delta\mathbf{b}=(\Delta b_{i}),\Delta\mathbf{c}=(\Delta c_{j}),\mathbf{v}_{(\nu)}^{(\mu)}=(v_{i(\nu)}^{(\mu)}),\mathbf{h}_{(\nu)}^{(\mu)}=(h_{j(\nu)}^{(\mu)}),\mathbf{p}_{(\nu)}^{(\mu)}=(p_{j(\nu)}^{(\mu)})と置いたならば、

\begin{align*} \Delta \mathbf{W}&=\frac{1}{N}\sum_{\mu=1}^{N}\left[\mathbf{v}_{(1)}^{(\mu)}(\mathbf{h}_{(1)}^{(\mu)})^{\top}-\mathbf{v}_{(0)}^{(\mu)}(\mathbf{p}_{(0)}^{(\mu)})^{\top}\right]\\ \Delta \mathbf{b}&=\frac{1}{N}\sum_{\mu=1}^{N}\left(\mathbf{v}_{(1)}^{(\mu)}-\mathbf{v}_{(0)}^{(\mu)}\right)\\ \Delta \mathbf{c}&=\frac{1}{N}\sum_{\mu=1}^{N}\left(\mathbf{h}_{(1)}^{(\mu)}-\mathbf{p}_{(0)}^{(\mu)}\right)\\ \end{align*}

と書けます。こちらの表現のほうがPyTorchでの実装上イメージしやすいでしょうか?ここでは\mathbf{v}_{(\nu)}^{(\mu)}(\mathbf{p}_{(\nu)}^{(\mu)})^{\top}という表現が行列になることを意識しましょう。

隠れ層の状態のサンプリング方法

とりあえず更新方法は明確になりました。これをコードとして実装すればいいわけですが、その前にどうやってサンプリングすればよいでしょうか?コードとしては以下のようになります。

def sample_hidden_vec(self, v: Tensor) -> tuple[Tensor, Tensor]:
    _lambda = torch.matmul(v, self.W.t()) + self.c
    prob_h = torch.sigmoid(_lambda)
    h_sample = torch.bernoulli(prob_h)
    return prob_h, h_sample

引数のvp(\mathbf{h}|\mathbf{v},\bm{\theta})\mathbf{v}にあたります。CD法を意識するならば\mathbf{v}_{(0)}^{(\mu)}が適切でしょうか。2行目の_lambda\lambda_{j}です。\mathbf{v}_{(0)}^{(\mu)}をもとにサンプリングしたならば、3行目のprob_h\mathbf{p}_{(0)}^{(\mu)}、4行目のh_sample\mathbf{h}_{(0)}^{(\mu)}です。

h_sampleで隠れ層の状態をサンプルする際にはベルヌーイ分布[4]に従ってサンプリングしています。これはtorch.bernoulliで0か1の値に確率的にサンプルしています。

可視層の状態のサンプリング方法

隠れ層のときと同じです。コードを示します。

def sample_visible_vec(self, h: Tensor) -> tuple[Tensor, Tensor]:
    _lambda = torch.matmul(h, self.W) + self.b
    prob_v = torch.sigmoid(_lambda)
    v_sample = torch.bernoulli(prob_v)
    return prob_v, v_sample

引数のhp(\mathbf{v}|\mathbf{h},\bm{\theta})\mathbf{h}にあたります。CD法に則るならば\mathbf{h}_{(0)}^{(\mu)}になります。2行目の_lambda\lambda_{i}です。3行目のprob_vはこれまでの説明では出ませんでしたがp(v_{i}=1|\mathbf{h},\bm{\theta})という確率の要素が格納されているベクトルです。4行目のv_sample\mathbf{v}_{(1)}^{(\mu)}です。

パラメータ更新の実装

早速コードを示します。Nは訓練データのサンプル数です(正確にはバッチ内のサンプル数)。

def cd_method(self, v_0: Tensor, N: int) -> None:
    # 1回サンプリングののCD法
    prob_h_0, h_0 = self.sample_hidden_vec(v_0)
    _, v_1 = self.sample_visible_vec(h_0)
    _, h_1 = self.sample_hidden_vec(v_1)

    # 勾配算出
    self.W.grad = (torch.matmul(h_1.t(), v_1) - torch.matmul(prob_h_0.t(), v_0)) / N
    self.b.grad = torch.sum(v_1 - v_0, dim=0) / N
    self.c.grad = torch.sum(h_1 - prob_h_0, dim=0) / N

1回だけサンプリングするCD法のフローは以下の通りです。

  1. 訓練データv_0から隠れ層の状態h_0とその状態に至った確率prob_h_0をサンプリング。
  2. 隠れ層の状態h_0から新たな可視層の状態v_1をサンプリング。
  3. 最後に可視層の状態v_1から隠れ層の状態h_1をサンプリング。

各状態はv_0, h_0, prob_h_0, v_1, h_1はそれぞれ\mathbf{v}_{(0)}^{(\mu)},\mathbf{h}_{(0)}^{(\mu)},\mathbf{p}_{(0)}^{(\mu)},\mathbf{v}_{(1)}^{(\mu)},\mathbf{h}_{(1)}^{(\mu)}に対応します。これで勾配算出に必要な材料はすべて揃いました。後は提示した勾配算出の式に従ってそのままパラメータを更新するだけです。

全ソースコード

実装したコードは以下の通りです。この制限ボルツマンマシンはオートエンコーダーのように入力した画像を再構成するモデルとして働きます。

https://github.com/torataro-tiger/RBM/blob/main/RBM_run.py

学習結果

学習に用いた訓練データは次の通りです。

  • 28x28 MNIST画像 60000枚(2値化して使用)

学習時のパラメータは次の通りです。

パラメータ
バッチサイズ 32
可視層ノード数 784 (=28x28)
隠れ層ノード数 128
エポック数 10
学習率 0.01

まず各エポックごとのRMSE(二乗平均平方根誤差)を見てみましょう。
エポックごとのRMSE
エポックごとのRMSE(縦軸: RMSE、横軸: エポック)
エポックごとにRMSEは小さくなっているので、うまく学習できるていることが伺えます。テストデータの再構成結果を見てみましょう。
テストデータ再構成結果
テストデータ再構成結果(上段: テストデータ、下段: 再構成結果)
いい感じですね!完全ではないものの再構成できていることがわかります。こんな感じで制限ボルツマンマシンは学習し情報を再構成できることがわかりました!

最後に

こうして理論ベースで話を展開すると、「なぜこんな更新式になるんだろう?」とか「なぜ学習ができるんだろう」といった疑問が解消されるかなと思います。というのも筆者自体が数式で理解しないと完全に理解した気になれないというのがあって、自分が書く記事は濃厚な内容になりがちです。もし何か間違いとかあったら優しく教えてください!では!

参考文献

(1) 一般社団法人日本物理学会「[解説追加しました(10/10)]2024年ノーベル物理学賞は、「人工ニューラルネットワークによる機械学習を可能にした基礎的発見と発明に対する業績」によりJohn J. Hopfield 氏(プリンストン大学、アメリカ)、Geoffrey E. Hinton 氏(トロント大学、カナダ)の2氏が受賞することに決定。」https://www.jps.or.jp/information/2024/10/2024novelprize.php (2025/01/18 閲覧)
(2) 安田宗樹「統計的機械学習理論とボルツマン機械学習」https://aics.riken.jp/labs/cms/workshop/20170322/presentation/yasuda.pdf (2025/01/18 閲覧)
(3) 麻生英樹, 安田宗樹, 前田新一, 岡野原大輔, 岡谷貴之, 久保陽太郎, ボレガラ・ダヌシカ(著),神嶌敏弘(編), 人工知能学会(監)「深層学習 ― Deep Learning」近代科学社 (2017)
(4) 田中和之(編)「確率的情報処理と統計力学 様々なアプローチとそのチュートリアル」サイエンス社 (2006)

脚注
  1. ベクトルといえば列ベクトルを指していると思ってください。 ↩︎

  2. \topは転置で列ベクトルから行ベクトルに変換しています。また\mathbf{W}\in\mathbb{R}^{m\times n}と書いたならばmn列の行列を指していると思ってください。 ↩︎

  3. サンプリングしたデータは\mathbf{v}_{(i)}のように下添え字で、観測したi番目のデータまたは訓練データは\mathbf{v}^{(i)}のように上添え字で表すことにします。 ↩︎

  4. ここの主張は参考文献(3)の3.4.3に記載してます。 ↩︎

Discussion

shaobashaoba

普段研究でAIを使っているのにも関わらず、ボルツマンマシンについては一切知らなかったため非常に興味深かったです。
1点質問があるのですが、通常パラメータ更新でよく用いられる誤差逆伝播法と今回紹介されている制限ボルツマンマシンの更新方法は本質的に全く異なるのでしょうか?
自分は研究でよく学習過程にサンプリングが絡む機械学習モデルを扱っており、誤差逆伝播ができない構造をどう回避するかを日頃から考えていたので気になりました。

寅太郎寅太郎

ご質問ありがとうございます!
まず筆者は面白いから勉強しただけで、専門ではないことにご留意いただければと思います。

まず誤差逆伝搬法が必要なDNNとこの記事で紹介したボルツマンマシンのパラメータ更新は両方とも\bm{\theta}_{t+1}\leftarrow\bm{\theta}_{t}-\eta\nabla Lで更新していくので本質の部分は同じととらえて問題ないと思います。

ただDNNに関してはシグモイド関数のような非線形関数の合成関数で多層化されているため、出力層から入力層に勾配を伝搬するには誤差逆伝搬法の通り、微分の連鎖律の原理を使って計算する必要があります。対してこの記事で紹介したボルツマンマシンはそのような構造になっておらず、むしろエネルギーが小さくなる方向にパラメータをいっぺんに更新していきます。

DNNと似たボルツマンマシンとしてディープボルツマンマシン(DBM)というのがあります。こちらは隠れ層が多層化されている構造ですが、この複雑な構造においてどのようにサンプリングしパラメータ更新を行っているか理解するとshaoba様のご研究に役立てるかもしれません。