A Blog Entry on Bayesian Computation by an Applied Mathematician
$$
$$
1 導入
1.1 はじめに
MP は情報理論では低密度パリティ検査符号の文脈で (Gallager, 1962) が，因果推論の文脈で (Pearl, 1982) が独立に開発している．
信念伝搬法は変分手法とは違い，ある種のモデル（木やランダムグラフなど）では正確な推論が可能になる上に，一般に速い．1 当然 Monte Carlo 法よりも速い．
そのため，変分推論の改良の源泉が，メッセージ伝搬を通じて探求されている (Winn and Bishop, 2005)．
1.2 ランダムグラフ上のスピン系の分配関数
頂点数 \(N\)，辺数 \(E\) の Erdös-Renyi ランダムグラフ \(\mathcal{G}(N,E)\) 上のスピン \((\sigma_i)_{i\in[N]}\) と相互作用 \((\psi_{ij})_{(i,j)\in E}\) を考える： \[ \psi_{ij}(\sigma_i,\sigma_j)=\exp\biggr(-\beta J_{ij}\sigma_i\sigma_j\biggl). \]
ランダムグラフでは平均的なループの長さは \(\log N\) であり，局所的に木の構造を持つため，このことを利用した近似を使って信念伝搬を行う．
\(Z_{i\to j}(\sigma_i)\) を，\(i\) を根とする木から \(j\) を根とする木を除去し，\(i\) におけるスピンの値が \(\sigma_i\) であると条件付けた場合の部分木上の分配関数，\(Z_i(\sigma_i)\) を \(i\) におけるスピンの値が \(\sigma_i\) であると条件付けた場合の木全体の分配関数とすると，次の関係がある： \[ Z_i(\sigma_i)=\prod_{j\in\partial i}\left(\sum_{\sigma_j}Z_{j\to i}(\sigma_j)\psi_{ij}(\sigma_i,\sigma_j)\right). \]
するとこのとき， \[\begin{align*} \eta_i(\sigma_i):=\frac{Z_i(\sigma_i)}{\sum_{\sigma'}Z_i(\sigma')}&=\frac{1}{\sum_{\sigma'}Z_i(\sigma')}\prod_{j\in\partial i}\left(\sum_{\sigma_j}Z_{j\to i}(\sigma_j)\psi_{ij}(\sigma_i,\sigma_j)\right)\\ &=:\frac{1}{z_i}\prod_{j\in\partial i}\left(\sum_{\sigma_j}\eta_{j\to i}(\sigma_j)\psi_{ij}(\sigma_i,\sigma_j)\right). \end{align*}\] は Boltzmann-Gibbs 分布の周辺分布になっている： \[ m_i=\langle\sigma_i\rangle=\sum_{\sigma}\eta_i(\sigma)\sigma. \]
このときの正規化定数は \[\begin{align*} z_i&=\sum_{\sigma_i}\prod_{j\in\partial i}\left(\sum_{\sigma_j}\eta_{j\to i}(\sigma_j)\psi_{ij}(\sigma_i,\sigma_j)\right) =\sum_{\sigma_i}\prod_{j\in\partial i}\left(\sum_{\sigma_j}\frac{Z_{j\to i}(\sigma_j)}{\sum_{\sigma'}Z_{j\to i}(\sigma')}\psi_{ij}(\sigma_i,\sigma_j)\right)\\ &=\sum_{\sigma_i}\frac{Z_i(\sigma_i)}{\prod_{j\in\partial i}\sum_{\sigma_j}Z_{j\to i}(\sigma_j)}=\frac{Z}{\prod_{j\in\partial i}\sum_{\sigma_j}Z_{j\to i}(\sigma_j)},\\ z_{j\to i}&=\frac{\sum_{\sigma_j}Z_{j\to i}(\sigma_j)}{\prod_{k\in\partial j\setminus i}\sum_{\sigma_k}Z_{k\to j}(\sigma_k)}. \end{align*}\]
と表せる．
\(z_i\) の分母には \(Z\) が現れていることを用いると \(Z\) は，任意の位置 \(i\in[N]\) を起点として，次のような展開表示ができる： \[\begin{align*} Z&=\sum_{\sigma_i}Z_i(\sigma_i)=z_i\prod_{j\in\partial i}\left(\sum_{\sigma_j}Z_{j\to i}(\sigma_j)\right)=z_i\prod_{j\in\partial i}\left(z_{j\to i}\prod_{k\in\partial j\setminus i}\sum_{\sigma_k}Z_{k\to j}(\sigma_k)\right). \end{align*}\] ここで，最右辺の \(\sum_{\sigma_k}Z_{k\to j}(\sigma_k)\) は \(z_{j\to i}\) を計算するのと同じ要領で計算されることになり，最終的に木の葉まで再帰的に計算することで，公式 \[ Z=z_i\prod_{j\in\partial i}\left(z_{j\to i}\prod_{k\in\partial j\setminus i}z_{k\to j}\cdots\right)=z_i\prod_{j\in\partial i}\left(\frac{z_j}{z_{ij}}\prod_{k\in\partial j\setminus i}\frac{z_k}{z_{jk}}\cdots\right)=\frac{\prod_{i\in[N]}z_i}{\prod_{(i,j)\in E}z_{ij}}. \]
よって，この再帰的計算により，自由エネルギー \[ F=-T\log Z=\sum_{i\in[N]}\biggr(-T\log z_i\biggl)-\sum_{(i,j)\in E}\biggr(-T\log z_{ij}\biggl). \] も得られることになる．
1.3 無限レンジ極限での高温解の安定性
(Thouless et al., 1977) は TAP 方程式により，(Sherrington and Kirkpatrick, 1975) 模型を解いた．
その際の結果が，同じく無限サイズのグラフである Bethe 格子 上でならば，信念伝搬法により議論でき，有効 cavity 場が次のように与えられるという： \[ \beta h_0=\sum_{i=1}^k\operatorname{atanh}\biggr(\tanh(\beta J_{ij})\tanh(\beta h_i)\biggl)+H. \]
このモデルでもスピングラス相が出現することが，信念伝搬法が失敗する（局所不安定になる）こととして現れる．これはグラフの木構造を破壊するような長距離の相関が出現することによる．スピングラス相と常磁性相の境界は Almeida-Thouless 線 (Thouless, 1986) という．
2 共同体抽出
2.1 はじめに
グラフの頂点をクラスタリングした際のクラスターに当たる概念を 共同体 (community) という．
最初におもつくような方法は，横断する辺の数が最小になるようなクラスター境界を決定する方法であるが，これは NP 困難である．
これを少し修正して，各辺に対して edge centrality を計算し，これを順に脱落させて edge centrality を計算しなおすというような反復を行う Girvan-Newman アルゴリズム (Girvan and Newman, 2002) が最も有名なものである（引用数２万！）．
2.2 modularity 最大化による方法
Girvan-Newman アルゴリズムの中では，分割の「良さ」の指標として modularity \[ Q=\frac{1}{2m}\sum_{ij}\left(A_{ij}-\frac{d_id_j}{2N}\right)\delta_{s_i}\delta_{s_j} \] を定義し，アルゴリズムの停止の指標としていた．
現在では，逆にこの modularity を最適化することでクラスタリングを達成する手法が主流の１つである．
辺を全て消去した状況から１つずつ追加していき，この modularity の値を最大化する階層的クラスタリングアルゴリズムとして (Newman, 2004), (Clauset et al., 2004) が提案された．
最適化問題として定式化された以上，模擬アニーリングが用いられることもある (Guimerà et al., 2004)．
2.3 スピングラス対応
じきに (Reichardt and Bornholdt, 2006) によって，共同体抽出の問題は Potts 模型の基底状態探索と等価であることが示された．
実際，各頂点に割り当てられた未知の Potts スピン \(\sigma_i\) に対して，同じ状態を持つ頂点同士は繋がろうとし，違う頂点を持つ場合は結合は疎になる傾向があるという状況は \[ H(\sigma)=-\sum_{i<j}J_{ij}\delta(\sigma_i,\sigma_j), \] \[ J_{ij}:=J\biggr(A_{ij}-\gamma p_{ij}\biggl) \] というハミルトニアンで表現できる．\(A_{ij}\) は隣接行列 \(A=(A_{ij})\in M_n(2)\) の成分，\(p_{ij}\) は辺の数の期待値，\(\gamma\) は推定されるクラスタ数を決定するハイパーパラメータ (resolution parameter とも表現するらしい) である．
こう捉えると，たしかに模擬アニーリングは１つの選択肢だ．
2.4 信念伝搬による復号
Potts モデルではないが，(Hastings, 2006) は同様にスピングラスモデルと同一視をし，こちらでは信念伝搬により基底状態探索を行った．
一方で (Newman and Leicht, 2007) は混合モデルとみなし，EM アルゴリズムによる方法を見出している．ここまで来ると確率的ブロックモデル (Snijders and Nowicki, 1997) による方法に非常に似通っており，(Decelle et al., 2011) はまさにこの見方をしている．
そもそも，EM アルゴリズムと西森ラインは深い関係がある．2 実際，EM アルゴリズムも，確率的ブロックモデルにおけるクラスの割り当ても，\(k\)-平均法のクラスタ数も，相転移を起こす．背後に何か物理過程と対応がつくものが存在するのかもしれない．
ブロックモデル (blockmodel / image graph) とはグラフの分割に関するモデルで，クラスタ数 \(q\)，\([q]\) 上の確率分布（頂点数の分布）\(\{n_\alpha\}_{\alpha\in[q]}\)，グループ間に辺が存在する確率を表す行列 \((p_{ab})_{a,b\in[q]}\) のパラメータからなる．3
2.5 スペクトルを通じた方法
モジュラリティ最適化と対照的な手法として，スペクトルを用いた方法がある．有名なサーベイには (Luxburg, 2007) がある．
最も直接的には，Laplacian 行列 の固有空間分解を通じて頂点集合を別の潜在空間に埋め込み，そこで \(\mathbb{R}^N\)-クラスタリングを行うという (Donath and Hoffman, 1973) 以来の方法である．この文献はグラフの分割を問題にしているが，(Donetti and Muñoz, 2004) はコミュニティ抽出に集中している．
こちらの方が数学的にグラフの構造を捉えられそうなものであるが，現状，物理学的な背景を持った最適化・ベイズ推論に基づいた手法の方が性能が良いようである．4
3 圧縮センシング
3.1 データ圧縮との違い
(Krzakala et al., 2015, p. 30) には大変魅力的な導入がなされている．
例えば JPEG 2000 ではデータ圧縮の技術が使われており，これは画像データが wavelet 基底表現では大変スパースなデータになることを用いている．
圧縮センシングは最初から観測がスパースであるとし，データの復元の代わりに観測がなんだったかを推定することを考える．この技術は MRI に応用される (Lustig et al., 2007) ことで有名である．
即ち，「データは正しい基底に関してはスパースになるはずである」という事前情報の下， \[ \vec{y}=G\vec{s},\qquad G\in M_{MN}(\mathbb{R}),M<N, \] という計画行列 \(G\) と低次元の観測のみを通じて，真の観測 \(\vec{s}\) を，ある行列 \(A\) に関して \[ \vec{s}=A\vec{x} \] として得られる \(\vec{x}\) はスパースになるという追加情報を通じて推定することが，圧縮センシングの問題になる．
3.2 スパースなデータの少数観測からの復元
\((A,\vec{x})\) の組について最適化を行い，最もスパースな \(\vec{x}\) を決定するという問題 \[ \operatorname*{argmin}_{\vec{x}}\|\vec{y}-F\vec{x}\|_{\ell_0} \] は NP-困難であるが，統計力学の方法により \(M/N=:\rho\) を一定にした比例的高次元極限 \(N\to\infty\) に関する漸近論が得られる．
結論として，\(\vec{s}\) の長さ \(N\) に対して，実際は \[ M=O\left(N^{1/4}\log^{5/2}N\right) \] の固定された (nonadaptive) 観測があれば十分な復元が可能であるという．
具体的には，特定の条件の下では，（適切な基底に関する）スパースな \(l^p\)-展開係数が大きい順に \(n\) 個わかれば，\(l^2\)-誤差が \(O(n^{1/2-1/p})\) のオーダーに従う (Pinkus, 1985)．
そして，\(M=O\left(n\log N\right)\) の観測があれば，この \(n\) 個の重要な係数が十分な精度で復元できる，というカラクリである．そのためのアルゴリズムには Basis Persuit (Chen et al., 1998) という凸最適化を通じた極めて効率的なアルゴリズムが利用可能である．
3.3 \(\ell_p\)-ノルム最適化の方法
代わりに \[ \operatorname*{argmin}_{\vec{x}}\|\vec{y}-F\vec{x}\|_{\ell_p} \] という最適化問題を考えるとこれは凸最適化問題であるため，突然線型計画法により多項式時間で解ける問題に変貌する．
\(p=1\) と取ることが実用上スパースな結果をうまく出してくれるようである．実際，一定の条件の下では \(p=0\) の場合（本当に解きたかった問題）の解に一致する (Candes and Tao, 2006)．
こうして LASSO 様のアルゴリズムにたどり着くが，このアプローチでは観測数の比 \(\rho\) は真の限界より大きく取る必要がある．一方で，後述の信念伝搬の方法だと復号限界ギリギリまで通用することが多い．
3.4 スピングラス対応
観測モデルが線型 Gauss であるとするとき，ランダムに観測されたスパース符号の復元の問題は，事前分布 \(p(x)=(1-\rho)\delta(x)+\rho\phi(x)\) が定める事後分布 \[ p(x|F,y)=\frac{1}{Z}\prod_{i=1}^N\biggr((1-\rho)\delta(x_i)+\rho\phi(x_i)\biggl)\prod_{\mu=1}^M\frac{1}{\sqrt{2\pi\Delta_\mu}}\exp\left(-\frac{\lvert y_\mu-F_{\mu-}\cdot x\rvert^2}{2\Delta_\mu}\right) \] を通じた Bayes 推論として定式化できる．するとこれは負の対数尤度 \[ H(x)=-\sum_{i=1}^N\log\biggr((1-\rho)\delta(x_i)+\rho\phi(x_i)\biggl)+\sum_{\mu=1}^M\frac{1}{2\Delta_\mu}\lvert y_\mu-F_{\mu-}\cdot x\rvert^2 \] が定める Boltzmann 分布と見ることができるが，これは長距離無秩序を持つスピングラス系の Hamiltonian となっており，この系の基底状態の特定として問題が捉え直せる．
3.5 平均場変分ベイズによる復号（一般論）
一般の事後分布が \(p(x|y)\,\propto\,p(y|x)p(x)\) の形で与えられているとし，分布 \(q\) に関するこの系の自由エネルギーは \[ \mathcal{L}=[\log p(y|x)p(x)]_q-\int_\Omega q(x)\log q(x)\,dx=[-E(y,x)]_q-\int_\Omega q(x)\log q(x)\,dx \tag{1}\] と表せる．ただし，\(E:=-\log p(y|x)p(x)\) は内部エネルギーとした．
すると，この \(\mathcal{L}\) を最小化する \(q\) が平衡分布となるわけである．これで自由エネルギーに関する変分原理に問題が定式化し直された．
平均場近似とは，この \(q\) に関して \[ q(x)=\prod_{i=1}^Nq_i(x_i) \] という積の形を仮定した場合の表現 \[ \mathcal{L}_{\mathrm{MF}}=[-E(y,x)]_q-\sum_{i=1}^N\int q_i(x_i)\log q_i(x_i)\,dx \] に関して最小化を考えることをいう．
実はこの場合の \(\mathcal{L}_{\mathrm{MF}}\) は KL-乖離度となっており，停留条件は \[ q_i(x_i)=\frac{e^{[-E(y,x)]_q}}{Z_i} \] で得られ，これは平均場方程式と呼ばれるものである．5
このことから，\(q\) が一般の分布族の場合でも，KL 乖離度を最小化することによって推論をするという変分推論の指針が示唆される．
3.6 平均場変分ベイズによる復号（圧縮センシングの場合）
この場合，\(R_i,\Sigma^2_i\) を用いれば \[ Q_i(x_i)=\frac{1}{Z_i}p(x_i)\frac{1}{\sqrt{2\pi\Sigma_i^2}}e^{-\frac{(x_i-R)^2}{2\Sigma_i^2}} \] というように，事前分布と Gauss 分布の積の形を持ち，\(\Sigma_i^2,R_i\) の値も順次決定される．
するとその結果を用いれば逐次的に計算することができる．実際，これを \[ \sum_{\mu}F^2_{\mu i}=1 \] の仮定の下で行ったものが，Iterative Threshouding (Maleki and Donoho, 2010) である．
3.7 信念伝搬による復号
式 (1) は \[ \mathcal{L}=[\log p(y|x)]_q-\operatorname{KL}(q,p) \] とも表せる．第一項は人工的な分布 \(q\) に関する平均であるから計算できるものとすると，\(\operatorname{KL}(q,p)\) も，平均場近似の下では，各 \(x_i\) は独立と仮定しているから各項ごとに計算できる： \[ \operatorname{KL}(q_i,p_i)=-\log\widetilde{Z}_i-\frac{c_i+(a_i-R_i)^2}{2\Sigma^2_i}. \] 右辺はノルム，\(q_i\) による \(x_i\) の平均と分散によって計算できる．
実際，メッセージ \[ m_{\mu\to i}(x_i)=\frac{1}{Z^{\mu\to i}}\int\prod_{j\ne i}dx_je^{-\frac{1}{2\Delta_\mu}\left(\sum_{j\ne i}F_{\mu j}x_j+F_{\mu i}x_i-y_\mu\right)^2}\prod_{j\ne i}m_{j\to\mu}(x_j) \] \[ m_{i\to\mu}(x_i)=\frac{1}{Z^{i\to\mu}}\biggr((1-\rho)\delta(x_i)+\rho\phi(x_i)\biggl)\prod_{\gamma\ne\mu}m_{\gamma\to i}(x_i) \] の伝播により計算できるところを，\(N\to\infty\) の極限では，密度 \(m_{\mu\to i},m_{i\to\mu}\) そのものではなく平均 \(a_{i\to\mu}\) と分散 \(v_{i\to\mu}\)，その局所的な近似 \(a_i,v_i\)（「信念」）を伝えるのみで近似計算が可能である．
その結果，とりわけ圧縮センシングの問題に関しては，正確性とスピードの面で極めて効率的な信念伝搬法が得られる (Kschischang et al., 2001), (Yedidia et al., 2003)．アルゴリズムが停止したと判断されたのち，\(a_i\) によって真の信号 \(x_i^*\) の推定とし，\(v_i\) は推定の不確実性を表す．\(N\to\infty\) の極限において一致性を持つ．
BP はあるクラスの分布（混合 Gauss など）に従う観測 \(x\) について，理論的な復号限界ピッタリまで復元可能である．
3.8 近似メッセージ伝搬 (AMP)
この信念伝搬方では \(2M\times N\) のメッセージが送信されるが，ある仮定の下では \(N+M\) のメッセージで十分である (David L. Donoho et al., 2009)．
これは物理的には TAP 方程式 (Thouless et al., 1977) に基づく消息であり，\(N\to\infty\) の極限で信念伝搬法に一致する．
その後，AMP は必ずしもスパースな観測に限らずとも，一般の設定における復号（特に \(M\)-推定）においても有用であることが報告されている (D. Donoho and Montanari, 2016)．
低ランク行列の推定問題に対する変種 (Antenucci, Krzakala, et al., 2019) も提案されている．
4 低ランク行列分解
低ランク行列を分解するという推定問題においては，情報理論的な復号限界の手前に，困難相 (hard phase) が生じ，そこではメッセージ伝搬法が無効になる．
低ランク行列分解における困難相は，AMP が適用可能な範囲にも出現し，そこでは AMP の効率が大きく落ちてしまうという (Antenucci, Franz, et al., 2019)．
この困難相は物理的にはスピングラス相に由来するものであるから，そのことを踏まえて AMP を修正した ASP (Approximate Survey Propagation) が提案されている (Antenucci, Krzakala, et al., 2019) が，それでも性能は改善しないという！
この困難相というものは，何か背後の物理的な理由（はたまたもっと深遠な理由？）によって，効率的なアルゴリズムは存在しないということが本当に保証されている領域なのかもしれない：
This result provides further evidence that the hard phase in inference problems is algorithmically impenetrable for some deep computational reasons that remain to be uncovered. (Antenucci, Franz, et al., 2019)
5 参考文献
5.1 (Krzakala et al., 2015)
Lecture note なので歩み寄りやすい．Bayes 推論とスピングラス系の対応，これを解くための信念伝搬法について非常に良い入門になる．
第四章でグラフのコミュニティ抽出，第五章で圧縮センシングについて扱っている．クラスタリングはとかく物理的な背景が深いことが，相転移に関する詳細な解析を通じてよくわかる．クラスタリング面白すぎる．\(k\)-means，EM アルゴリズム，Gibbs サンプラーはいずれもクラスタリングに使用可能で，いずれも物理を背景に持ち，拡散過程を極限に持ち，それぞれに相転移があるのかもしれない．
5.2 (Zdeborová and Krzakala, 2016)
In physics the origins of these methods trace back to the Bethe approximation (Bethe and Bragg, 1935), the work of Peierls (Peierls and Bragg, 1936) and the Thouless–Anderson–Palmer (TAP) equations (Thouless et al., 1977). Interestingly, these equations can be seen as “improved” variational mean-field equations, and this is at the roots of many perturbation approaches (such as (Plefka, 1982), (Georges and Yedidia, 1991) or (Kikuchi, 1951)). (Zdeborová and Krzakala, 2016, p. 471)
5.3 (Fortunato, 2010) 共同体抽出のレビュー
グラフによるモデリングとその共同体抽出が，社会学・疫学・生化学で活躍している様子がよくわかり，大変なモチベーションが上がるようなレビューである．
(村田剛志, 2009) では community detection は コミュニティ抽出 としている．
5.4 (D. L. Donoho, 2006) 圧縮センシングのレビュー
荘厳すぎる完成された理論の感じよ．
(池田思朗, 2012) が非常にわかりやすい入門的解説を日本語で与えている．
5.5 (Yedidia et al., 2003)
BP アルゴリズムに関する詳細な解説．特に Bethe 近似や変分推論との関係を強調しており，BP アルゴリズムは自由エネルギーに対する Bethe 近似の停留点に収束することが示されている．
なお BP アルゴリズムは (池田思朗 et al., 2004) は 確率伝搬法，(Kabashima, 2005) は 要約伝搬法 と訳している．
References
Footnotes
(Krzakala et al., 2015, p. 21) 4.1節，(Fortunato, 2010, p. 124) 9.2節．↩︎
(Krzakala et al., 2015, p. 29) “In fact previous methods like spectral methods are not so good as the algorithm proposed in this section.”↩︎
同様の事実を (池田思朗 et al., 2004, p. 396) 定理１は情報幾何学の言葉で述べている．↩︎