統計数理（2004）
第 52 巻 第 2 号 393–405
c 2004 統計数理研究所
［研究詳解］

情報幾何学に基づく確率伝搬法の解析
1

2

池田 思朗 ・田中 利幸 ・甘利 俊一

3（受付 2004 年 2 月 27 日；改訂 2004 年 6 月 29 日）

要旨1980 年代後半 Pearl が提案した確率伝搬法は，大規模なグラフィカルモデルに対する確率推論のための計算手法である．同等の手法は統計物理学，統計学，誤り訂正符号の復号法などにも存在し，広く用いられている．確率伝搬法は木の構造のグラフに対してはグラフの大きさに比例した計算量で厳密解が得られる．しかしループを持つグラフに対しては繰り返し計算の収
束性，および得られた結果の近似精度ともに理論的には十分理解されていなかった．一方で確率伝搬法は実用上有効な手法であり，その性質を理論的に明らかにすることは重要である．本
研究では情報幾何学に基づく枠組みにより確率伝搬法を表現し，収束性や近似精度を議論する．

キーワード： 確率伝搬法，情報幾何学，グラフィカルモデル．
1.

はじめに

グラフィカルモデル（Lauritzen and Spiegelhalter（1988）, Jordan（1999））では，複数の確率変数の同時分布をグラフによって表現する．グラフによって表現された複数の確率変数の一部のみが観測されたとき，その条件付確率分布から観測されていない確率変数の値を推論する問
題を考える．この問題は確率推論と呼ばれ，人工知能，統計物理学，情報理論など様々な分野で重要である．
Pearl（1988）はこの問題に対し，確率伝搬法（Belief Propagation）と呼ばれる簡便な繰り返し計算法を提案した．確率伝搬法は条件付き分布から各確率変数の周辺分布を求め，確率推論を得る．木のグラフに対しては収束が保証され，グラフの辺の数に比例した計算量で正しい推論
結果が得られることが分かっているが，グラフがループを持つ場合には必ずしも収束しないこと，収束した場合でも得られる解が一般に近似解となることが知られている．
同等の計算手法は他の分野でも広く用いられている．統計物理のベーテ（Bethe）
近似（Kabashima
and Saad（1998）），また，低密度パリティ検査符号（Gallager（1962）, MacKay（1999）），ター
ボ符号（Berrou et al.（1993）, McEliece et al.（1998））といった誤り訂正符号の復号法などは確率伝搬法と等しいことが知られている．統計物理や符号理論で扱う問題に対するグラフィカル
モデルは一般にはループを持ち，収束性や近似精度といった問題がある．
確率推論の問題は MCMC などの手法によっても解が得られる．その場合ループがあっても構わないが，精度の良い解を得るために必要な計算量は多い．前に挙げた誤り訂正符号ではグラフの確率変数の数が通常数百から千程度であり，実時間で復号するには計算量の少ない確率
1
2
3

統計数理研究所：〒106–8569 東京都港区南麻布 4–6–7
東京都立大学大学院 工学研究科：〒192–0397 東京都八王子市南大沢 1–1
理化学研究所 脳科学総合研究センター：〒351–0198 埼玉県和光市広沢 2–1

394

統計数理 第 52 巻

第 2 号 2004

伝搬法が適している．
これまで我々は情報幾何学（甘利・長岡（1993）, Amari and Nagaoka（2000））に基づき確率伝搬法を表現し，解の安定性，近似精度といった問題を扱ってきた（Ikeda et al.（2002）， 池田 他
（2002），Tanaka et al.（2002），Ikeda et al.（2003, 2004a, 2004b））．本稿ではこれまでの結果をまとめ，新たに得られた 3 次の摂動展開に基づく近似精度の評価について示す．
2.

確率推論の問題と情報幾何的枠組み

2.1 問題の表現
x = (x1 , . . . , xn )T を観測できない確率変数，y = (y1 , . . . , ym )T を観測された確率変数とする．本稿では簡単のため xi が 2 値変数，特に xi ∈ {−1, +1} の場合を考える．多値変数への拡張は簡単であり，連続値への拡張も場合によっては可能である（Ikeda et al.（2003, 2004b））．
確率推論の問題は y の条件付きでの x の分布 q(x|y) から（簡単のため以下では q(x|y) をq(x) とかくことにする），x に関する推論を得ることである．1 つの方法は q(x) を最大にする x（MAP 推論: maximum a posteriori）を用いることである．MAP 推論は推論結果が真のx と異なる確率を最小にするが，探索空間が n とともに指数関数的に増える．ここでは別の推論，周辺事後確率分布の最大化（MPM 推論: maximization of the posterior marginals）を考える．q(x) の周辺分布を q(xi ), i = 1, . . . , n とするとき MPM 推論では各成分の推論結果をq(xi = +1) ≥ q(xi = −1) ならば x̂i = +1，それ以外の場合は x̂i = −1 とする．この結果各 xiが誤って推論される確率は最小となる．ηi を q(x) による xi の期待値とする，
ηi = Eq [xi ] =xi q(xi ) .
xi

