無限の猿

勉強したことを書き残してます

主成分分析と独立成分分析とスパースコーディングの比較[python]

データを教師なしで変換する行列分解手法、主成分分析(PCA: Principle Component Analysis)、独立成分分析(ICA: Independent Component Analysis)、スパースコーディング(SC: Sparse Coding)の比較。


行列分解手法の明確な定義は知らないが、ここではデータを表すベクトル \bf{x}_n\in R^Dの集合を横に並べた行列を \bf{X}\in R^{D\times N}として、基底を表す行列\bf{U}\in R^{D\times K}と係数のを表す行列\bf{A}\in R^{K\times N}の積、
    \bf{X = UA}
に変換する手法とする。これはすなわち、元のデータ \bf{x}\bf{U}の列にあたる基底 \bf{u}_kの線形和、
    \bf{x} = \sum_k^K \bf{u}_k a_k = \bf{Ua}
で表現することを意味する。 \bf{a}はデータ \bf{x}の基底\bf{U}で表される空間での表現に相当することになる。


ここで\bf{U}が決まっていれば、\bf{A}を求めるのは線形の逆問題(\bf{U}が正則な正方行列であればその逆行列\bf{X}にかけてやれば良い)となるが、行列分解問題では双方を同時に求める問題となる。自由度も高くなり、\bf{U}\bf{A}は一意には決まらない。したがって、\bf{U}\bf{A}に対して他に何らかの拘束とか条件とかを与えてやる必用がある。この条件の違いが行列分解手法の結果の違いにつがなる。


今回興味あるのは基底\bf{U}の方なので、これの求め方をまとめてみる。

主成分分析(PCA: Principle Component Analysis)

    \bf{U} = \rm{argmin}_{\bf{U}}||\bf{X-UA}||^2,\\ s.t.\bf{U}^T\bf{U}=\bf{I}
