PRML 12章主成分分析を EM アルゴリズムで解いてみる

PCA を試す、PPCA を試すとくると、次は確率的主成分分析を EM アルゴリズムで解いてみよう。今回も R で実装。

さすがにそろそろコードが長くなってきたので、全体は github にて。

http://github.com/shuyo/iir/blob/master/pca/ema.r

E-step と M-step のところくらいは載せておこうか。
しかし R にはどうして行列のトレースをとる関数がないのだろう。sum(diag(〜)) なんて、微妙にカッコ悪いんだけど。そう言えばノルムを取る関数もなかったな……。

		### E-step:

		# E[z_n] = M^-1 W^T (x_n - x^bar) (PRML 12.54)
		Ez <- t(M_inv %*% t(W) %*% xn_minus_x_bar);

		# E[z_n z_n^T] = sigma^2 M^-1 + E[z_n]E[z_n]^T (PRML 12.55)
		Ezz <- list();
		sum_Ezz <- matrix(numeric(M*M), M);
		for(n in 1:N) {
			ezz <- sigma2 * M_inv + Ez[n,] %*% t(Ez[n,]);
			Ezz[[n]] <- ezz;
			sum_Ezz <- sum_Ezz + ezz;
		}

		### M-step:

		# W_new = {sum (x_n - x^bar)E[z_n]^T}{sum E[z_n z_n^T]}^-1 (PRML 12.56)
		W <- xn_minus_x_bar %*% Ez %*% solve(sum_Ezz);

		# sigma_new^2 = 1/ND sum{ |x_n-x^bar|^2 - 2E[z_n]^T W^T (x_n-x^bar) + Tr(E[z_n z_n^T] W^T W) } (PRML 12.57)
		sigma2 <- sum(xn_minus_x_bar^2) - 2 * sum(diag(Ez %*% t(W) %*% xn_minus_x_bar));
		for(n in 1:N) {
			sigma2 <- sigma2 + sum(diag(Ezz[[n]] %*% t(W) %*% W));
		}
		sigma2 <- sigma2 / N / D;

E[z_n z_n^T] が行列の列なので、ちゃんとループ回して list() に突っ込むしかなくて面倒。
それ以外は、数式を素直にコードに落とすだけ。
E[z_n] のような「ベクトルの列」は行列で持つようにしてしまえば、結構１行で計算が出来てしまう。
例えば σ^2 の更新式の中にある $\sum_{n=1}^N \mathbb{E}[\bf{z}_n]^T \bf{W}^T (\bf{x}_n - \bar{\bf{x}})$ も n で回してしまいそうになるが、Ez を各行が E[z_n] であるＮ×Ｍ行列とし、データ点の $\bf{x}_n - \bar{\bf{x}}$ は xn_minus_x_bar <- t(oilflow) - mu; でＤ×Ｎ行列として扱えば、sum(diag(Ez %*% t(W) %*% xn_minus_x_bar)) で計算できる。
ポイントは sum(diag(〜)) で対角成分のみの和を取って後は捨てているところ。つまり結構無駄な計算をさせてしまっているわけだが*1、単にコードが短いだけではなく、これでも自分でループ回すよりずっと速い。

収束の様子を GIF アニメでどうぞ。

収束すると、回転反転を除いて普通の PCA とほぼ同じ結果になる。
回転反転の具合は初期値依存。W と σ^2 の初期値は正規乱数〜 N(0,1) で埋めている。

EM のいいところは、データ空間が高次元になったときに計算量がお得な場合がある、とのことだが、oil flow データは１２次元しかない。
EM のいいところ・その２は、ランダム欠損値を扱うことが出来ること。PRML 図 12.11 にそれを表しているという図が載っている。これはぜひやってみたいと思い、ちょっと考えた。

PRML 9.3章(p157) には「データ集合中の欠損値を非観測変数と見なして EM アルゴリズムを用いる」とある。
例えば x_ni (n 番目のデータ点の i 番目のフィールド) が欠損している場合、この x_ni を z_n 達と一緒に潜在変数として扱って、EM アルゴリズムの更新式を導出しなくちゃならない。
そのためには事後確率 p(x_ni | X) が必要になるが(本当は条件から x_ni を除かなければならない)、それは p(x_ni | X) = Σ p(x_ni|z_n) p(z_n|X) で、x_ni と z_n は線形の関係があるので計算できて、正規分布になる。
あとは同時分布の対数期待値の式を、x_ni を潜在変数として展開してあげると、E-step では E[x_ni], E[x_ni^2], E[x_ni z_n] を求めれば良さそう。
と、いうところまではあたりをつけたが、ランダム欠損した x_ni ごとに E[x_ni], E[x_ni^2], E[x_ni z_n] を求めて、M-step もそれらにあわせて計算式を実装して……って涙出そうになるほどめんどくさい……。

というわけで日和ってしまいました。ごめんね。次はベイズ PCA。