MPM 推論は x̂i = sgn ηi とする．これは周辺分布の積ni=1 q(xi ) あるいは x の期待値

η = Eq [x] ,

が分かればただちに計算できる．仮に n とともに指数的に増加しない，簡単な計算によって ηあるいは ni=1 q(xi ) が計算できれば MPM 推論は有効な推論手法となる．本稿で扱う確率伝
搬法は簡単な計算によって ni=1 q(xi ) の近似を得ようというものである．
今 xi は 2 値なので，全ての q(x) に対し ln q(x) は {xi } の高々 n 次関数として表現できる．
グラフィカルモデルで表現できる問題の多くは確率変数間の相互関係が限られており，低次の関数で表現できる場合が多い．そこで ln q(x) を次のように表現する



L

(2.1)



ln q(x) = h · x +cr (x) − ψq ,r=1

ここで h · x = i hi xi は xi の線形項，cr (x), r = 1, . . . , L，は単項式でも多項式でもよいが高次項を表し，ψq は規格化定数の対数である．cr (x) はクリーク関数とも呼ばれる．ボルツマン
マシンや古典的なスピングラスのモデルでは cr (x) は xi の二次関数である．
cr (x) = wij xi xj ,ここで r は辺の番号を表すもので，xi と xj 間の辺を示す．以下ではこのモデルをボルツマン
マシンと呼ぶ（図 1）．
一般のグラフィカルモデル，特に無向グラフでは q(x) を以下のようにクリーク関数の積として定義することが多い．

nq(x) =1
φi (xi )φr (xr ) ,Zq i=1r∈C395

情報幾何学に基づく確率伝搬法の解析

図 1.

ボルツマンマシン．

ここで C はクリークの集合である．(2.1) 式の表現と φi (xi ), φr (xr ) とは次の関係がある．
hi =1 φi (xi = +1)ln,
2 φi (xi = −1)

cr (x) = ln φr (xr ) ,ψq = ln Zq .

ボルツマンマシンでは φr (xr ) は φr (xi , xj ) となる．本論文では cr (x), r = 1, . . . , L は一次独立であるとする．なお，本稿では無向グラフのみを扱うが全ての有向グラフは無向グラフで表現できることから一般性は失っていない（Lauritzen and Spiegelhalter（1988））．
2.2 情報幾何学と多様体
本節では情報幾何学に基づく議論のための準備を行う．
まず，確率伝搬法を考えるために必要な確率分布全体の族 S を考える．
(2.2)

S = {p(x; θ, v)

p(x; θ, v) = exp[h · x + θ · x + v · c(x) − ψ(θ, v)], θ ∈ n , v ∈ L } ,(θ, v) は自然パラメータであり，θ = (θ1 , . . . , θn )T , v = (v1 , . . . , vL )T , c(x) = (c1 (x), . . . , cL (x))T ,v · c(x) = L
r=1 vr cr (x) である．S は定義より指数型分布族である．q(x) が S に含まれることは θ = 0, v = 1L = (1, . . . , 1)T とおけば簡単に確かめられる．次に v = 0 となる S の部分多様体を M0 と呼ぶことにする．
M0 = {p0 (x; θ) = exp[h · x + θ · x − ψ0 (θ)] | θ ∈ n } .

M0 の分布は各成分が独立であり，自然パラメータは θ である．さらに，各成分が独立な分布は全て M0 に含まれることから q(x) の周辺分布の積 ni=1 q(xi ) も M0 に含まれる．すなわち
ni=1 q(xi ) に対応する M0 の座標 θ を見つければ MPM 推論が得られる．
次に S の中の e-平坦，m-平坦な部分多様体を定義する．
e-平坦な部分多様体: 部分多様体 M ⊂S は次の r(x; t) が全ての t ∈ [0, 1], q(x), p(x) ∈ M に対して M に含まれるとき e-平坦である．



ln r(x; t) = (1 − t) ln q(x) + t ln p(x) + c(t) ,c(t) は規格化定数である．e-平坦な部分多様体が 1 次元の曲線のとき，特に e-測地線と呼ぶ．
{r(x; t) | t ∈ [0, 1]} は p(x) と q(x) を結ぶ e-測地線である．
m-平坦な部分多様体: 部分多様体 M ⊂S は次の r(x; t) が全ての t ∈ [0, 1], q(x), p(x) ∈ M
に対して M に含まれるとき m-平坦である．
r(x; t) = (1 − t)q(x) + tp(x) .
m-平坦な部分多様体が 1 次元の曲線のとき，特に m-測地線と呼ぶ．{r(x; t) | t ∈ [0, 1]} は p(x)と q(x) を結ぶ m-測地線である．

396

統計数理 第 52 巻

第 2 号 2004

