MCMCの画像のノイズ除去への応用

adventar.org
この記事はKyoto Universuty Advent Calendar 2021の23日目の記事です.

元々は逆強化学習について記事を書こうと考えていたのですが,
時間が無かったので別の内容について雑な記事を書きました. 申し訳ございません.
授業か何かでやった内容が面白かったので記事にまとめました.

はじめに
背景というか定義
MCMC
- Metropolis-Hastings法
- Gibbs sampling
MCMCを応用した画像のノイズ除去
- 原理
- 結果と実装
雑な感想と近況報告

はじめに

詳細釣り合いと極限分布の関係を雑に雑に述べて, 次にMCMCの代表的な手法を説明する.
最後にMCMCの応用として画像のノイズ除去を述べる.

背景というか定義

概要しか書けていません. 申し訳無いです...

・離散時間マルコフ連鎖
定義 : 状態空間 $\mathbb{S}$ 上の確率過程 $\{X_n:n \in \mathbb{Z}\}$ が $\forall n \in \mathbb{N}$ と状態列 $j_0,...,j_n \in \mathbb{S}$ に対して,
$P(X_{n+1}=j | X_n=i, X_{n-1}=j_{n-1},...,X_0 = j_0) = P(X_{n+1}=j|X_n = i)$ を満たす時,
$\{X_n\}$ を状態空間 $\mathbb{S}$ 上の離散時間マルコフ連鎖と呼ぶ.

※ マルコフ性 : 次の状態 $X_{t+1}$ が現在の状態 $X_t$ にのみ依存する事

・推移確率行列
定義 : 行列 $\mathbf{P}=(P_{i,j})_{i,j \in \mathcal{S}} \geq \mathbf{O},\ \ P_{i,j}=P(X_{n+1}=j|X_n = i)$ をマルコフ連鎖 $\{X_n\}$ の推移確率行列と呼ぶ.
要は $i \rightarrow j$ に遷移する確率を要素 $P_{i,j}$ とする行列の事.

推移確率行列 $\mathbf{P}$ については以下が成立する事が知られている.
1. $\mathbf{P}$ の各行の和は $1$ となる. ( $\mathbf{Pe}=\mathbf{e}$ )
2. $\forall n, m \in \mathbb{Z}_+$ に対して, $\mathbf{P}^{n+m} = \mathbf{P}^n\mathbf{P}^m$ が成立

2つめはチャップマン=コルモゴロフ方程式と呼ばれる. これはマルコフ性から証明できる.

・既約(連結クラス)
定義 : 状態部分集合 $\mathbb{C \subset S}$ が既約とは, $\forall (i,j)\in \mathbb{C}^2$ に対して, $i$ から $j$ に到達可能という事.
雑にまとめるとどの状態間も行き来できるという事.

・再帰的
定義 : $\tau_j$ を状態 $j$ への初到達時刻と置く. $P(\tau_i < \infty | X_0 = i) = \sum_{n=1}^{\infty} P(\tau_i = n|X_0 = i)=1$ が成立する時,
状態 $i$ は再帰的であると言う. 任意の状態が再帰的な時, マルコフ連鎖は再帰的と言う.

※ 再帰的かつ $\mathbb{E}[\tau_i|X_0 = i$ < \infty]も成立する時, 状態 $i$ は正再帰的であると言う.
※ 規約なマルコフ連鎖に置いては, 「ある状態 $i$ が(正)再帰的 $\Rightarrow$ 全ての状態が(正)再帰的」が成立.
(互いの状態の到達可能性から示せる.)

・周期性
定義 : 状態 $i$ に対して, $\{ n \in \mathbb{N}\ |\ P(X_n = i|X_0 = i) = \mathbf{P}^{(n)}_{i,i} > 0\}$ の最大公約数 $d_i$ を状態 $i$ の周期と呼ぶ.

※ 周期 $d_i$ が $1$ に等しい時, 非周期的であると言う.
※ 互いに到達可能な状態は周期が同じになる.

・定常分布
定義 : $\mathbf{\pi P } = \mathbf{\pi}$ を満たす確率ベクトル $\mathbf{\pi} = (\pi_i)_{i \in \mathbb{S}}$ の事.

※ 常に存在するとは限らない.
※ 規約かつ再帰的なマルコフ連鎖は定常分布を持つ.
※ 定常分布と極限分布は異なる(定常分布に到達するとは限らない)

