EMアルゴリズムから導くRichardson-Lucy法

## はじめに

本記事では、尤度関数がポアソン分布に従うときのEMアルゴリズムとRL法（Richardson-Lucy法）の関係について説明します。

RL法は比較的シンプルな画像デコンボリューション法として知られています。一方で、EMアルゴリズムから導出しようとすると、多くの概念を整理する必要があり、筆者自身、試行錯誤を重ねました。

本記事では、RL法のEMアルゴリズムによる導出を自分なりに整理した内容を共有します。まだ理解しきれていない部分もあるため、もし間違いや改善点があれば、ぜひコメントで教えていただけると幸いです。

EMアルゴリズムのRL法の導出には、以下の論文を参考にしました。

[Shepp, L. A., & Vardi, Y. (1982). Maximum Likelihood Reconstruction for Emission Tomography.](https://courses.grainger.illinois.edu/npre435/fa2021/notes/shepp_vardi.pdf)

## Richardson-Lucy法とは

RL法は、ポアソン統計に基づく画像復元手法の一つで、特に天文学や医療画像の分野で利用されています。RL法はEMアルゴリズムの枠組みの中で導出されます。つまり、反復的な方法を用いる手法であり、**最尤解に収束**することが知られています。

RL法の詳細な解説については、以下の記事が参考になります。

https://qiita.com/yusuke_s_yusuke/items/df45b1210573408a6fda

## EMアルゴリズムとは

EM（Expectation-Maximization）アルゴリズムは、**隠れた（未観測の）変数を含む確率モデルの最尤推定を行うための手法** です。特に、データの一部が観測できない場合や、直接パラメータを推定するのが難しい場合に有効な方法として知られています。

EMアルゴリズムは以下の2つのステップを交互に実行することで、推定値を改善していきます。

1. **Eステップ（期待値計算）**: 現在の推定値をもとに、隠れた変数（潜在変数）の期待値を計算します。
2. **Mステップ（最大化）**: Eステップで得られた期待値を使って、パラメータを更新します。

この2つのステップを繰り返すことで、最尤推定が収束していきます。

EMアルゴリズムの詳細な解説については、以下の記事が参考になります。

https://qiita.com/kenmatsu4/items/59ea3e5dfa3d4c161efb



## EMアルゴリズムの適用に必要な要素
ここでは、観測された光子カウントデータから真の画像を推定する問題を対象とします。
EMアルゴリズムを適用するために、以下の3つの変数を設定します。

- **推定パラメータ**: 真の画像 $W(\xi)$
- **観測量**: 観測画像 $H(x)$
- **潜在変数**: 実際に放射された光子数 $w(\xi)$

ここで、$W(\xi)$ は位置 $\xi$ における真の放射強度であり、$H(x)$ は観測位置 $x$ における観測データです。

ポアソン分布の性質より、放射光子数の期待値は以下のように表されます。

```math
E[w(\xi)] = W(\xi)
```

また、位置 $\xi$ から放射された光子のうち、観測位置 $x$ で検出された数を $w(x|\xi)$ とすると、その期待値は以下のようになります。

```math
E[w(x|\xi)] = W(\xi) P(x|\xi)
```

ここで、$P(x|\xi)$ は位置 $\xi$ から放射された光子が観測位置 $x$ に到達する確率を表します。

:::note
**一般的な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]$ を次のように求めます。

```math
\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)$ は、周辺のどこかで観測されるとすると次の条件を満たします。

```math
w(\xi) = \int w(x|\xi) dx
```

この両辺の期待値を取ると、

```math
E[w(\xi)|H, W] = \int E[H(x|\xi)|H, W] dx
```

先ほど求めた $E[H(x|\xi)|H, W]$ を代入すると、

```math
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)$と置くと次のように表せます。

```math
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ステップで推定した隠れ変数の期待値を、そのまま新しい推定値として採用することで最尤推定と一致し、以下のように表されます。

```math
W_{\text{new}}(\xi) = E[w(\xi) | H, W_{\text{old}}]
```

このため、一般的なEMアルゴリズムにおけるMステップの「最尤推定の最適化処理」が不要となり、よりシンプルな更新則が得られます。

結果的に
```math
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アルゴリズムを完全に理解できているわけではありませんが、これからも色々な手法に触れながら、また記事にまとめていきたいと思います。