はじめに
本記事では、尤度関数がポアソン分布に従うときのEMアルゴリズムとRL法（Richardson-Lucy法）の関係について説明します。
RL法は比較的シンプルな画像デコンボリューション法として知られています。一方で、EMアルゴリズムから導出しようとすると、多くの概念を整理する必要があり、筆者自身、試行錯誤を重ねました。
本記事では、RL法のEMアルゴリズムによる導出を自分なりに整理した内容を共有します。まだ理解しきれていない部分もあるため、もし間違いや改善点があれば、ぜひコメントで教えていただけると幸いです。
EMアルゴリズムのRL法の導出には、以下の論文を参考にしました。
Shepp, L. A., & Vardi, Y. (1982). Maximum Likelihood Reconstruction for Emission Tomography.
Richardson-Lucy法とは
RL法は、ポアソン統計に基づく画像復元手法の一つで、特に天文学や医療画像の分野で利用されています。RL法はEMアルゴリズムの枠組みの中で導出されます。つまり、反復的な方法を用いる手法であり、最尤解に収束することが知られています。
RL法の詳細な解説については、以下の記事が参考になります。
EMアルゴリズムとは
EM（Expectation-Maximization）アルゴリズムは、隠れた（未観測の）変数を含む確率モデルの最尤推定を行うための手法 です。特に、データの一部が観測できない場合や、直接パラメータを推定するのが難しい場合に有効な方法として知られています。
EMアルゴリズムは以下の2つのステップを交互に実行することで、推定値を改善していきます。
- Eステップ（期待値計算）: 現在の推定値をもとに、隠れた変数（潜在変数）の期待値を計算します。
- Mステップ（最大化）: Eステップで得られた期待値を使って、パラメータを更新します。
この2つのステップを繰り返すことで、最尤推定が収束していきます。
EMアルゴリズムの詳細な解説については、以下の記事が参考になります。
EMアルゴリズムの適用に必要な要素
ここでは、観測された光子カウントデータから真の画像を推定する問題を対象とします。
EMアルゴリズムを適用するために、以下の3つの変数を設定します。
- 推定パラメータ: 真の画像 $W(\xi)$
- 観測量: 観測画像 $H(x)$
- 潜在変数: 実際に放射された光子数 $w(\xi)$
ここで、$W(\xi)$ は位置 $\xi$ における真の放射強度であり、$H(x)$ は観測位置 $x$ における観測データです。
ポアソン分布の性質より、放射光子数の期待値は以下のように表されます。
E[w(\xi)] = W(\xi)
また、位置 $\xi$ から放射された光子のうち、観測位置 $x$ で検出された数を $w(x|\xi)$ とすると、その期待値は以下のようになります。
E[w(x|\xi)] = W(\xi) P(x|\xi)
ここで、$P(x|\xi)$ は位置 $\xi$ から放射された光子が観測位置 $x$ に到達する確率を表します。
一般的なEMアルゴリズムの問題設定との違い
EMアルゴリズムでは、潜在変数の事後分布を推定し、それを用いてパラメータを更新するために、KLダイバージェンスなどを考慮することがあります。しかし、RL法では潜在変数 $w(\xi)$ の期待値が $W(\xi)$ に直接一致するため、複雑な確率分布の推定が不要になります。
その結果、Eステップでは単純に期待値を求め、Mステップではその期待値に基づいてスケーリングするだけのシンプルな更新式が導かれます。次章では、具体的にEステップとMステップの計算方法を説明します。
Eステップ（期待値の計算）
Eステップでは、観測データ $H(x)$ のもとで、各位置 $\xi$ における放射光子数 $w(\xi)$ の期待値を求めます。
各位置 $\xi$ で放射される光子数 $w(\xi)$ は独立なポアソン分布に従い、各光子が観測位置 $x$ で検出される確率は $P(x|\xi)$ で決まります。
このとき、ポアソン分布の加法性により、観測データ $H(x)$ もポアソン分布に従い、異なるピクセル $\xi$ からの寄与が線形に分解できるという性質を持ちます。
したがって、各位置 $\xi$ から放射された光子 $W(\xi)$ が観測位置 $x$ で検出される確率 $P(x|\xi)$ を用いて、条件付き期待値 $E[w(x|\xi)|H, W]$ を次のように求めます。
\begin{align}
E[w(x|\xi)|H, W] &= E[w(x|\xi)|H(x), W]\\
&= H(x) \frac{W(\xi) P(x|\xi)}{\int W(\xi) P(x|\xi) d\xi}
\end{align}
ここで、「$w(x|\xi)$ の寄与は、$H(x)$だけでなく、他の観測点 $H(x')$ によって影響を受けるのではないか？」という疑問が生じるかもしれません。しかし、観測の性質上、各 $w(x|\xi)$ は確率的に $P(x|\xi)$ に基づいて決まり、他の $H(x')$ の値とは独立に決まります。
また、放射光子数 $w(\xi)$ は、周辺のどこかで観測されるとすると次の条件を満たします。
w(\xi) = \int w(x|\xi) dx
この両辺の期待値を取ると、
E[w(\xi)|H, W] = \int E[H(x|\xi)|H, W] dx
先ほど求めた $E[H(x|\xi)|H, W]$ を代入すると、
E[w(\xi)|H, W] = \int H(x) \frac{W(\xi) P(x|\xi)}{\int W(\xi) P(x|\xi) d\xi} dx
このようにして、各位置 $\xi$ における放射光子数の期待値を計算できます。
推定パラメータを$w_{\text{old}}(\xi)$と置くと次のように表せます。
E[w_{\text{old}}(\xi)|H, W] = \int H(x) \frac{w_{\text{old}}(\xi) P(x|\xi)}{\int w_{\text{old}}(\xi) P(x|\xi) d\xi} dx
Mステップ（期待値の最大化）
Mステップでは、Eステップで得た期待値を用いて、推定パラメータ $W(\xi)$ を更新します。
一般に、EMアルゴリズムでは Mステップで「Eステップで計算した潜在変数の期待値を用いて、完全データの対数尤度を最大化する」という処理を行います。
しかし、ポアソン分布の場合、最尤推定値は期待値と一致する という特性があります。これは、ポアソン分布の対数尤度関数を最大化する際に、最適な $W(\xi)$ が $w(\xi)$ の期待値と一致するためです。
したがって、Eステップで推定した隠れ変数の期待値を、そのまま新しい推定値として採用することで最尤推定と一致し、以下のように表されます。
W_{\text{new}}(\xi) = E[w(\xi) | H, W_{\text{old}}]
このため、一般的なEMアルゴリズムにおけるMステップの「最尤推定の最適化処理」が不要となり、よりシンプルな更新則が得られます。
結果的に
W_{\text{new}}(\xi) = W_{\text{old}}(\xi) \int \frac{H(x)}{\int W_{\text{old}}(\xi) P(x|\xi) d\xi} P(x|\xi) dx
この式は、RL法の更新式と一致します。したがって、EMアルゴリズムからRL法が導かれることがわかります。
まとめ
本記事では、EMアルゴリズムを用いることでRL法の更新式を自然に導けることを説明しました。EMアルゴリズムとRL法が密接に関係し、数学的な裏付けのもとで定式化できるのはとても興味深いですね。
ちなみに、RL法は1972年、1974年に考案され、EMアルゴリズムは1977年に汎用的に導入されました。これを考えると、人類は数学的な理論が確立される前から、このようなアルゴリズムを自然と活用していたのかもしれません。それを知ると、改めて数学の面白さを感じます。
まだまだEMアルゴリズムを完全に理解できているわけではありませんが、これからも色々な手法に触れながら、また記事にまとめていきたいと思います。