定義よりただちに指数型分布族が e-平坦であることがわかる．したがって S と M0 は e-平坦である．次に m-射影を定義する．本稿では e-平坦な多様体への m-射影を考えれば十分であるので，次のように定義する．
Definition 1. M を S の e-平坦な部分多様体とし，q(x)∈S とする．q(x) から M 上の点
p(x) への m-測地線が Fisher 情報行列を Riemann 計量として M と直交する点を q(x) から
M への m-射影と呼ぶ．
Theorem 1. M が S の e-平坦な部分多様体であるとき q(x)∈S から M への m-射影は唯一である．今 q(x)∈S から M への m-射影を ΠM ◦q(x) とかくと，ΠM ◦q(x) は M においてq(x) からの Kullback-Leibler（KL）divergence を最小にする点である．すなわち(2.3)

ΠM ◦q(x) = argmin D[q(x); p(x)] = argmin
p(x)∈M

p(x)∈Mq(x) ln
x

q(x)
.
p(x)

ここで D[·; ·] は KL-divergence である．
2.3 MPM 推論
本節では MPM 推論が q(x) から M0 への m-射影と等しいことを示す．(2.3) 式の定義から
q(x) から M0 への m-射影がパラメータ θ ∗ に対応するならば
p0 (x; θ ∗ ) = ΠM0 ◦ q(x)とかける．これより M0 への m-射影により求まるパラメータ θ を次のように書くことにする．
θ ∗ = πM0 ◦ q(x) = argmin D[q(x); p0 (x; θ)] .
θD[q(x); p0 (x; θ)] を θ で微分しxq(x) − ∂θ ψ0 (θ∗ ) = 0 .

(2.4)
x

∂θ は θ による微分を表す．指数型分布族の定義から(2.5)

exp(h · x + θ · x) =∂θ ψ0 (θ) = ∂θ lnxxp0 (x; θ) .
x

ここで M0 の期待値パラメータ η0 (θ) を次のように定義する．
η0 (θ) =(2.6)

xp0 (x; θ) = ∂θ ψ0 (θ) .
x

M0 の 2 つの座標系 θ と η 0 とは 1 対 1 に対応するので，(2.4), (2.5), (2.6) 式から m-射影がq(x) の周辺化と等しく，MPM 推論と等価であることがわかる．
3.

確率伝搬法の情報幾何

3.1 確率伝搬法の情報幾何的表現
本節では確率伝搬法の情報幾何的な表現を与える．グラフィカルモデルにおける確率伝搬法はメッセージと呼ばれる変数の更新規則を定義するのが通常である（Pearl（1988）, Weiss（2000））．
まず，ボルツマンマシンに対する確率伝搬法を Yedidia et al.（2001）にしたがって定義する．確率伝搬法ではメッセージ mji (xi ) を次の式で更新し，全てが収束した後，得られたメッセージ

397

情報幾何学に基づく確率伝搬法の解析m∗ji (xi ) を用いてビリーフ bi (xi ) を求める．
mt+1
ji (xi ) =(3.1)

1
φj (xj )φij (xi , xj )Z xj1bi (xi ) =  πi (xi )Z

mtkj (xj )k∈N (j)\im∗ki (xi ) .
k∈N (i)N (j) は j 番目の節につながっている節の集合，Z ，Z  はそれぞれ mt+1
ji (xi ) および bi (xi ) のxi に関する和を 1 に規格化する．木のグラフではビリーフ bi (xi ) は q(x) の xi に関する周辺
分布 q(xi ) と一致するが，ループのある場合には一般に一致せず，近似となる．
以下ではこの確率伝搬法の情報幾何学に基づく表現を与える．確率伝搬法では q(x) と p0 (x; θ)
さらに次式の pr (x; ζ r ), r = 1, . . . , L, を用いる．
pr (x; ζ r ) = exp[h · x + cr (x) + ζ r · x − ψr (ζ r )] ,

ζ r ∈ n ,r = 1, . . . , L .

pr (x; ζ r ) はクリーク関数を 1 つだけ含んだ確率分布である．cr (x) は一般に {xi } の高次の関数であるが，{cr (x)} の全てを含む場合に比べて 1 つのみを含む場合は扱い易い．他のクリーク関数は ζ r · x で代替されている．pr (x; ζ r ) は指数型分布族であり，以下で定義される Mr はe-平坦である．
Mr = {pr (x; ζ r ) | ζ r ∈ n } ,r = 1, . . . , L .

自然パラメータは ζ r である．Mr の期待値パラメータを ηr (ζ r ) と定める．
(3.2)

η r (ζ r ) = ∂ζr ψr (ζ r ) =xpr (x; ζ r ) ,r = 1, . . . , L .

