この記事は 
Sansan Advent Calendar 2021の25日目 
の記事です。
クリスマスなので、例によって物思いにふける訳ですが、今回は学生時代のある思い出を回想してみようと思います。
この記事で解説する情報統計力学は、2024年ノーベル物理学賞と関連しています 
https://www.nobelprize.org/prizes/physics/2024/popular-information/
情報統計力学との出会い
物理学を専攻していた学生時代、大学の図書館が好きでいろんな棚を物色していた時、西森秀稔氏のスピングラス理論と情報統計力学という本に出会いました。
情報統計力学って何？って感じで、情報理論と統計力学（ミクロな粒子の熱力学、確率論的な扱いによって現象を説明する物理理論）がどう混ざるのか気になったのが始まりです。
その本では、統計力学と量子力学を使ってデジタル画像を復元するという話が展開されていました。
もうかなり衝撃的です。
画像って自然のものではなくて、人工物じゃないですか？
それを自然を記述する物理理論で取り扱うって意味がわからなかったです。
読み進めていくうちに、なるほど。
実際には画像を物理量として直接扱うのではなく、量子力学で頻出するスピンという電子などの粒子が持つ性質に対して行われるモデル化の方法・考え方を応用して、理論展開するというものでした。
とはいえ、物理で学んだことが情報分野で応用されるというものは、自分が報われたというか、何だか嬉しい気持ちになりました。
もう10年以上前の話ですが、今でもこのときの衝撃が鮮明に記憶に残ってます。
物理の話で恐縮ですが、ランダウ本を初めて読んで最小作用の原理を知ったときに匹敵する衝撃です。
そこで、今回はその物理理論が画像修復にどう応用されているのか、ブラウザアプリを通してお伝えできればと思います。
テーマ
オリジナルの画像（元画像）が、何らかの原因によって劣化してしまった状態を考えます。
例えば、以下のように。
| 元画像 | 劣化画像 |
|---|---|
 |
 |
