[ベイズ推論による機械学習入門]numpyだけで線形次元削減を実装してみた。

須山敦志著の『ベイズ推論による機械学習入門』の第5章を読みました。5.1章 線形次元削減をpythonで実装してみました。全コードは私のgitにあります。

この投稿では以下の次元削減アルゴリズムを用いて、入力画像のデータ量を圧縮します。

###記号の使い方

小文字の太字は縦ベクトル、大文字の太字は行列は表します。$p$次元の単位行列を${\bf I}_p$であらわします。
$D$次元ガウス分布を、平均パラメータベクトル$\boldsymbol{\mu}\in \mathbb{R}^D$, 共分散行列${\bf \Sigma}\in \mathbb{R}^{D \times D}$を用いて
$$
{\mathcal N}({\bf x}|{\bf \mu},{\bf \Sigma})=\frac{1}{\sqrt{(2\pi)^D|{\bf \Sigma}|}} {\rm exp} \left\{ -\frac{1}{2}({\bf x}-\boldsymbol{\mu})^{\rm T}{\bf \Sigma^{-1}({\bf x}-\boldsymbol{\mu})} \right\}
$$
で定義します。

#理論

詳細は『ベイズ推論による機械学習入門』の5.1章を読んでください。ここでは、前提と結論のみを書くに留めます。

##前提

観測されたデータ$Y=\{{\bf y}_1,\cdots,{\bf y}_N\}$を低次元の潜在変数$X=\{{\bf x}_1,\cdots,{\bf x}_N\}$で表現します。ここで、${\bf y}_n\in \mathbb{R}^D$,${\bf x}_n\in \mathbb{R}^M$で$M<D$とします。ここで、行列パラメータ${\bf W} \in \mathbb{R}^{M \times D}$およびバイアスパラメータ$\boldsymbol{\mu}\in\mathbb{R}^D$を用いて

\begin{align}
p({\bf y}_n|{\bf x}_n,{\bf W},\boldsymbol{\mu}) = {\mathcal N}({\bf y}_n|{\bf W}^{\rm T}{\bf x}_n+\boldsymbol{\mu},\sigma_y^2{\bf I}_D)
\end{align}

という確率分布を考えます。つまり、データ量$MN$から$MD+DN$となるデータ圧縮を行うわけです。$\sigma_y^2$は分散パラメータです。${\bf W},{\boldsymbol \mu},{\bf X}$に関しても、平均がゼロベクトルになるようなガウス分布を仮定します。

\begin{align}
p({\bf W}) = \prod_{d=1}^D {\mathcal N}({\bf W}_d|{\bf 0},{\bf \Sigma}_{\omega}) \\
p(\boldsymbol{\mu}) = {\mathcal N}(\boldsymbol{\mu}|{\bf 0},{\bf \Sigma}_\boldsymbol{\mu}) \\
p({\bf x}_n)={\mathcal N}({\bf x}_n|{\bf 0},{\bf I}_D)
\end{align}

ここで共分散行列${\bf \Sigma}_\boldsymbol{\mu},{\bf \Sigma}_{\omega}$はハイパーパラメータです。また、${\bf W}_d\in\mathbb{R}^M$は行列${\bf W}$の$d$番目の列ベクトルです。

##結論

観測データ${\bf Y}$が与えられたあとの事後分布は平均場近似を用いて

\begin{align}
p({\bf X},{\bf W},{\boldsymbol\mu}|{\bf Y}) \approx q({\bf X})q({\bf W})q({\boldsymbol \mu})
\end{align}

と書き直します。ここで、頑張って計算すると

\begin{align}
q({\boldsymbol\mu})={\mathcal N}(\boldsymbol{\mu}|\hat{{\bf m}}_\boldsymbol{\mu},\hat{{\bf \Sigma}}_\boldsymbol{\mu}) \\
q({\bf W}) = \prod_{d=1}^D {\mathcal N}({\bf W}_d|\hat{\bf{m}}_{\omega_d},\hat{{\bf \Sigma}}_{\omega}) \\
q({\bf X}) = \prod_{n=1}^N {\mathcal N}({\bf x}_n|\hat{\boldsymbol{\mu}}_{{\bf x}_n},\hat{{\bf \Sigma}}_x)
\end{align}

を得ます。ここで、

\begin{align}
\hat{{\bf m}}_{\boldsymbol \mu} = \sigma_y^2{\bf \Sigma}_{\boldsymbol \mu}\sum_{n=1}^N({\bf y}_n-\hat{{\bf m}}_{\omega}^{\rm T}\hat{{\boldsymbol \mu}}_{{\bf x}_n}) \\
\hat{\Sigma}_{\boldsymbol \mu}^{-1} = N\sigma_y^2{\bf I}_D+{\bf \Sigma}_{\boldsymbol \mu}^{-1} \\ \\