xη r (ζ r ) を求めるための計算量は pr (x; ζ r ) の周辺化と同等である．このためには cr (x) に含まれる xi の全ての組み合わせについて計算しなければいけないが，この計算は可能だとする．全ての r に対して η r (ζ r ) を求めるほうが q(x) の周辺化を直接行うより圧倒的に計算量は少ない．なお，ζ r から η r (ζ r ) を求めるのは簡単だが，一般に逆は簡単でないことを注意しておく．
確率伝搬法では q(x) の代わりに pr (x; ζ r ) の周辺化を全ての辺の数だけ行い，全体としてq(x) の近似を行う．つまり pr (x; ζ r ), r = 1, . . . , L によって各 cr (x) の影響を表現し，全ての影
響を θ にまとめ，p0 (x; θ) によって i q(xi ) を近似する．その際，繰り返し計算によって {ζ r }および θ を更新する．
我々は，誤り訂正符号で用いられる特殊なグラフィカルモデルに対しては池田 他（2002）および Ikeda et al.（2004a）で，一般のグラフに対しては Ikeda et al.（2003, 2004b）で確率伝搬法の情報幾何的な表現を得ている（Richardson（2000）も同様の枠組みを与えている）．ここではその結果を示す．
p0 (x; θt ) および pr (x; ζ tr ) を時刻 t におけるそれぞれ M0 , Mr , r = 1, . . . , L での q(x) の近似とする．確率伝搬法は ζ tr および θ t を以下のように更新する．



確率伝搬法の情報幾何的表現（1） 初期値を t = 0, ξ tr = 0, ζ tr = 0, r = 1, . . . , L とする．
（2） t を１つずつ増加させ ξt+1
, r = 1, . . . , L を次のように更新する．
rξ t+1= πM0 ◦pr (x; ζ tr ) − ζ tr .
r

(3.3)

（3） θt+1 と ζ t+1
を以下のように更新する．
r
(3.4)

ζ t+1
=
r
r  =rξt+1
r ,θ t+1 =

ξt+1
=
r
r

1
L−1

ζ t+1
.
r
r

398

統計数理 第 52 巻

第 2 号 2004

（4） 2 と 3 を {ξ tr } が収束するまで繰り返す．
θ t = r ξ tr , θt = ξtr + ζ tr は常に成り立つ．確率伝搬法が収束した時点でのパラメータを θ ∗ ,{ζ ∗r }, {ξ∗r } とする．直感的に確率伝搬法を理解するには q(x), p0 (x; θ∗ ), pr (x; ζ ∗r ), r = 1, . . . , Lに関する以下の関係を見るのが良い．
q(x) = exp[h · x + c1 (x) + · · · + cr (x) + · · · + cL (x) − ψq ]p0 (x; θ ∗ ) = exp[h · x + ξ∗1 · x + · · · + ξ ∗r · x + · · · + ξ∗L · x − ψ0 (θ ∗ )]p1 (x; ζ ∗1 ) = exp[h · x + c1 (x) + · · · + ξ∗r · x + · · · + ξ ∗L · x − ψ1 (ζ ∗1 )]..
.
pr (x; ζ ∗r ) = exp[h · x + ξ∗1 · x + · · · + cr (x) + · · · + ξ∗L · x − ψr (ζ ∗r )]..
.
pL (x; ζ ∗L ) = exp[h · x + ξ∗1 · x + · · · + ξ ∗r · x + · · · + cL (x) − ψL (ζ ∗L )] .

この式から確率伝搬法では各 cr (x) を ξ∗r · x によって表現していることがわかる．
情報幾何的表現で用いた ξ r と (3.1) 式で用いたメッセージとの関係をボルツマンマシンの場合に示す．辺 r が i と j の節を結ぶものとすると，
ξr,i =



1 mji (xi = +1)ln,
2 mji (xi = −1)

ξr,j =

1 mij (xj = +1)ln,
2 mij (xj = −1)

ξr,k = 0 for k = i, j .

この関係を用いると (3.1) 式のメッセージの更新則と (3.3) 式が対応すること，p0 (x; θ ∗ ) がi bi (xi ) と等しいことが分る．
3.2

停留点の持つ性質

Theorem 2. 確率伝搬法の停留点は以下の 2 つの条件を満たす（Ikeda et al.（2004a））．
m-条件: θ∗ = πM0 ◦ pr (x; ζ ∗r ).
e-条件:

θ∗ =

L

1
ζ∗.
L − 1 r=1 r

確率伝搬法の停留点で m-条件が満たされていることは (3.3) 式と θ ∗ = ζ ∗r + ξ∗r から簡単に確かめられる．e-条件は (3.4) 式からただちに確かめられる．
情報幾何的な意味を明らかにするため，S の 2 つの部分多様体 M ∗ と E ∗ を定義する．
∗

M =
(3.5)

∗

E =




p(x)




p(x) ∈ S,

xp(x) =x
∗ t0p(x) = Cp0 (x; θ )


L

r=1




∗

∗




xp0 (x; θ ) = η 0 (θ ) ,x

pr (x; ζ ∗r )trL

tr = 1, tr ∈  ,r=0

C : 規格化定数.
∗

M は m-平坦な部分多様体であり M0 および Mr , r = 1, . . . , L と直交する．一方，E ∗ は e平坦な部分多様体である．また定義より，M ∗ は p0 (x; θ∗ ) を，E ∗ は p0 (x; θ∗0 ) と pr (x; ζ ∗r ),r = 1, . . . , L を含む．
前の 2 つの条件はこれらの部分多様体を用いて次のようにかける．
m-条件: M ∗ が pr (x; ζ ∗r ), r = 1, . . . , L を含む．