・極限分布
定義 : 任意の初期分布 $\mathbf{\pi}^(0)$ に対して, $\lim_{n \rightarrow \infty}\mathbf{\pi}^{(n)}$ が同じ確率ベクトル $\mathbf{p} = (p_s)_{s \in \mathbb{S}}$ に収束する時,
$\mathbf{p}$ を極限分布と呼ぶ. 推移確率行列がエルゴート的(規約かつ非周期的で再帰的な事)な時には極限分布と定常分布が一致する. ( $\lim_{n\rightarrow \infty}P^{(n)}_{i,j} = \pi_j,\ \ \ \forall (i,j) \in \mathbb{S}$ )

・詳細釣り合い
定義 : 推移確率行列 $\mathbf{P}$ に対して, $\mu_i P_{i,j} = \mu_j P_{j,i}$ なる方程式の事.
これが成立する時, $\sum_{i \in \mathbb{S}} \mu_i P_{i,j} = \sum_{i \in \mathbb{S}} \mu_j P_{j,i} = \mu_j$ が成立するので
$\mathbf{\mu}$ は $\mathbf{\mu P} = \mathbf{\mu}$ を満たす. (これは定常測度である事に注意)

これは, 詳細釣り合いを満たす様にマルコフ連鎖を構成してやれば目的の分布を定常分布に持つ様にできる事を意味している. MCMCはこれを利用している.

MCMC

求めたい分布を定常分布に持つマルコフ連鎖を用いて, 分布からのサンプルを生成する手法のこと.
十分に分布が収束するまでの期間(バーンイン)のサンプルデータは破棄する.

Metropolis-Hastings法

===================================
提案分布を $h$ とし, 近似したい分布の密度関数を $g$ とする. エルゴード性を仮定.
0. $\mathbf{X}^{(0)} := \mathbf{x}$ : 初期値
1. $\mathbf{V}^{(t)}|\mathbf{x}^{(t)} \sim h(\mathbf{v}^{(t)}|\mathbf{x}^{(t)})$
2. $a^{(t)}=\min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}$
3. $\left\{\begin{array}{lll} P(\mathbf{X}^{(t+1)} = \mathbf{V}^{(t)}) = a^{(t)} \\ P(\mathbf{X}^{(t+1)} = \mathbf{x}^{(t)}) = 1 - a^{(t)} \end{array}\right.$ として, $\mathbf{V}^{(t)}$ を移動先の点として受け入れるか確率的に決定する. -> 1に戻る.
===================================

詳細釣り合いの成立について,
1. $\mathbf{x}^{(t)} = \mathbf{x}^{(t+1)}$ の時, 自明
2. $\mathbf{x}^{(t)} \neq \mathbf{x}^{(t+1)}$ の時, $p(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) = \min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})$ と書ける.
なので,
$\begin{equation}\begin{split}p(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)}) &= \min\left\{1, \frac{h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})}{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)})g(\mathbf{X}^{(t)})} \right\}h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) g(\mathbf{X}^{(t)}) \\ &= \min\{h(\mathbf{V}^{(t)}|\mathbf{X}^{(t)}) g(\mathbf{X}^{(t)}),\ h(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)})\} \\ &= p(\mathbf{X}^{(t)}|\mathbf{V}^{(t)})g(\mathbf{V}^{(t)}) \end{split}\end{equation}$
が成立する.
以上より詳細釣り合いが示された. また, この時の分布は $\mathbf{X}^{(t)}\sim g(\mathbf{x})$ に従う.

Gibbs sampling

===================================
$\mathbf{x} = [(\mathbf{x})_1,...,\ (\mathbf{x})_m$ ], $(\mathbf{x})_{-i} := [(\mathbf{x})_1,...,(\mathbf{x})_{i-1},(\mathbf{x})_{i+1},...,(\mathbf{x})_m$ ]と置く. エルゴード性を過程
1. 適当に $i$ を選択する. (ランダムとか順番とか)
2. $(\mathbf{X}^{(t+1)})_i | \mathbf{x}^{(t)}\ \sim\ \ g((\mathbf{x})_i|(\mathbf{x}^{(t)})_{-i}), \ \ \ (\mathbf{x}^{(t+1)})_{-i} = (\mathbf{x}^{(t)})_{-i}$ として定める. --> 1.に戻る
(条件付き分布 $g((\mathbf{x})_i|(\mathbf{x})_{-i}) := \frac{g(\mathbf{x})}{\int g(\mathbf{x})d(\mathbf{x})_i}$ である. )
===================================