基底\bf{U}は、元のデータを再現でき、なおかつその列 \bf{u}_kは互いに直交すべし、という式。これにより、許された基底の数Kで最大限データを再現できるような直交基底を得ることができるはずである。実はこれだけではK>1のとき\bf{U}は一意には決まらない。 \bf{u'}_k=\sum^K_i\beta_i \bf{u}_iとなるような直交基底を新たに用いても、それが同じ空間を張るために誤差項は変化しないためである。
実装においては、データ\bf{X}の共分散行列\bf{X}^T\bf{X}/N固有ベクトル固有値の大きい順にK個、基底として採用することで実現できる。この方法の場合、固有値の大きい基底ほど誤差項の減少に寄与するように求められるので、上記の問題は起こらない。さらにこのとき、係数\bf{a}_nの各要素は無相関になる。機械学習分野で有名なビショップさんの黄色い本パターン認識と機械学習 下 (ベイズ理論による統計的予測)やWeb上の様々なところで解説されていると思います。

独立成分分析(ICA: Independent Component Analysis)

    \bf{U} = \rm{argmin}_{\bf{U}}||\bf{X-UA}||^2 - \lambda H(\bf{A;U})
ただし、 \lambdaは独立性の重みパラメータ、 H(\bf{A;U})\bf{A}エントロピーを表す。\bf{a}エントロピーが高いほど、\bf{a}の要素間の独立性が高まる(理解が正しいか微妙)。具体的に言えば、\bf{a}の確率密度分布p(\bf{a})と、要素の周辺分布の積\Pi_k^K p(a_k)がカルバック・ライブラー距離の意味で近くなる(p(\bf{a}) = \Pi_k^K p(a_k)のとき\bf{a}の要素は独立となる)。コスト関数は誤差項を小さくしつつこのエントロピーを大きくすることで小さくなるので、得られる基底\bf{U}は元のデータの再現性を保ちつつ、変換後の係数が独立になるようなものになるはずである。
しかしながら、独立性の測り方は色々あるらしいので詳しくは詳解 独立成分分析―信号解析の新しい世界あたりで。

スパースコーディング(SC: Sparse Coding)

    \bf{U} = \rm{argmin}_{\bf{U}}||\bf{X-UA}||^2 - \lambda \sum^N_n|\bf{a}_n|_1
独立成分分析とは違い、第二項がL1ノルムになっている。Lasso正則化とも呼ばれるらしい。このL1ノルムの項は、係数\bf{a}がスパース(ほとんどの要素がゼロ)になることを要請する。したがって、データごとに少数の基底のみを使って元のデータをよく表現できるような基底が得られるはずである。ここで「データごとに少数の基底のみ」というのは別のデータを表現するのにはまた違った「少数の基底」を使っても良いということである。ややこしく言えば、任意の少数の基底ではられる部分空間になるべく全てのデータが乗るような基底ということになる。多分、基底は多ければ多いほどデータを少数の基底で表現できる可能性が高くなるからだと思うのですが、スパースコーディングを行う際にはデータの次元Dに対して基底の数Kを多くとる(多くの場合K>D)オーバーコンプリートということがよく行われる。
独立成分分析と同様、実装方法は色々ある。同じくスパースコーディングを行うL0正則化ではOrthogonal Matching Pursuit(OMP)、L1正則化ではLars-Lassoとか。イラストで学ぶ 機械学習 最小二乗法による識別モデル学習を中心に (KS情報科学専門書)では解析解が求まるL2正則化で上界を抑える方法を解説している。画像処理の分野でデノイズや超解像なんかで流行っている気がする。あと、詳しく解説してる日本語の参考書が見つからない。

実装

実装においては、主成分分析は固有値計算するのみで計算が早い一方、独立成分分析とスパースコーディングが目的関数の最適化に繰り返しが必要となり、比較的計算が遅い。自分のにわか知識で実装しても計算速度が残念になりそうなので、pythonのライブラリを利用させていただくことにする。

scikit-learn

行列分解手法(主成分分析、独立成分分析他)を含む機械学習手法のライブラリ。

from numpy import *
from matplotlib.pyplot import *
from sklearn.decomposition import PCA
from sklearn.decomposition import FastICA
from sklearn.decomposition import DictionaryLearning

# 混合ガウス分布に従う乱数の生成
def rand_gauss_mix(mu, sigma, p, N):
	d, K = mu.shape
	p_cumsum = cumsum(array(p) / sum(p))
	X = zeros([d, 0])
	for n in range(N):
		p_ = p_cumsum - random.rand()
		p_[p_<0] = 1
		k = p_.argmin()
		x = random.multivariate_normal(mu[:,k], sigma[:,:,k]).reshape(-1,1)
		X = hstack((X, x))
	return X

# 回転行列
Rot = lambda rad: array([[cos(rad), sin(rad)], [-sin(rad), cos(rad)]])


# データの生成
N = 1000
distribution_type = 'rectangle'

if distribution_type == 'gauss':
	mu = array([0, 0])
	sigma = array([[0.1,0.1],[0.2,0.3]])
	X = random.multivariate_normal(mu, sigma, N).T
elif distribution_type == 'rectangle':
	rad = (1./6)*pi
	ext = [1, 0.4]
	X = 2*random.rand(2,N)-ones([2,N])
	X = dot(Rot(rad), X)
	X = dot(diag(ext), X)
elif distribution_type == 'closs':
	mu = zeros([2,2])
	p = [0.5, 0.5]
	sigma0 = diag([0.3,0.003])
	rad = [1./8*pi, 7./8*pi]
	sigma = zeros([2, 2, 2])
	for i in range(2):
		sigma[:,:,i] = dot( Rot(rad[i]), dot(sigma0, Rot(rad[i]).T) )
	X = rand_gauss_mix(mu, sigma, p, N)

# 主成分分析
decomposer = PCA()
decomposer.fit(X.T)
Upca = decomposer.components_.T
Apca = decomposer.transform(X.T).T

# 独立成分分析
decomposer = FastICA()
decomposer.fit(X.T)
Uica = decomposer.mixing_ 
Aica = decomposer.transform(X.T).T

# スパースコーディング
decomposer = DictionaryLearning()
decomposer.fit(X.T)
Usc = decomposer.components_ .T
Asc = decomposer.transform(X.T).T


axis('equal')
plot(X[0],X[1],'xc')
Upca = Upca / sqrt((Upca**2).sum(axis=0))
Uica = Uica / sqrt((Uica**2).sum(axis=0))
Usc = Usc / sqrt((Usc**2).sum(axis=0))
for i in range(2):
    p_pca = plot([0, Upca[0,i]], [0, Upca[1,i]], '-r')
    p_ica = plot([0, Uica[0,i]], [0, Uica[1,i]], '-b')
    p_sc = plot([0, Usc[0,i]], [0, Usc[1,i]], '-g')
legend(('data', 'PCA', 'ICA', 'SC'))
legend(loc="best", prop=dict(size=12))
show()

subplot(1,3,1)
plot(Apca[0], Apca[1], 'xc')
title('PCA')
subplot(1,3,2)
plot(Aica[0], Aica[1], 'xc')
title('ICA')
subplot(1,3,3)
plot(Asc[0], Asc[1], 'xc')
title('SC')
show()

実験

十字の分布の時

細いガウス分布を交差させて作った十字型の分布(distribution_type = 'closs')から2次元データをサンプリングして、主成分分析、独立成分分析、スパースコーディングをそれぞれ行った結果。
f:id:shiten4:20140110122514p:plainf:id:shiten4:20140110122524p:plain
上図赤線、青線、緑線はそれぞれから得られた2つの基底ベクトルの方向を表す。青線と緑線は左上、左下で重なっているので注意。下図はデータをそれぞれの基底を用いて変換した結果。すなわちデータを基底の線形結合で表した場合の係数\bf{a}_nのプロット。
まず主成分分析(赤)は固有値計算の結果に従って、第一基底を1番分散の大きい方向、第二基底をそれに直交する方向にとっている。分布はガウス型にはなっていないが、この×型の分布は横長なので、ほぼ真横に第一基底が取られている。下の図でも変換後にはデータは無相関になっていることが確認できる。
独立成分分析(青)は、×の細いラインに沿った方向に基底をとっている。変換されたデータも、縦横の軸上の周りにのみ分布するようになっている。×字に重なっているそれぞれの分布を独立した発生源とみなして分離していると考えられる。
スパースコーディング(緑)も、独立成分分析と同様に×の細いラインに沿った方向に基底をとっている。これはそれぞれのデータをできるだけ少ない基底(この場合は1つだけ)で表現しようとした結果だと思う。変換後の結果は独立成分分析以上にデータが縦横軸上に集中した、スパースなものになっていることが分かる。

平行四辺形型の分布の時

f:id:shiten4:20140110122530p:plainf:id:shiten4:20140110122536p:plain
主成分分析はデータを無相関なものに変換しているのみであるが、独立成分分析は菱型の辺と平行に基底を取り、データを辺が軸と平行な四角形の分布に変換している。様々な信号が重なりあうほど中心極限定理により分布がガウス分布に近づく傾向にあるので、データの独立性の測り方の一つとして、その分布なるべくガウス分布と異なる形であるかを用いることができる。独立成分分析の結果は、変換されたデータの各要素が一様分布になっているという意味で、きちんと独立な成分に分離されているようです。一方スパースコーディングはこの実験では独立成分分析の結果とは全く異なる基底をとっている。この菱型の分布はそれぞれのデータを一つの基底のみで表現するのは厳しい。頑張ってなるべくデータを再現できるように基底をとっているのであろうが、この結果では基底の軸上にしか再構成できないので、適切な行列分解もとい基底の取得ができていないと思われる。