399

情報幾何学に基づく確率伝搬法の解析

e-条件: E ∗ が q(x) を含む．
m-条件が上のように書き直せることは定義より明らかである．e-条件については (3.5) 式で∗
t0 = −(L − 1), t1 = · · · = tL = 1 とおけば θ∗ = L
r=1 ζ r /(L − 1) と同値なことがわかる．
確率伝搬法が収束した点では e-条件と m-条件が同時に満たされるが，だからといって p0 (x; θ∗ )が真の周辺分布 ni=1 q(xi ) であるわけではない．これは M ∗ と E ∗ との間の差から生じる．M ∗と E ∗ はともに p0 (x; θ∗ ) と pr (x; ζ ∗r ) を含んでいるが，他の点については明らかではない．本
来ならば q(x) が M ∗ に含まれればよいのだが，計算量を減らすためにこれを E ∗ で置き換えるのが確率伝搬法である．同様の仕組みは平均場近似などでも用いられている．
特殊なのはグラフが木の構造の場合である．木のグラフでは確率伝搬法によって必ず収束し，正しい周辺分布が求まることが分っている（Pearl（1988））．これより次の結果が得られる．



Proposition 1. q(x) が木のグラフで表現できるとき q(x), p0 (x; θ∗ ), pr (x; ζ ∗r ), r = 1, . . . , Lが M ∗ と E ∗ に含まれる．



木のグラフの場合には p0 (x; θ∗ ) = i q(xi ) であり，確率伝搬法によって真の周辺分布が求まるが，ループのあるグラフの場合には一般には q(x) ∈/ M ∗ である．すなわち確率伝搬法の収束点は特殊な場合を除き，真の周辺分布を与えない．
3.3 停留点の安定性
ここでは線形近似に基づき，確率伝搬法の停留点の局所的な安定性を考える．停留点 {ζ ∗r }に摂動 ∆ζ r を加え ζ r = ζ ∗r + ∆ζ r とする．確率伝搬法のステップ（2）の結果を ξr = ξ ∗r + ∆ξrとすると ∆ξr は次のように表される．
∆ξr = I0 (θ ∗ )−1 Ir (ζ ∗r )∆ζ r − ∆ζ r =(I0 (θ∗ )−1 Ir (ζ ∗r ) − En )∆ζ r .

ここで I0 (θ∗ ) は p0 (x; θ) の，Ir (ζ ∗r ) は pr (x; ζ r ) の Fisher 情報量行列であり En は n 次の単位行列である．確率伝搬法の計算を一ステップ行った後のパラメータを ζ r = ζ ∗r + ∆ζ r とかくと，∆ζ r =

L

(I0 (θ∗ )−1 Ir (ζ ∗r ) − En )∆ζ r .

r  =r

この結果から次の条件が成り立つとき，確率伝搬法の停留点が安定であることがわかる．
Theorem 3. （Ikeda et al.（2004a））停留点の周りでの確率伝搬法の性質は次のように近似
できる．




ここで

T =





∆ζ 1
..
.
∆ζ L

=T




∆ζ 1
..
.
∆ζ L

O

I0−1 I2 − En

I0−1 I1 − En..
.
I0−1 I1 − En

O

···
..

···

,

.

···

I0−1 IL − En
..
.
..
.
O

,

である．ただし I0 = I0 (θ ∗ ), Ir = Ir (ζ ∗r ) とした．T の全ての固有値 λi , i = 1, . . . , nL が |λi | < 1
を満たすとき，停留点は安定である．

400

統計数理 第 52 巻

第 2 号 2004

上の結果は全てのパラメータを同時に更新する場合の停留点の条件である．実際にはパラメータを一つずつ更新する場合やランダムに更新するなど様々な方法があり，それぞれで安定性は多少異なる．いずれの場合も同様な解析が可能である．
3.4 解の近似精度
本節では摂動展開に基づく近似精度の評価，すなわち確率伝搬法で求まった結果と真の周辺分布との差について述べる．2 次の摂動展開の結果については Ikeda et al.（2004a）にあるが，まずその結果を一般のグラフィカルモデルに対して導き，その後ボルツマンマシンに対する 3
次の摂動展開の結果について述べる．
以下では (2.2) 式の p(x; θ, v) を用いて議論する．定義より明らかに
p0 (x; θ) = p(x; θ, 0) ,pr (x; ζ r ) = p(x; ζ r , er ) ,q(x) = p(x; 0, 1L )である．ただしer = (0, · · · , 0, 1, 0, · · · , 0)T ,
↑
r

 

1L = (1, · · · , 1)T =
L

L

er .
r=1

さらに p(x; θ, v) の期待値パラメータを次のように定義する
η(θ, v) = ∂θ ψ(θ, v) =xp(x; θ, v) .
x