※ これはMH法の特殊なケースに対応.
$h_i([(\mathbf{x})_1,...,(\mathbf{x})_{i-1},(\mathbf{v})_i,(\mathbf{x})_{i+1},...,(\mathbf{x})_m]|\mathbf{x}) = g((\mathbf{v})_i | (\mathbf{x})_{-i}) = \frac{g(\mathbf{v})}{\int g(\mathbf{x})d(\mathbf{x})_i}$
これに対して, $\frac{h_i(\mathbf{x|v})g(\mathbf{v})}{h_i(\mathbf{v|x})g(\mathbf{x})} = \frac{g((\mathbf{x})_i|(\mathbf{v})_{-i})g(\mathbf{v})}{g((\mathbf{v})_i|(\mathbf{x})_{-i})g(\mathbf{x})} = \frac{\frac{g(\mathbf{x})}{\int g(\mathbf{x})d(\mathbf{x})_i}g(\mathbf{v})}{\frac{g(\mathbf{v})}{\int g(\mathbf{x})d(\mathbf{x})_i}g(\mathbf{x})} = 1$ より,
選択確率が1のMH法に対応している事が言える.

※ 実際の応用では1の $i$ の選択は番号順にやる事が多いらしい.

MCMCを応用した画像のノイズ除去

膨張・収縮を使った方が実用的レベルでは速いと思います.

原理

ノイズが乗る前の画像(未知)を $H$ , ノイズが乗った後の画像(既知)を $V$ と置く.
尤度関数 $P(H|V)$ に従って画像を生成.

同時確率は $P(H, V) = \frac{1}{Z}\exp \left( \alpha \sum_{i} \sum_j v_{ij} h_{ij} + \beta \sum_{i',j' \in N(i,j)} h_{ij}h_{i'j'} \right)$ として定義する.
ただし, $Z$ は正規化項, $N(i,j)$ は $(i,j)$ の近傍. $\alpha, \beta$ はハイパーパラメータ.

画像データの配列の要素は $\pm 1$ が入るので, この様に定義すると, 周囲の点と同じ値を持つほど確率が上がり,
持たないほど下がる様に同時確率を構成できる.

結果と実装

ギブスサンプリングを用いて尤度関数を近似する.
結果は以下の通り.

↑ノイズを載せた画像(よく見るアレ)

↑ノイズ除去後の画像

実装は以下.

import numpy as np
from tqdm import tqdm
import matplotlib.pyplot as plt

def load_image(filename):
	tmp = plt.imread(filename)
	img_gray = 0.299*tmp[:,:,2] + 0.587*tmp[:,:,1] + 0.114*tmp[:,:,0]  # grayscale
	img = np.where(img_gray>0.5, 1, -1)
	fixed_img = np.zeros([img.shape[0]+2, img.shape[1]+2])
	fixed_img[1:-1, 1:-1] = img
	return fixed_img

# V:ノイズ有りの画像, H:元の画像
# 点(i,j)でのP(H|V) := 1/Z exp(v[i,j]*h[i,j] + \sum_{k,l \in N(i,j)}h[i,j]h[k,l] )  ただしN(i,j)は(i,j)の近傍
# H[i,j]をh_ijと書く事にすると, P(h_ij=1 |V,H_N(i,j))が求まる.

def denoise(noised_pict, total_step, burnin_step, alpha, beta):
	ans = np.random.choice([1, -1], size=noised_pict.shape)
	x_max, y_max = ans.shape[0], ans.shape[1]
	for _ in tqdm(range(burnin_step+total_step)):
		for i in range(1, x_max-1):
			for j in range(1, y_max-1):
				tmp = beta*np.sum([ans[i-1+l][j-1+m] for l in range(3) for m in range(3)]) + alpha*noised_pict[i][j]
				prob_h = np.exp(tmp)/(np.exp(tmp)+np.exp(-tmp))  # H[i,j]=1となる条件付き確率
				ans[i][j] = 2*(np.random.random_sample() < prob_h) -1.0
	return ans

if __name__ == "__main__":
	file_name = "./img_noised.png"
	tmp = load_image(file_name)
	plt.imshow(tmp, cmap="gray")
	plt.show()
	ans = denoise(tmp, 1000, 100, 1.0, 0.5)
	plt.imshow(ans, cmap="gray")
	plt.show()

雑な感想と近況報告

詳細釣り合いと定常分布の関係を書いた記事があまり多くなかった気がするので書きました.

2021年はずっとロボコンに追い回されていました. 夏休みを返して...(恨み)

Emileの備忘録

色々な事をだらだらと