PRML 12章ベイズ的主成分分析を R で

はてなダイアリーがリニューアルしたらしいので、R で主成分分析を実装してみよう。
PCA を試す、PPCA を試す、EMアルゴリズムでPCAを解く、まで済んだので、次はベイズ的主成分分析。

コード全体は github にて。

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

幸いというか、EM アルゴリズムとの違いは、成分毎の精度を制御するパラメータ alpha が増えたのので、W^new の更新式がちょっと変わったくらい。

		### M-step:

		# W_new = {sum (x_n - x^bar)E[z_n]^T}{sum E[z_n z_n^T] + sigma^2 A}^-1 (PRML 12.63)
		W <- xn_minus_x_bar %*% Ez %*% solve(sum_Ezz + diag(sigma2 * alpha));

		# 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(t(W) %*% xn_minus_x_bar %*% Ez));
		for(n in 1:N) {
			sigma2 <- sigma2 + sum(diag(Ezz[[n]] %*% t(W) %*% W));
		}
		sigma2 <- sigma2 / N / D;

		# alpha_i = D / w_i^T w_i (PRML 12.62)
		alpha <- D / diag(t(W) %*% W);

こんな感じ。

ベイジアン PCA のメリットはなんと言っても ARD(関連度自動決定)。
なんと主部分空間の次元 M を自動的に決定してくれるのだ！

では、M を大きくしていって見よう。
まずは M=2。

PCA / PPCA / PCA with EMA とちゃんと同じ絵が出ている。

では M=3。
３次元の内、軸を２つずつとって３面図を書いてみた。

うーん。縮退している様子はない。αも 18 前後
M=2 とも似ているようで似ていない(理由は後述)。

Ｍとして「データ空間の次元マイナス１」をとれば、完全なサンプル共分散を再現する(PRML p294 の訳注参照)ということであり、そこではきっと ARD も確認できるだろう、ということでいっきに M=11 へ。
同じように１１次元から２つずつとって散布図を描いたけど、５５面図になってしまって多すぎるので、とりあえず一部だけ。

うーん、つぶれてないなあ……
αは次のような値に収束していた。

M=11, I=100000, alpha=( 364.59, 12.14, 17.72, 91.93, 4066.70, 195.10, 29.29, 365.44, 737.34, 68.07, 1175.07)

精度パラメータαが無限大になれば、その要素はつぶれたことになるのだが、最大でも 4066 程度なので、全然つぶれていない。

また M=11 くらいになるとめちゃめちゃ収束が遅くて、数万回くらいイテレーションしないと行けない。
上の例では 10万回回している(尤度の計算がめんどくさいのでやってない*1……モンテカルロ積分するのもあり？)。

まあこれはきっと ARD が効かない、というより、ベイズっつっても線形に過ぎないわけで、そんな線形変換で次元削減できるなんて甘い話はそうそうない、ということなんじゃあないかな。
Oil Flow のようなデータは、当然きっとデータ空間の非線形な部分多様体に乗っていて、線形な切り口ではどこから見てもつぶれてないんだろう。
というわけでここは最終回、非線形なカーネル PCA に期待。

【追記】
M=15(>D=12) でやってみたけど、やっぱり発散しなかった。さすがにこの条件では発散して欲しいよね。うーんまた実装どっか間違ってんのかな……。何度か見直したんだけど。
【／追記】

【さらに追記】
a_i の初期値を大きく(1000〜10000)、M=20 とかにすると縮退するようになった！
が、実行する度に毎回つぶれず残る要素数が異なる。全然つぶれないこともあれば、２個、４個、１０個、１５個などバラバラ。初期値 10000 だと全部つぶれる（もちろんエラーで止まるｗ）こともめずらしくなく。
どうすればもう少し安定してくれるか、模索中。
【／さらに追記】

ダメだったか〜で終わるのも何なので。
さっきの M=3 に戻って、どうして M=2 と似ていないのかというと、回転の自由度があがっていて、元の軸と平行でなくなっているから。
それを見るために、３次元の散布図を描いて、さらにそれを回してみた。700KB over のアニメGIF。