近似精度の解析のため p(x; θ, v) を M ∗ に拘束する．すなわち
η(θ, v) = η(θ∗ , 0) = η(θ∗ )が常に満たされるようにする．このとき θ は v の関数である．必要があれば θ(v) とかく．以
下では v が 0 から 1L まで変化するにつれて θ が θ ∗ からどのように変化するかを評価する．
まず θ(v) の 2 次までの摂動展開の結果を示す．M ∗ に拘束されていることから
0=

(3.6)

d
∂η ∂θ
∂η
η(θ, v) =
+
.
dv
∂θ ∂v
∂v

∂η/∂θ および ∂η/∂v は p(x; θ, v) の Fisher 情報量行列の一部である．これ以降 i, j, k は θ のr, s, t は v の成分を表すことにする．Gθθ = (∂η/∂θ), Gθv = (∂η/∂v) と定める．ここで任意のv と θ = θ(v) に対して Gθθ = I0 (θ) であることに注意する．(3.6) 式から(3.7)

0 = I0 (θ)∂θ+ Gθv (θ) ,∂v∂θ= − I0−1 (θ)Gθv (θ) .
∂v

これにより θ(v) の v に関する 1 次の微分 ∂θ/∂v が得られた．また G̃θv を −∂θ/∂v と定義す
る．2 次の微分係数はd2
η(θ, v) = 0 ,
dvdv

より
(3.8)

I0 (θ)

∂2θ
∂θ ∂θ
∂θ
∂θ
= − Tθvv  − Tθθθ
− Tθθv  − Tθθv
∂v∂v ∂v ∂v ∂v
∂vとなる．ここでTθθθ =∂3ψ,∂θ∂θ∂θ

Tθθv =∂3ψ,∂θ∂θ∂v

Tθvv  =∂3ψ∂θ∂v∂v 401

情報幾何学に基づく確率伝搬法の解析である．∂θ/∂v = −G̃θv (θ) および ∂ 2 θ/∂v 2 を (θ, v) = (θ∗ , 0) で評価し，θ(v) をまずは 2 次のテーラー展開により v で近似する．(θ, v) = (θ∗ , 0) におけるオペレータ d/dv を以下のように
定義すれば
∂
∂
d
=B=
− G̃θv (θ ∗ )
.
dv
∂v
∂θ




∂ 2 θ/∂v∂v は∂2θ
= −I0 (θ∗ )−1 B 2 η(θ∗ ) ,∂v∂v v=0となる．簡単のため B 2 の (r, s) 成分を Brs = Br Bs とかくことにする．なお (d2 /dvdv)η(θ∗ ) = 0
であっても一般に B 2 η(θ∗ ) = 0 である．
θ(v) を (θ∗ , 0) の周りで v の 2 次までで近似した結果は以下の通りである．
θ(v) = θ ∗ +

(3.9)







∂θ
1
∂2θ
v + vTv + o(||v||3 )∂v v=02∂v∂v v=0 θ ∗ − G̃θv (θ∗ )v −1 T −1 ∗v I0 (θ )(B 2 η(θ∗ ))v .
2

ここまで m-条件を考えていたが，e-条件を考慮する必要がある．e- 条件はθ∗ = −(3.10)L(ζ ∗r − θ ∗ ) ,r=1とかける．次の分布を考える(3.11)

p(x; ζ r , er ) = exp[h · x + ζ r · x + cr (x) − ψ(ζ r , er )] .

ここで p(x; ζ r , er ) =1 = pr (x; ζ r ) である．p(x; ζ r , er ), r = 1, . . . , L が M ∗ に含まれるように拘束すると (3.9) 式の結果から ζ r − θ∗ は によって以下のように近似できる．
ζ r − θ∗  −G̃θv (θ∗ )er −1 −1 ∗I (θ )Brr η(θ ∗ ) 2 .
2 0

以下の議論では θ に対する cr (x) の影響は小さいと仮定する．この仮定の下では がある程度大きな値を取っても近似が成り立つことから → 1 とした点を考える．この意味での近似を≈ によって表すこととする． → 1 ととれば ζ r − θ ∗ は ζ ∗r − θ∗ となりζ ∗r − θ∗ ≈ −G̃θv (θ∗ )er −1 −1 ∗I (θ )Brr η(θ ∗ ) .
2 0(3.10) 式から θ ∗ は次式を満たす．
(3.12)

θ∗ = −L(ζ ∗r − θ∗ ) ≈ G̃θv (θ∗ )1L +r=1

1 −1 ∗I (θ )
2 0

Brr η(θ∗ ) .
r

ここで次の分布を考える．
(3.13)

p(x; u, 1L ) = exp[h · x + u · x + 1L · c(x) − ψ(u, 1L )] .

この分布を M ∗ に拘束し を 0 から 1 まで増やすと u は θ(1L ) となるとする．一般に θ(1L ) = 0であることから p(x; θ(1L ), 1L ) は q(x) と等しくない．(3.9) 式の結果から
θ(1L ) − θ ∗ ≈ −G̃θv (θ ∗ )1L −1 −1 ∗I (θ )
Brs η(θ∗ ) .
2 0
r,s