{\hat {\bf m}}_{\omega_d} = \sigma_y^2 {\hat {\bf \Sigma}}_{\omega}\sum_{n=1}^N(y_{n,d}-\hat{m}_{\mu,d}{\hat{\boldsymbol \mu}}_{{\bf x}_n}) \\
{\hat {\bf \Sigma}}_{\omega}^{-1} = \sigma_y^2\sum_{n=1}^N(\hat{{\boldsymbol \mu}}_{{\bf x}_n}\hat{{\boldsymbol \mu}}_{{\bf x}_n}^{\rm T}+\hat{\bf \Sigma}_x)+{\bf \Sigma}_{\omega}^{-1} \\ \\

\hat{{\boldsymbol \mu}}_{{\bf x}_n} = \sigma_y^2 \hat{{\bf \Sigma}}_x {\bf m}_{\omega}({\bf y}_n-\hat{{\bf m}}_{\boldsymbol \mu}) \\
\hat{\Sigma}_x^{-1} = \sigma_y^2\sum_{d=1}^D(\hat{{\bf m}}_{\omega_d}\hat{{\bf m}}_{\omega_d}^{\rm T}+\hat{{\bf \Sigma}}_{\omega}) + {\bf I}_{M}

\end{align}

です。

#実装

上記のアルゴリズムを用いて、入力画像${\bf Y}$の次元削減を行った後に、元画像を復元するプログラムを書きました。復元画像$\bar{\bf Y}$は

\begin{align}
\bar{\bf y}_n = {\hat{\bf m}}_{\omega}^{\rm T} {\boldsymbol \mu}_{{\bf x}_n} + {\hat {\bf m}}_{\boldsymbol \mu}
\end{align}

で求まります。

以下、ハイパーパラメータは著者のgitのコードを参考に決めました。

## Input format: code input_image
from tqdm import tqdm
import numpy as np 
import cv2
import sys

if __name__ == "__main__":

	Y = cv2.imread(sys.argv[1],0)
	Y = cv2.resize(Y,(512,512))

	D = Y.shape[0]
	N = Y.shape[1]
	M = int(64*4)

	MAXITER = 20 #How many times will be repeted

	#Hyper Parameters
	Sigma_omega = np.identity(M)*0.1
	Sigma_mu = np.identity(D)

	sigma_y = 31

	#Initial State
	m_mu = np.random.rand(D).reshape((D,1))
	m_omega = np.zeros((M,D))
	mu_x = np.random.rand(M,N)

	Sigma_mu_hat = np.random.rand(D,D)
	Sigma_oemga_hat = np.random.rand(M,M)
	Sigma_x_hat = np.random.rand(M,M)

	for i in tqdm(range(MAXITER)):

		#mu is reloated

		sum1 = np.zeros((D,1))
		for n in range(N):
			sum1 = sum1 + Y[:,n].reshape((D,1)) - np.dot(m_omega.T,mu_x[:,n].reshape((M,1)))

		m_mu = sigma_y**(-2) * np.dot(Sigma_mu_hat , sum1)
		Sigma_mu_hat = np.linalg.inv(N * sigma_y**(-2) * np.identity(D) + np.linalg.inv(Sigma_mu))

		#W is reloated

		for d in range(D):
			sum2 = np.zeros((M,1))
			sum3 = np.zeros((M,M))
			for n in range(N):

				sum2 = sum2 + (Y[d,n] - m_mu[d])*mu_x[:,n].reshape((M,1))
				sum3 = sum3 + np.dot(mu_x[:,n].reshape((M,1)),mu_x[:,n].reshape((1,M))) + Sigma_x_hat
			
			m_omega[:,d] = sigma_y**(-2) * np.dot(Sigma_oemga_hat , sum2).reshape((1,M))

		Sigma_oemga_hat = np.linalg.inv(sigma_y**(-2) * sum3 * np.linalg.inv(Sigma_omega))

		#X is reloated
	     
		sum4 = np.zeros((M,M))
		for d in range(D):
			sum4 = sum4 + np.dot(m_omega[:,d].reshape((M,1)),m_omega[:,d].reshape((1,M))) + Sigma_oemga_hat
		Sigma_x_hat = np.linalg.inv(sigma_y**(-2)*sum4 + np.identity(M))

		for n in range(N):		
			mu_x[:,n] = ( sigma_y**(-2)*Sigma_x_hat @ m_omega @ (Y[:,n].reshape((D,1)) - m_mu) ).reshape((1,M))
		
	reductioned_Y = m_omega.T @ mu_x + m_mu
	cv2.imwrite("reductioned_img.jpg",reductioned_Y)

#実験結果

入力画像${\bf Y}$に対して$M$を変化させていったときの復元画像の様子は以下の通りになりました。