そして、私達はその劣化した画像だけしかない状態（オリジナル画像を知らない状態）でどれだけオリジナル画像に近い画像を再現できるか？
というテーマになります。
オリジナル画像を知らない状態というのは実際によくあることで、誰かがアップロードした後の画像とか、USBメモリに保存された遠い過去のデジカメ写真とかFaxで受信した手書きの絵とかです。
通信ネットワーク環境の不備やハードウェア的な欠陥、熱や放射線による影響でデジタルの0 or 1の値が書き換わってしまう、みたいな状況です。
実際はエラー訂正技術によって回復されるので、今回のテーマは身の回りで実際に応用されている技術ではないです。
画像修復の仕組み
以下の画像を見た時、人間にはきっと「2つの三角形が合体したような、砂時計が傾いたような画像」に見えます。
なぜでしょうか？
考え方は非常にシンプルです。
白が多めのエリアにある規則性のない黒いドットがノイズで、反対に黒い中にポツンとある白もノイズ、であると認識しているからです。
これはつまり、あるドットだけ見た時、周りと同じ色じゃなければ、そのドットは間違っていると考えられます。
そして、すべてのドットに対して、そのような判断を実行して、間違ってそうなら塗り替える、という作業を繰り返す、という工程を複数回反復して行うと、ある時点でそれ以上は結果が変わらないだろうと考えられます。
その時点の画像を修復結果として採用する、という考え方で修復してみる、というのが本テーマのソリューションになります。
実は、あるドットの周りを見て判断する、というのが物理学でも実際に使う計算手法（平均場近似）になっています。
物理学では、そのドットを0 or 1のような離散的な2つの値を持つ、自然界に存在する量（スピン）に対して行います。
そのスピンに対する計算を画像修復で応用するのです。
シミュレーター
詳細に入る前に、今回の理論で実際に修復された画像を確認できるアプリを用意しました。
オリジナル画像はこちらで用意したテンプレしかないです。すみません。
demo
Sandbox
プレビューを開始できたら、画面右上のOpen in a new tabをクリックすると別タブで操作できます 
See the Pen Image Restoration by hiroshi ueda (@hermiltonian) on CodePen.
遊び方
簡単2ステップ
- 劣化させるボタンクリック
- 修復するボタンクリック
また、それぞれのパラメータを変化させながら、修復画像がどう変わるかお楽しみいただければと思います。
| パラメータ | 意味 |
|---|---|
| 劣化確率 | 各画素が劣化する確率です |
| 西森パラメータ | 劣化確率$p$を用いて、$\frac{1}{2},\log \frac{1-p}{p}$と表されるものです。物理で言う逆温度（$1/(kT)$）に対応する量です。 |
| 修復逆温度 | 画像修復過程で使う(逆)温度の初期値です。大きい値ほど隣り合う画素同士が同じ状態(色)である確率が高くなります。 |
| 交換定数 | 交換相互作用の文脈で登場する定数です。修復逆温度同様、大きければ大きいほど隣の画素と同じ状態になりやすくなります。 |
| 誤差 | 数値計算上の許容誤差です |
| 反復回数 | アニーリングによる反復計算の回数です |
理論のコンセプト
理論の詳細を語るには、背景となる物理知識や数式が頻出するので、詳細は割愛して以下ではコンセプトだけまとめておきたいと思います。
詳細が気になる方は、以下のリファレンスから当たるといいと思います。
- 書籍：スピングラス理論と情報統計力学、西森秀稔
- 有限温度の画像修復(多体問題としての情報処理-統計力学と情報科学の接点-,研究会報告)、西森秀稔
- 平均場近似の数理 -確率伝搬法という名の物理と情報の交差点-、 田中和之
物理学の予備知識
磁気モーメントというものがあります。
これは磁石などの磁力を表す量で、その強さと方向を併せ持つベクトルになります。
下記のリンクの図が分かりやすいのですが、棒磁石をそれ以上分割できなくなるまで分割していくと、原子磁石とも言われる磁力を持つ原子に至ります。
その原子磁石が持つ磁力の源の一つとして、電子の自転（スピン）による影響があります。
詳しくはこちらをどうぞ。
磁気モーメント自体は、コンパスと同じ様に磁場（磁界）の影響を受け、磁場と同じ向きに並ぶ性質があります。
これは、磁場に抗うことによる無駄なエネルギー消費を抑えたい、という自然が要求する法則の結果です。
基本的に物理学では、全体のエネルギーが極力小さくなる状態に自然と移行するようにできています。
そこで、状態を記述するのにエネルギーをよく用います。
$n$という状態のエネルギーを$E_{n}$と表した時、原子のような微視的な世界でエネルギー$E_{n}$が実現する確率は以下のカノニカル分布によって表現されます。
\Pr [n] = \frac{e^{-\beta E_{n}}}{\sum\limits_{n} e^{-\beta E_{n}}}
ここで、$\beta$はボルツマン定数$k$と温度$T$を使って、$\beta = 1/(kT)$と表される量（逆温度）です。
電子のスピンは上向き（↑）と下向き（↓）と呼ばれる2つの状態があります。
これが以下で導入する白黒画像の白と黒にちょうど対応してきます。
2つの水素原子を考えた時、それぞれの原子核の周りを飛び回っていた電子は、もう一方の電子と入れ替わるような動きをすることがあります。
これを交換相互作用といいます。この相互作用の強さを表す量を$J$としたとき、2つのスピン$s_{1}$, $s_{2}$の間のエネルギー（交換エネルギー）は
E = -2 J \mathbf{s}_{1} \cdot \mathbf{s}_{2}
となります。上記はベクトルの内積で表現していますが、スピンのようにz軸方向に上向き、下向きの2つの状態しかとれないと仮定すると、単に以下のように表すことができます。
E_{1,2} = -2 J s_{1} s_{2}
このような単純なモデルとイジングモデルといいます。
原画像と劣化画像
状態が2つしかない情報を伝達したとき、ある確率で状態がもう一方に遷移するような通信路を二元対称通信路といい、最も簡単な誤り通信路のモデルとなります。
白黒画像を二元対称通信路を経由して送信したとき、各々の画素はある確率でその白黒を反転させ、受信者のもとへ送られます。
元々の画像内の画素$i$の白黒を$\xi = \pm 1$で表現することにします。
この表現はイジングスピンのものと同じです。
そして、すべての画素つまり画像を
\mathbf{\xi} = (\xi_{1}, \xi_{2},\dots, \xi_{N})
で表すことにします。ここで、$N$は画素数を表しています。
また、画像$\mathbf{\xi}$を受信した結果を
\mathbf{\tau} = (\tau_{1}, \tau_{2},\dots, \tau_{N})
とします。
各$\tau_{i}$は$\pm \xi_{i}$の値をとり、確率$p$で$\tau_{i}=-\xi_{i}$、すなわち反転しているとします。
したがって、$\tau_{i}=\xi_{i}$となる確率は$1-p$です。
仮定から、原画像$\mathbf{\xi}$の約$pN$個の画素が反転して画像が乱れることから、$\mathbf{\tau}$を劣化画像と呼びます。
今回の目的は、劣化画像$\mathbf{\tau}$から原画像$\mathbf{\xi}$を推定することです。
画像におけるカノニカル分布
原画像と劣化画像の、同じ一つの画素に注目します。
原画素$\xi_{i}$が反転しなければ、その画素は送受信で安定的であったとみなせ、一方反転するとその画素、すなわちイジングスピンは交換エネルギーだけ高い状態へと遷移したとみなすことができます。
電子のスピンは上向きが$1/2$、下向きが$-1/2$という値を取るので、画素とは$s=\xi /2$という関係があります。
よって、交換エネルギーは$-J\tau_{i}\xi_{i}/2$と書けるので、このスピンを使ったカノニカル分布は
\Pr[\tau_{i} | \xi_{i}] = \frac{e^{\beta J\tau_{i}\xi_{i}/2}}{\sum\limits_{\tau_{i}\xi_{i} = \pm 1} e^{\beta J\tau_{i}\xi_{i}/2}}
= \frac{e^{\beta J\tau_{i}\xi_{i}/2}}{2 \cosh{(\beta J/2)}}
となります。
ここで、$\tau_{i}\xi_{i} = \pm 1$の状態に関して和をとったのは、起こりうる状態が反転状態$\tau_{i}\xi_{i} = -1$かそうでない状態$\tau_{i}\xi_{i} = 1$の2つしかないからです。
上の式で表された確率は、$\beta=1/(kT)$という形で温度$T$や交換積分$J$などといった物理量が含まれ、それらの画像に対する定義がない今の段階ではあまり意味のない表式となってしまいます。
そこで、反転確率$p$を含んだ新たなパラメータを導入します。
先程の式内の$\beta J/2$を無次元のパラメータ$\beta_{p}$に置き換え、
\beta_{p} = \frac{1}{2}\,\log{\frac{1-p}{p}}
と定義します。
これを上式に代入すると、以下のようになります。
\Pr[\tau_{i} | \xi_{i}]
= \frac{e^{\beta_{p} \tau_{i}\xi_{i}}}{2 \cosh{\beta_{p}}}
このカノニカル分布を使って、反転する確率$\Pr [\tau_{i}=-1 | \xi_{i} = +1]$を求めてみます。
e^{2\beta_{p}} = \frac{1-p}{p},\qquad e^{-2\beta_{p}} = \frac{p}{1-p}
であることに注意すると、
\Pr[\tau_{i}=-1 | \xi_{i} = +1]
= \frac{e^{-\beta_{p}}}{2 \cosh{\beta_{p}}}
= \frac{e^{-\beta_{p}}}{e^{\beta_{p}} + e^{-\beta_{p}}}
= \frac{1}{e^{2\beta_{p}} + 1}
= \frac{1}{\frac{1-p}{p} + 1}
= p
となります。一方、
\Pr[\tau_{i}=+1 | \xi_{i} = +1]
= \frac{1}{1 + e^{-2\beta_{p}}}
= \frac{1}{1 + \frac{p}{1-p}}
= 1-p
であり、パラメータ$\beta_{p}$を用いたことでカノニカル分布から定義のない物理量を消すことができ、さらに反転確率を用いて式を書き換えることができました。
式$\beta_{p} = \frac{1}{2} \log{\frac{1-p}{p}}$で定義されたパラメータ$\beta_{p}$は西森（Nishimori）パラメータと呼ばれます。
以上の一つの画素についての考察から、画像に対しての知見を得ることができます。
原画像$\mathbf{\xi}$を送信して劣化画像$\mathbf{\tau}$を得る条件付きカノニカル分布$\Pr[\mathbf{\tau} | \mathbf{\xi}]$は、
同一画像内で画素同士の反転が独立であるならば、
\Pr[\mathbf{\tau} | \mathbf{\xi}]
= \prod_{i} \Pr[\tau_{i} | \xi_{i}]
= \frac{e^{\beta_{p} \sum\limits_{i} \tau_{i}\xi_{i}}}{(2 \cosh{\beta_{p}})^{N}}
となります。
ベイズの定理
事象$A$, $B$が同時に起こる確率$P(A, B)$は、$B$が起こったときの$A$が起こる条件付き確率$P(A|B)$を用いて
P(A, B) = P(A|B) P(B)
と書けます。
この表式は、$A$, $B$を入れ替えても同じですから、
P(A|B) P(B) = P(B|A) P(A)
\iff P(A|B) = \frac{P(B|A) P(A)}{P(B)}
が成り立ちます。
ここで、$P(B)$は
P(B) = \sum_{A} P(B|A) P(A)
と表せることに注意すると、
P(A|B) = \frac{P(B|A) P(A)}{\sum\limits_{A} P(B|A) P(A)}
となります。この式をベイズ（Bayes）の公式と呼びます。
この式は、事象$A$の下で$B$が起こる事象を観測後、その確率を以って$A$であった確率を求める式です。
すなわち、事前分布$P(A)$ののち、$B$が起こることで$P(A)$は観測者にとって$P(A|B)$という事後分布へと更新されます。
この操作をベイズ改訂といいます。
これを利用すると、事前分布$P(A)$を知らずとも初期値$P(A)$を適当に定めることで$P(B|A)$を求め、これを$P(A)$の推定に用いることができます。
この方法をベイズ推定といいます。
画像修復に以上のベイズの定理を応用して、
\Pr[\mathbf{\xi} | \mathbf{\tau}]
= \frac{\Pr[\mathbf{\tau} | \mathbf{\xi}] \Pr[\mathbf{\xi}]}{\sum\limits_{\mathbf{\xi}} \Pr[\mathbf{\tau} | \mathbf{\xi}] \Pr[\mathbf{\xi}]}
となりますが、事前分布$\Pr[\mathbf{\xi}]$は受信者には知り得ないため、その推定値$\Pr[\mathbf{\sigma}]$で置き換えざるを得ません。
このようにして導入された、画像$\mathbf{\xi}$の推定値$\mathbf{\sigma}$を修復画像といいます。
\Pr[\mathbf{\sigma} | \mathbf{\tau}]
= \frac{\Pr[\mathbf{\tau} | \mathbf{\sigma}]\Pr[\mathbf{\sigma}]}{\sum\limits_{\mathbf{\sigma}} \Pr[\mathbf{\tau} | \mathbf{\sigma}] \Pr[\mathbf{\sigma}]}
ここで、$\Pr[\mathbf{\tau} | \mathbf{\sigma}]$は前節最終式で与えられる$\Pr[\mathbf{\tau} | \mathbf{\xi}]$の推定値ですが、この式をそのまま使うことはできません。
というのも、前節最終式中の西森パラメータ$\beta_{p}$は反転確率$p$を含んでいますが、これは受信者にはわからない情報です。
そこで、$q$を$p$の推定値として、
\beta_{q} = \frac{1}{2}\,\log \frac{1-q}{q}
とし、これを西森パラメータと置き換えます：
\Pr[\mathbf{\tau} | \mathbf{\sigma}]
= \frac{e^{\beta_{q} \sum\limits_{i} \tau_{i}\sigma_{i}}}{(2 \cosh{\beta_{q}})^{N}}
事前分布の推定値$\Pr[\mathbf{\sigma}]$は依然不明なままであるので、経験則的な分布を導入することがよくなされます。
事前分布にしばしば用いられるモデルは、局所的な一様性を要求する方法です。
これは、画像の黒い領域の中に一つだけ白い画素はないであろうという、冒頭に触れた経験に基づくモデルです。
式内には$\beta J = J/(kT) = K$なるパラメータが現れますが、この値は交換相互作用の強さを制御するパラメータとみなせ、これを大きな値にすることで経験則の影響を大きくすることができます。
以上から、事前分布を
\Pr[\mathbf{\sigma}] = \frac{e^{K \hspace{-.5em} \sum\limits_{<i,j>} \hspace{-.5em} \sigma_{i} \sigma_{j}}}{\sum\limits_{\mathbf{\sigma}} e^{K \hspace{-.5em} \sum\limits_{<i,j>} \hspace{-.5em} \sigma_{i} \sigma_{j}}}
と仮定することができます。
これで事後分布をイジングスピン変数を用いて表すことことができます。
\Pr[\mathbf{\sigma} | \mathbf{\tau}]
= \frac{e^{\beta_{q} \sum\limits_{i} \tau_{i}\sigma_{i} + K \hspace{-.5em} \sum\limits_{<i,j>} \hspace{-.5em} \sigma_{i} \sigma_{j}}}{\sum\limits_{\mathbf{\sigma}} e^{\beta_{q} \sum\limits_{i} \tau_{i}\sigma_{i} + K \hspace{-.5em} \sum\limits_{<i,j>} \hspace{-.5em} \sigma_{i} \sigma_{j}}}
この式が、局所一様的事前分布モデルでの事後分布である修復画像のベイズ推定値となります。
具体的な解法
この式を具体的に解く方法がさまざま研究され、今回は平均場アニーリング反復法という手法で修復画像を推定することにしました。
アニーリングとは、温度の高い状態から徐々に下げていき、エネルギーがより小さい状態に落ち着くのを期待する手法です。
話が飛びますが、世界初の商用量子コンピューターのD-Waveは、量子アニーリングという方式の量子コンピュータです。
量子アニーリングの理論を構築したのがまさに西森氏です。
シミュレーターの実装
今回、シミュレーターをSvelteで実装しました。
自分はコンポーネント指向のフロントエンドフレームワークを実務で使ったことがなかったのですが、それでも非常に簡単に実装することができました。
以下にソースコードを共有しておきます。
※後に、SvelteKitによる実装に移行しました
実装で参考にしたリンク
さいごに
今日は2021/12/25 クリスマス ですね。
少しでもその気分を味わおうと、シミュレーターに小細工を仕込みました、というほど小細工ではなくプルプル主張してますが、お楽しみいただければ嬉しいです 