402

統計数理 第 52 巻

第 2 号 2004

(3.12) 式から
1
θ(1L ) ≈ − I0−1 (θ∗ )Brs η(θ∗ ) .
2
r=s(3.14)

期待値パラメータはテーラー展開に基づき(3.15)

η(0, 1L )  η(θ(1L ), 1L ) − ∇θ η(θ∗ )θ(1L ) ≈ η(θ∗ ) +1
Brs η(θ∗ )2 r=sと近似できる．η(0, 1L ) は q(x) による x の期待値である．したがって (3.15) 式は q(x) による x の期待値と確率伝搬法による結果との差を示している．
Theorem 4. 真の条件付確率 q(x) による x の期待値を ηM P M = η(0, 1L ) とし，確率伝搬法によって求まった期待値を η(θ∗ ) とする．2 次までの摂動展開の結果 ηM P M と η(θ ∗ ) の差は以下のように近似できる．
1(3.16)
ηM P M ≈ η(θ∗ ) +Brs η(θ∗ ) .
2 r=s

この定理に示された近似誤差は M ∗ の e-曲率と関係している（Ikeda et al.（2004a））．グラ
フィカルモデルが与えられればこの項を具体的に計算できる．具体的な形はモデルによって異なるが，cr (x) が {xi } の単項式の場合には次の定理が得られている (Tanaka et al.（2002））．



Theorem 5. cr (x), cr (x), r = r  が共通な xi を 2 つ以上含まないとき，摂動展開による 2
次の誤差項 r=s Brs η(θ ∗ )/2 は 0 となる．

ボルツマンマシンでは cr (x) = wij xi xj となるため，cr (x), cr (x), r = r  は共通の xi を持っていたとしても一つのみである．したがってこの近似誤差項は 0 となることがわかる．
ボルツマンマシンの誤差項を評価するためには少なくとも 3 次の摂動展開を行う必要がある．ここでは最近得た結果を示す．3 次までの近似は(3.17)

θ(v)  θ∗ +







∂2θ
∂3θ
∂θ
1
1
v + vTv+
vr vs vt∂v v=02∂v∂v v=0
6 r,s,t ∂vr ∂vs ∂vt v=0






 


から求まるわけだが，(3.14) 式の誤差項については 2 次までの項が 0 となることから
1
θ(1L ) 
6∂3θ
−∂vr ∂vs ∂vt v=0
r,s,t

r

∂3θ
∂vr3 v=0となる．∂ 3 θ/∂vr ∂vs ∂vt はd3
η(θ, v) = 0 ,dvdvdvから多少の計算をすれば求まる．ボルツマンマシンに対して具体的に計算をすると，次の結果を得る．
2θi (1L ) ≈ ηi∗

2

wij wjk wki (1 − ηj∗ )(1 − ηk∗ ) .
j=k

これは一番短いループのみが θi (1L ) に寄与することを示している（図 2）．

403

情報幾何学に基づく確率伝搬法の解析

図 2. 誤差項へ寄与する短いループ．

Theorem 6. 3 次までの摂動展開により，ボルツマンマシンに対する η M P M と η(θ∗ ) の差は以下のように近似できる．
22

ηM P M,i ≈ ηi∗ − ηi∗ (1 − ηi∗ )(3.18)

2

wij wjk wki (1 − ηj∗ )(1 − ηk∗ ) .
j=k

現在は 3 次までの項の結果しかないが，さらに高次の摂動展開はより長いループを反映した項が現われるものと考えられる．
(3.16), (3.18) 式の結果を用いると MPM 解の精度が向上できる期待がある．しかし，例えばボルツマンマシンの場合に計算してみると，必ずしも精度が向上するとも限らない．原因としてはいくつか考えられる．摂動展開の仮定としている θ への cr (x) の影響が小さいとするが，wij がある程度大きい場合にはこの近似が成りたたない，また，ボルツマンマシンなどではより高次の項の影響があり，3 次の摂動展開の項のみを考えることはあまり意味がない，などである．精度は確率伝搬法の一つの重要な問題である．今後，どうやって 2 次，3 次，あるいはより高次の摂動展開と精度向上を結びつけるかが課題の一つである．
4.

まとめ

確率伝搬法は統計物理学や誤り訂正符号など，応用上重要な手法である．本稿では確率伝搬法を理解するための情報幾何学的な枠組みについて述べた．
一般に確率伝搬法で重要な問題は収束性と近似精度である．我々は提案した枠組みに基づき，
局所的な安定性の条件を示し，近似誤差についても摂動展開から主要項を示した．近似誤差については 3 次までの漸近展開の結果では最も短いループが主要項を構成しているが，今後より高次の展開についても調べる必要がある．
確率伝搬法は様々な分野で独自の研究がなされており，関連するアルゴリズムや近似についても広く論じられている．情報幾何学による枠組みではそれらを全て表現し，等しく扱うことが可能である．我々は CCCP, TRP, GBP といった関連するアルゴリズムがこの枠組みで表現できることを示した（Ikeda et al.（2004b））．今後は各分野で得られた結果を理解し，新たなア
ルゴリズムの提案に結び付けたいと考えている．
謝

辞

査読者には有益な御指摘を頂きました．感謝致します．

404

統計数理 第 52 巻

第 2 号 2004

参 考 文 献甘利俊一，長岡浩司（1993）. 『情報幾何の方法』，岩波講座 応用数学［対象 12］，岩波書店，東京．
Amari, S. and Nagaoka, H.（2000）. Methods of Information Geometry, AMS and Oxford University
Press, Providence, Rhode Island.
Berrou, C., Glavieux, A. and Thitimajshima, P.（1993）. Near Shannon limit error-correcting coding
and decoding: Turbo-codes, Proceedings of IEEE International Conference on Communications, Geneva, Switzerland.
Gallager, R. G.（1962）. Low density parity check codes, IRE Transactions on Information Theory,
IT-8, 21–28.
池田思朗，田中利幸，甘利俊一（2002）. ターボ復号の情報幾何，電子情報通信学会論文誌， J85-D-II（5）,
758–765.
Ikeda, S., Tanaka, T. and Amari, S.（2002）. Information geometrical framework for analyzing belief
propagation decoder, Advances in Neural Information Processing Systems, Vol. 14（eds. T. G.
Dietterich, S. Becker and Z. Ghahramani）, 407–414, MIT Press, Cambridge, Massachusetts.
Ikeda, S., Tanaka, T. and Amari, S.（2003）. Stochastic reasoning, free energy and information geometry, Research Memorandum, No. 890, The Institute of Statistical Mathematics, Tokyo.
Ikeda, S., Tanaka, T. and Amari, S.（2004a）. Information geometry of turbo codes and low-density（6）, 1097–1114.
parity-check codes, IEEE Transactions on Information Theory, 50
Ikeda, S., Tanaka, T. and Amari, S.（2004b）. Stochastic reasoning, free energy and information geometry, Neural Computation, 16（9）, 1779–1810.
Jordan, M. I.（1999）. Learning in Graphical Models, MIT Press, Cambridge, Massachusetts.
Kabashima, Y. and Saad, D.（1998）. Belief propagation vs. TAP for decoding corrupted messages,
Europhysics Letters, 44（5）, 668–674.
Lauritzen, S. L. and Spiegelhalter, D. J.（1988）. Local computations with probabilities on graphicalstructures and their application to expert systems, Journal of the Royal Statistical Society B,
50, 157–224.
MacKay, D. J. C.（1999）. Good error-correcting codes based on very sparse matrices, IEEE Transactions on Information Theory, 45（2）, 399–431.
McEliece, R. J., MacKay, D. J. C. and Cheng, J.-F.（1998）. Turbo decoding as an instance of Pearl’s“belief propagation” algorithm, IEEE Journal on Selected Areas in Communications, 16（2）,
140–152.
Pearl, J.（1988）. Probabilistic Reasoning in Intelligent Systems: Networks of Plausible Inference,
Morgan Kaufmann, San Mateo, California.
Richardson, T. J.（2000）. The geometry of turbo-decoding dynamics, IEEE Transactions on Infor（1）, 9–23.
mation Theory, 46
Tanaka, T., Ikeda, S. and Amari, S.（2002）. Information-geometrical significance of sparsity in Gallager codes, Advances in Neural Information Processing Systems, Vol. 14（eds. T. G. Dietterich,
S. Becker and Z. Ghahramani）, 527–534, MIT Press, Cambridge, Massachusetts.
Weiss, Y.（2000）. Correctness of local probability propagation in graphical models with loops, Neural（1）, 1–41.
Computation, 12
Yedidia, J. S., Freeman, W. T. and Weiss, Y.（2001）. Bethe free energy, Kikuchi approximation,and belief propagation algorithms, Technical Report 2001-16, Mitsubishi Electric ResearchLaboratories, Cambridge, Massachusetts.

Proceedings of the Institute of Statistical Mathematics Vol. 52, No. 2, 393–405（2004)

405

Information Geometrical Framework to Analyze Belief Propagation
Algorithm
Shiro Ikeda1 , Toshiyuki Tanaka2 and Shun-ichi Amari3
1 The Institute of Statistical Mathematics

2 Department of Electronics and Information Engineering, Tokyo Metropolitan University
3 Brain Science Institute, RIKEN

Belief propagation (BP) is a universal method of stochastic reasoning. It gives exact
inference for stochastic models with tree interactions, and works well even if the models
have loopy interactions. Its performance has been analyzed separately in many fields,such as, AI, statistical physics, information theory, and information geometry. The present
paper provides a unified framework for understanding BP. The stability of BP is analyzed
from this framework, and its approximation accuracy is investigated.

Key words: Belief propagation, information geometry, graphical model.