情報論的学習理論テクニカルレポート 2009
Technical Report on Information-Based Induction Sciences 2009 (IBIS2009)

確率伝搬法による確率的画像処理における統計的性能評価

Statistical Performance Analysis in Probabilistic Image Processing byBelief Propagation
片岡 駿∗
Shun Kataoka

安田 宗樹 ∗
Muneki Yasuda

田中 和之 ∗
Kazuyuki Tanaka

Abstract: We propose new schemes for evaluating statistical performances of image
restoration based on Bayesian statistics. The schemes are constructed by means of belief propagation, which is a powerful approximate methods in probabilistic inference. Theyare reduced to problems involving the solution of simultaneous integral equations for distributions of messages in belief propagation. Our schemes are proposed for two cases. Oneof them is for the case in which the original image is given explicitly. The other one isfor the case in which the prior probability for the original image is given. We show somenumerical results of the statistical performance analysis for the probabilistic restoration ofbinary images.
Keywords: Bayesian statistics, Markov random fields, belief propagation, probabilistic
information processing, statistical mechanical informatics, Bayesian network

1

まえがき
確率的画像処理 (Probabilistic Image Processing) はベイズ統計 (Bayesian Statistics) にもとづいて確率モデル化が行われ, データの不確定性を確率モデルという形で系統的に取り扱うことのできる理論的枠組みとして知られている. 従来, 画像処理に対する確率モデルは大規模なものとなり, 計算困難の問題がつきまとうことから,現実の応用という点では, それほど注目されてきたわけではなかった. 最近, 確率伝搬法 (Belief Propagation:BP) とよばれる近似計算技法が提案されることで, 計算困難を回避することが比較的容易となり, 急速にその応用への道が拡がりつつある [1, 2, 3].
このような中で, 確率的画像処理の性能を統計的に評価する方法に対する要請が高まりつつある. 統計的性能評価は確率的に定式化されたシステムの性能を多くのサンプルに対する数値実験を通して得られた結果に対する標本平均を通して評価するのではなく, 解析的にその性能を計算してしまうものである. 確率的画像処理に対する統計的性能評価は事前確率分布をガウシアングラフィ∗ 東北大学大学院情報科学研究科, 980-8579 仙台市青葉区荒巻字青葉 6-3-09, tel. 022-795-5885, e-mail {xkataoka,muneki,kazu}@smapip.is.tohoku.ac.jp,
Graduate School of Information Sciences, Tohoku University,63-09,Aramaki-Aza-Aoba, Aobaku , Sendai 980-8579

カルモデルなどの可解モデルに仮定した場合に対してはいくつかの研究成果がすでに報告されている [4, 5]. しかし, 一般の確率モデル, 特に離散確率場に対してはまだ具体的な方法は提案されていない.
誤り訂正符号における低密度パリティ検査符号 (LowDensity Parity Check: LDPC) 符号, 移動体通信における符号分割多重接続 (Code Division Multiple Access:CDMA) 復調方式に対しての統計的性能評価法はレプリカ法 (Replica Method) とよばれる統計力学の方法をもとにして提案されている [6, 7, 8]. しかしレプリカ法が適用できるのは完全グラフ, 疎結合ランダムグラフなどのある特殊な構造をもつグラフ上の確率モデルに限定される. 確率的画像処理に対しては原画像の事前確率分
布を完全グラフ上で定義した場合についてはレプリカ法による統計的性能評価についての研究が報告されている[6, 9, 10, 11]. しかし, 実際の画像処理の応用において取り扱かわれるグラフは正方格子となる場合が多く, レプリカ法の適用範囲からは外れてしまう.
本研究では確率伝搬法をもちいた確率的画像処理に対する統計的性能評価法を提案する. 第 2 節では本論文で取り扱う確率的画像処理の定式化をベイズ統計の枠組みにもとづいて与える. 第 3 節では確率的画像処理における確率伝搬法の定式化とそれを拡張した統計的性能評価法の提案を行う. 第 4 節では第 3 節で提案した評価法についてのいくつかの数値実験の結果を与える. 第 5 節はまとめである.

画像 ⃗h を推定するためにもちいる事後確率分布は
()∑∑1exp αhi hj + βgi hi P ⃗h | ⃗g , α, β =Zh{i,j}∈E

確率的画像処理2

本節では, 確率的画像処理のベイズの公式にもとづく定式化を与え, 定式化した確率モデルによる画像修復のi∈V(3))(と導かれる.Zh は事後確率分布 P ⃗h | ⃗g , α, β の規格化定数である.

統計的性能を測るための枠組みを与える. 取り扱う原画
像は白黒の 2 値画像に, ノイズは加法的白色ガウスノイ2.2ズに限定して説明する.

2.1

最大周辺事後確率推定

本論文では式 (3) の条件付確率分布からの推定には各画素 i ごとに周辺事後確率分布を最大化する階調値をそ

画像のモデル化の画素の推定値 hi とする最大事後確率推定 (Maximum

原画像を白黒の 2 値画像とし, 各画素は正方格子上に配列しているとする. 各画素にラベル i を割り当てることとし, すべての画素の集合を V = {1, 2, · · ·, |V |} とあらわすこととする. また, すべての最近接画素対 {i, j} のPosterior Marginal: MPM) 推定を採用する. 画素 i の原画像の推定値 hi に対する周辺事後確率分布は)∑ (Pi (hi | ⃗g , α, β) =P ⃗h | ⃗g , α, β

(4)

⃗h\hi

集合を E とあらわす. 画素 i の階調値に対する確率変数を fi とし, i 番目の画素が白の状態をとるとき fi = +1,である. 本論文では 2 値画像を扱っているため, この周黒の状態をとるとき fi = −1 とあらわすことにする. 原
画像の確率場は f⃗ = (f1 , f2 , · · ·, f) により表現される.

辺事後確率分布を最大化する階調値 h∗i は hi の期待値

原画像に対する事前確率分布は次のように仮定する.
 J ∑()1P f⃗ | J =(fi − fj )2(1)
exp − 2
Zf{i,j}∈E

ここで J は正値のみをとるハイパパラメータであり, Zfは規格化定数である.
劣化画像は原画像において各画素 i ごとに独立にガウスノイズが fi に加えられることで生成されたものであるとし, 劣化画像の各画素 i の階調値を gi とあらわすことにする. すると原画像 f⃗ が与えられたときの劣化画像⃗g = (g1 , g2 , · · ·, g|V | ) が生成される条件付き確率分布は()
P ⃗g | f⃗, σ{}()|V |
11 ∑2= √exp − 2(gi − fi )(2)2σ2πσ 2i∈VTとあらわされる. ここで σ 2 は加法的白色ガウスノイズの分散である.
式 (1) および式 (2) をベイズの公式に代入することで劣化画像 ⃗g が与えられたときの原画像 f⃗ に対する事後確率分布を導出することができる. このとき原画像と区別するために原画像の推定画像を ⃗h = (h1 , h2 , · · ·, h|V | )であらわし, 推定者にとって未知となるパラメータ J とσ −2 をそれぞれ α と β であらわし, 劣化画像 ⃗g から推定

⟨hi | ⃗g , α, β⟩ をもちいて()h∗i = arg max Pi (hi | ⃗g , α, β) = sgn ⟨hi | ⃗g , α, β⟩hi(5)により与えられ，劣化画像 ⃗g からの原画像の推定画像はh∗ を成分とする ⃗h∗ で与えられる．
i2.3

最大周辺事後確率推定の統計的評価MPM 推定による画像の修復度をはかる尺度として,本論文では原画像 f⃗ と推定画像 ⃗h∗ の間の規格化されたハミング距離を採用する. これは)1 ∑(Dh (f⃗, ⃗h∗ ) ≡1 − δfi ,h∗i
|V |i∈V
11 ∑= −fi sgn(⟨hi |⃗g , α, β⟩)
2 2|V |

(6)i∈Vによりあらわされる.

MPM 推定における原画像と推定画像の間の修復度の統計的性能は以下のように Dh (f⃗, ⃗h∗ ) の統計平均により定義される. まず Dh (f⃗, ⃗h∗ ) を原画像 f⃗ と劣化画像 ⃗g に関して平均した量() ()∑∫[Dh ]f,g ≡d⃗g Dh (f⃗, ⃗h∗ )P ⃗g | f⃗, σ P f⃗ | Jf⃗∫1 ∑∑
1= −d⃗g fi sgn(⟨hi |⃗g , α, β⟩)
2 2|V |i∈V f⃗() ()× P ⃗g | f⃗, σ P f⃗ | J(7)は事前分布 (1) に従い生成されるような全ての画像に関式 (6) の確率伝搬法による近似表式は式 (9) によりする修復度の期待値をあらわしており, 本論文であつか

Dh (f⃗, ⃗h∗ )う確率的画像修復モデルの性能評価となる. 式 (7) にお∑∑
1
1λij ≃ −fi sgn βgi +2 2|V |いて原画像はその分布を除いて特定されていないことに注意. 他方, Dh (f⃗, ⃗h∗ ) を劣化画像 ⃗g に関して平均した量[Dh ]g∫()≡ d⃗g Dh (f⃗, ⃗h∗ )P ⃗g | f⃗, σ∫()
11 ∑= −d⃗g fi sgn(⟨hi |⃗g , α, β⟩)P ⃗g | f⃗, σ2 2|V |

i∈V

(11)

j∈∂iとあらわされる. すなわち固定点方程式 (10) の解を式(11) に代入することにより, 原画像 f⃗ と観測された劣化
画像 ⃗g から確率伝搬法により得られた推定画像 ⃗h∗ の間のハミング距離が計算されることとなる.

i∈V

(8)

3.2

はある特定の（分布をもたない）画像 f⃗ に対する画像修復に関して期待される性能をあらわしている.
式 (8)(の [D)h ]g は式 (7) の [Dh ]f g において画像の事前分布 P f⃗ | J をある特定の画像 f⃗ に関するデルタ関数とすることにより得られることに注意.

メッセージの分布関数による最大周辺事
後確率推定の統計的評価

ここでは式 (7) および式 (8) により定義された原画像と推定画像の間のハミング距離の統計的評価の確率伝搬法による近似表式を与える. 原画像やノイズによりこの確率モデルはランダムネスをもつ. したがってそれらのランダムネスによる分布を考慮した確率伝搬法を考える

確率伝搬法と統計的性能評価

3

本節では, 式 (7) および式 (8) により定義された原画像と推定画像の間の修復度の統計的評価を近似的に計算する方法を確率伝搬法をもとに与える. ここで述べる方法は, 統計力学において確率伝搬法に対応する方法をス
ピングラスモデルに適用した研究 [12] を確率的画像処理に拡張することで与えられる.1

3.1のが自然である. そのためには確率伝搬法におけるメッセージ λij の分布関数 Pij1 (λij ) というものが導入され,その分布関数をもちいて上記のハミング距離の統計的評価があらわされる形となる.
まず原画像が未知の場合, 原画像が未知であるがそれを生成した事前確率分布 (1) のハイパパラメータ J と原画像が与えられた時の劣化画像を生成する事後確率分
布 (2) のハイパパラメータ σ, 未知パラメータ α, β を既知とすると確率伝搬法におけるメッセージ λij の分布関

確率伝搬法

確率伝搬法とは与えられた確率分布の周辺確率分布を計算する近似アルゴリズムである. 式 (3) の確率分布に対する周辺確率分布 (4) の確率伝搬法における近似表式¯{}
はメッセージとよばれる量 λij , λji ¯{i, j}∈E をもちいて次のようにあらわされる.
 
∑
1λij  hi(9)exp βgi +Pi (hi | ⃗g , α, β) ≃Zi

数 Pij1 (λij ) は式 (10) より次の積分方程式の解として与えられる.

∫Pij1 (λij ) =dgj Pj (gj )[{∏k∈∂j\{i}1dλjk Pjk(λjk )× δ λij − tanh−1 tanh(α) tanh βgj +∑}]λjk k∈∂j\{i}(12)j∈∂i{ ¯}
ここで Zi は規格化定数であり, ∂i = j ¯{i, j}∈E は画素 i のすべての隣接画素からなる集合である. メッセー¯{}
ジ λij , λji ¯{i, j}∈E は以下の固定点方程式の解として与えられる.




∑
λij = tanh−1 tanh(α) tanh βgj +λjk  (10)k∈∂j\{i}

ここで δ[· · · ] は Dirac のデルタ関数であり, Pi (gi ) は次の式で定義される.
∫) ()∏∑ (Pi (gi ) ≡dgj
P ⃗g | f⃗, σ P f⃗ | J (13)j∈V \{i}f⃗

ただし Pi (gi ) を厳密に求めることはやはり難しいので,Pi (gi ) も確率伝搬法をもちいて計算する必要があるが,ここではメッセージの分布を考慮する必要がないので通常の確率伝搬法の範囲内で計算可能である.

1 確率伝搬法は統計力学においてベーテ近似 (Bethe Approxima-tion) とよばれているものと等価である．

積分方程式 (12) の解として得られるメッセージに対する分布関数をもちいて式 (7) の [Dh ]f,g は近似的に次のようにあらわされる.

[Dh ]f,g ≃11−2 2|V |∑∑∫i∈Vfidgi Pi (gi )× Pij1 (λij )fi sgn βgi +(α, β) = (0.465, 1) において最小値 0.0466 をとっている.
これは α = J, β = σ −2 で最小値となることを意味して∏dλijj∈∂i∑λij 

おいて不等式により与えられている最適修復条件と一致

(14)

j∈∂i

原画像 f⃗ が特定されている場合のメッセージの分布関数 Pij2 (λij ) は次の積分方程式の解として与えられる.

∫Pij2 (λij ) =

∏dgj

[{}
(gj − fj )2exp −2σ 22πσ 2{× δ λij − tanh

図 1 の結果は画像の事前確率分布のハイパパラメータ{J, σ} と一致するモデルのパラメータ {α, β} をもちいたときが最適であることを示しており, これは一見自明な結果のように考えられる. しかしながらここでの統計にはこの最適性は自明ではない (最適点のシフトが起こる可能性がある). したがってここでの結果は上記の最1−1する.

評価は確率伝搬法による “近似的” 評価であるため一般

2(λjk )
dλjk Pjkk∈∂j\{i}×√おり, 文献 [6, 9] における 2 値画像の確率的画像処理にtanh(α) tanh βgj +∑}]

適性に関して本論文で扱う画像処理モデルに対して確率伝搬法が統計的に良好な結果を与えることがわかる.

λjk k∈∂j\{i}(15)
この積分方程式の解として得られるメッセージに対する分布関数をもちいて式 (8) の [Dh ]g は次のようにあらわされる.

∫∏
11 ∑1
√
[Dh ]g ≃ −dgi
dλij Pij2 (λij )22 2|V |2πσi∈Vj∈∂i{}∑(gi − fi )2× fi sgn βgi +λij exp
(16)
2σ 2

最小値
[Dh]f,g0.16
0.120.08
0.04
1.801.40.20.4
α 0.6 0.8

10.6

β1 0.2j∈∂i

数値実験4本節では第 3 節で与えられた確率伝搬法をもちいた

図 1: J = 0.465 の事前確率分布 (1) に対して分散 σ 2 = 1の加法的白色ガウスノイズのもとでの統計的性能 [Dh ]f gの (α, β) 依存性.

Dh の統計的性能評価の近似式 (14), (16) に α, β の値を具体的に与えて数値実験を行う. ここで劣化画像を生成0.16する条件付確率分布 (2) のハイパパラメータ σ の値を 10.144.1

原画像が未知の場合

0.12[Dh]f,gとした. なお式 (12)，(15) の積分方程式はモンテカルロ積分を用いて評価した．

0.1

0.08

原画像が特定されていないが, 原画像の事前確率分布

0.06(1) の J の値が与えられ, その事前確率分布に従って生0.04

成されるすべての原画像に対する統計的性能 [Dh ]f,g に00.20.4 α0.60.81ついての数値実験を与える. 統計的性能 [Dh ]f,g は前節

図 2: J = 0.465 の事前確率分布 (1) に対して分散 σ 2 = 1の積分方程式 (12) を数値的に解き, 得られたメッセージの加法的白色ガウスノイズのもとで β = 1 と固定しての分布関数 Pij1 (λij ) を式 (14) に代入することで計算さ場合の統計的性能 [Dh ]f g の α-依存性.

れる. さまざまの α, β の値に対して得られる [Dh ]f,g の値を図 1 に与える. ここで原画像に対する事前確率分布(1) のハイパパラメータ J の値を 0.465 としている. また図 2 は図 1 における β = 1 の [Dh ]f,g のふるまいを抜き出したものである. 図 1, 図 2 において [Dh ]f,g は原画像が既知の場合4.20.160.14的に固定した場合の, 加法的白色ガウスノイズに対する0.12

統計的性能 [Dh ]g についての数値実験を与える. 統計的

0.1[Dh]g

原画像を図 3(画像サイズ 128 × 128) として 1 つ具体

性能 [Dh ]g は前節の積分方程式 (15) を数値的に解き, 得0.08られたメッセージの分布関数 Pij2 (λij ) を式 (16) に代入0.06することで計算される. さまざまの α, β の値に対して得0.04られる [Dh ]g の値を図 4 に与える. 図 4 から β = 1 の部

0.02

00.20.4 α

0.60.81分を抜き出したものを図 5 に与える. [Dh ]g は α = 0.56付近で最小値をとり α = 0.56 での [Dh ]g の値は 0.0311

図 5: 図 3 に対して分散 σ 2 = 1 の加法的白色ガウスノである. この結果から図 3 の画像を修復する際の最適パイズのもとでの統計的性能 [Dh ]g の β = 1 の場合の αラメータは (α, β) = (0.56, 1) 付近であることがわかる.

依存性.

第 2 節で述べたように [Dh ]g は [Dh ]f,g において原画像 f⃗ をある画像に固定することにより得られる. 実際

5まとめの画像処理の性能評価を行う場合には原画像 f⃗ は具体

本論文では原画像が既知である場合と未知である場合的に与えた上で行われるので, [Dh ]f,g よりむしろ [Dh ]gに対して確率伝搬法のメッセージの分布関数を考えるこの方が実際には重要となる場合が多いと考えられる. 図とで確率的画像処理における統計的性能評価を行った.

4 はハミング距離の意味で統計的に最適な点が存在することを示しており, これは確率的画像処理においてパラその結果, ループのないグラフ構造で定式化された確率メータ推定の重要性を主張するための動機となる. 本提した近似的統計評価法においても確率的画像処理の統計案手法は各画像について最適点が存在するか否かの判定的性能評価が可能であることを示した. 今回は簡単のたツールとしての発展が期待される.

めに 2 階調の場合を考えたが, 本論文の評価法は階調値伝搬法による統計評価法をループのある正方格子に適用が多値の場合にも拡張が可能であり, これは今後の課題である．また期待値最大化 (Expectation Maximization:EM) アルゴリズムによるハイパパラメータ推定法に対する統計的評価 [5, 10, 11] への本論文の提案法の拡張も今後に残された課題である.

謝辞
図 3: 第 4.2 節の数値実験でもちいる原画像

本研究の一部は文部科学省科学研究費補助金 (No.21700247,No.18079002) およびグローバル COE プログラム “Center of Education and Research for Information Electronics Systems” の補助を得て行われたものである.

最小値[Dh]g0.14
0.1
0.06
0.02

参考文献
[1] A. S. Willsky: “Multiresolution Markov models for1.801.4
0.20.4α10.6

0.80.6
1β0.2

図 4: 図 3 に対して分散 σ 2 = 1 の加法的白色ガウスノ
イズのもとでの統計的性能 [Dh ]g の (α, β) 依存性.

signal and image processing”, Proceedings of theIEEE, vol.90, no.8, pp.1396-1458 (2002).
[2] K. Tanaka, “Statistical-mechanical approach to image processing”, Journal of Physics A: Mathematical and General, vol.35, no.37, pp.R81-R150(2002).

[3] 田中和之, 確率モデルによる画像処理技術入門, 森北出版 (2006)[4] K. Tanaka and J. Inoue, “Maximum likelihood hyperparameter estimation for solvable Markov random field model in image restoration”, IEICE
Transactions on Information and Systems, vol.E85D, no.3, pp.546-557 (2002).
[5] K. Tanaka and D. M. Titterington, “Statisticaltrajectory of approximate EM algorithm for probabilistic image processing”, Journal of Physics
A: Mathematical and Theoretical, vol.40, no.37,
pp.11285-11300 (2007).
[6] H. Nishimori, Statistical Physics of Spin Glass andInformation Processing, Oxford University Press(2001).
[7] Y. Kabashima and D. Saad, “Statistical mechanics of low density parity check codes”, Journal ofPhysics: Mathematical and General, vol.37, no.6,
pp.R1-R43 (2004).
[8] T. Tanaka, “A statistical-mechanics approach tolarge-system analysis of CDMA multiuser detectors” IEEE Transactions on Information Theory,
vol.48, no.11, pp.2888-2910 (2002).
[9] H. Nishimori and K. Y. M. Wong, “Statistical mechanics of image restoration and error-correcting
codes”, Physical Review E, vol.60, no.1, pp.132144 (1999).
[10] J. Inoue and K. Tanaka, “Dynamics of the maximum marginal likelihood hyperparameter estimation in image restoration: Gradient descent versus
expectation and maximization algorithm”, Physical Review E, vol.65, no.1, Article No.016125(2002).
[11] J. Inoue and K. Tanaka, “Mean field theory of EM
algorithm for Bayesian gray scale image restoration”, Journal of Physics A: Mathematical andGeneral, vol.36, no.43, pp.10997-11010 (2003).
[12] T. Morita, “Spin-glass and ferromagnetic phasesof the random-bond Ising model on the Bethelattice”, Physica A, vol.125, no.2/3, pp.321-343(1984).