カーネル主成分分析をゆるふわに理解する(2/2)
「カーネル多変量解析」をゆるふわに読んでいて、第3章に入ってきました。第3章の最初のテーマとして、カーネル主成分分析を取り上げています。
前回は、カーネル主成分分析の概要を紹介した上で、その手順をステップバイステップで導出しました。しかしその過程で、一つの誤魔化しをしました。この誤魔化しがない形に修正することが、今回の記事の目的です。
また、その修正が終わった上で、カーネル主成分分析を実装してみます。
前回のふりかえり
カーネル主成分分析は、そのままでは線形分離しにくいk次元データに対して、超高次元のm次元に射影することで良い感じの主成分軸を取り出し、主成分をd個選ぶことでk>dへ圧縮する前処理でした。その中で一度m次元を通るはずなのですが、カーネルトリックでその計算をスキップしてしまうことが面白いところでした。
前回のカーネル主成分分析手順の導出では、あるk次元のサンプルデータ$${x^{(i)}}$$のn個の集まりに対して、$${\phi(x)}$$でm次元に射影後に主成分に射影する関数$${g(x)}$$を求め、$${g(x)}$$から重要な主成分d個だけを抽出することでカーネル主成分分析$${f(x)}$$を得ました。
その途中で、不完全な$${g(x)=\sum^n_ia_ik(x^{(i)}, x)}$$が得られました。あとは$${g(x)}$$の中にある$${a=(a_1,a_2,…,a_n)}$$を求めれば、$${g(x)}$$を同定でき、そこから簡単に$${f(x)}$$を求められます。
この$${a=(a_1,a_2,…,a_n)}$$を求める方法として、$${g(x)}$$を用いてサンプルデータを主成分に射影した場合に、その結果の分散が最大になるようにすべきだと考えました。
$${g(x)}$$を用いてサンプルデータを主成分に射影した結果は（カーネル行列$${K}$$を用いて）$${Ka}$$とかけます（詳細は前回参照）。この分散を求めたい。その際に、分散は$${1/n|Ka|^2}$$と計算して後続の説明を続けました。しかし、$${x}$$が標準化されていて平均値が0だとしても、$${Ka}$$の平均値が0であるとは限りません。
もし、$${Ka}$$の平均値が0でなかった場合、果たして$${f(x)}$$をどうやって求めればいいのでしょうか？
カーネル行列の中心化
では改めて、$${Ka}$$の分散を求めてみましょう。
まず、普通の分散の式を使いやすい形に変換します。
$${1/n\sum^n_i(x_i-\bar{x})^2}$$
= $${1/n\sum^n_i(x_i^2-2x_i\bar{x}+\bar{x}^2)}$$
= $${1/n\sum^n_i(x_i^2)-1/n\sum^n_i(2x_i\bar{x})+1/n\sum^n_i(\bar{x}^2)}$$
= $${1/n\sum^n_i(x_i^2)-2\bar{x}/n\sum^n_i(x_i)+\bar{x}^2}$$
= $${1/n\sum^n_i(x_i^2)-\bar{x}^2}$$
つまり、分散は二乗和平均から平均の二乗を引けば出ます。
ということで、$${Ka}$$の分散も、二乗和平均から平均の二乗を引いて求めます。
二乗和平均は$${1/n|Ka|^2=1/n(a^TK^2a)}$$でした。
平均の二乗を計算してみます。わかりやすさのため$${Ka=b}$$とおいて$${b=(b_1,b_2,…,b_n)^T}$$と書くと、平均の二乗は$${(b_1+b_2+…+b_n)^2/n^2}$$です。行列で書くと
$${1/n^2\begin{pmatrix}b_1&b_2&…&b_n\end{pmatrix}\begin{pmatrix}1\\1\\…\\1\end{pmatrix}\begin{pmatrix}1&1&…&1\end{pmatrix}\begin{pmatrix}b_1\\b_2\\…\\b_n\end{pmatrix}}$$
ここで、1だけでできたnxn行列を$${1_n}$$と書くとすると、上の式はこうなります($${Ka=b=(b_1,b_2,…,b_n)^T}$$も利用しながら）
$${1/n^2(Ka)^T1_nKa = 1/n^2a^TK^T1_nKa=1/n^2a^TK1_nKa}$$
ということで$${Ka}$$の分散は
$${1/n(a^TK^2a)-1/n^2a^TK1_nKa=1/na^TK(Ka-1/n1_nKa)=1/na^TK(I_n-1/n1_n)Ka}$$
長々と書きましたが、まとめると、$${Ka}$$の平均値が0と仮定したときの分散は$${1/n(a^TK^2a)}$$だったけど、平均値が0でないなら$${1/na^TK(I_n-1/n1_n)Ka}$$というだけの話です。
この分散が最大になる$${a}$$を求めるためにラグランジュの未定乗数法を使うと、前回とちょっとだけ違う感じになります。
$${(I_n-1/n1_n)Ka=\lambda a}$$
つまり$${(I_n-1/n1_n)K}$$の固有ベクトルを$${a}$$としたとき、主成分上のサンプルデータ写像の分散が最大になります。これで$${g(x)}$$がもとまるので、後は固有値が大きい順にd本$${a}$$をとれば、カーネル主成分分析$${f(x)}$$が求まります。
これで完全に手順が導出されました。
更新された手順を実装してみる
ということで更新された手順はこうなります。3)の$${K}$$が$${(I_n-1/n1_n)K}$$になっただけです。
1) k次元の$${x}$$とその結果の$${y}$$について、n個のサンプル$${x^{(1)},x^{(2)},…,x^{(n)}}$$を用意する
2) カーネル関数を決め（例えばガウスカーネル）、n個のサンプルから$${K}$$を計算する
3) $${(I_n-1/n1_n)K}$$の固有ベクトル$${a}$$を計算する
4) 固有値が大きい順にd個の固有ベクトルを取得して並べる
5) 求めたいデータに対して各サンプルデータと$${k(x^{(i)}, x)}$$を計算し、4)の行列をかければd次元の主成分上の点を示す値が得られる
導出長かった……。けどなんか理解できた感があります。ということで実装してみましょう。
まずクラスを定義します。カーネル関数はデフォルトでガウスカーネルを使うようにしますが、後でも変えられます。
import numpy as np
def gaus(Xi, Xj, gamma):
dif = Xi - Xj
return np.exp(-gamma * dif.T.dot(dif))
class KCPA():
def __init__(self):
self.k_func = lambda Xi, Xj: gaus(Xi, Xj, gamma=15)$${f(x)}$$をfit_transformという関数として定義します。サンプルデータXと射影したい次元dを受け取ります。内部では、$${K}$$を計算した後に中心化、固有ベクトルをd個取得して、Xの主成分上の値を返します。
def fit_transform(self, X, d):
K = self._K(X)
K = self._centralize(K)
alphas = self._alphas(K, d)
return K.dot(alphas)Kを求めるのは簡単。n x n回カーネル関数を計算して、reshapeしてあげるだけです（なんかもっと上手く書けそう）。
def _K(self, X):
row = X.shape[0]
ks = []
for i in range(row):
for j in range(row):
ks.append(self.k_func(X[i], X[j]))
ks = np.array(ks)
return ks.reshape(row, -1)中心化（$${(I_n-1/n1_n)K}$$）。
def _centralize(self, K):
dim = K.shape[0]
j = np.identity(dim) - 1/dim * np.ones((dim, dim))
return j.dot(K)最後に$${a}$$を求める。
def _alphas(self, K, d):
eigen_vals, eigen_vecs = np.linalg.eig(K)
indices = np.argsort(-eigen_vals)#降順にするために-をつけた
return eigen_vecs[:, indices[0:d]]上の5つを並べればクラスができます。一応まとめると
import numpy as np
def gaus(Xi, Xj, gamma):
dif = Xi - Xj
return np.exp(-gamma * dif.T.dot(dif))
class KCPA():
def __init__(self):
self.k_func = lambda Xi, Xj: gaus(Xi, Xj, gamma=15)
def fit_transform(self, X, d):
K = self._K(X)
K = self._centralize(K)
alphas = self._alphas(K, d)
return K.dot(alphas)
def _K(self, X):
row = X.shape[0]
ks = []
for i in range(row):
for j in range(row):
ks.append(self.k_func(X[i], X[j]))
ks = np.array(ks)
return ks.reshape(row, -1)
def _centralize(self, K):
dim = K.shape[0]
j = np.identity(dim) - 1/dim * np.ones((dim, dim))
return j.dot(K)
def _alphas(self, K, d):
eigen_vals, eigen_vecs = np.linalg.eig(K)
indices = np.argsort(-eigen_vals)#降順にするために-をつけた
return eigen_vecs[:, indices[0:d]]じゃあ、前回出たmake_moonsをカーネル主成分分析し、第一主成分と第二主成分の値を二次元上に描画してみましょう。
第一主成分（CP1）で綺麗に分離できています。
なお、この図をみて「カーネル主成分分析がクラスラベルを綺麗にわけた」と勘違いしやすいです。カーネル主成分分析はクラスラベル$${y}$$を受け取っていないので、クラスラベルは知りません。ただ、分散が大きくなるような主成分を選んだら右の図のようにデータが散りばり、クラスラベルを振るとちょうど綺麗に分かれているのです（元図で青丸同士、赤バツ同士が近い結果です）。
カーネル関数による結果の違い
さて実装も終わって結果もみて、めでたしめでたし……とは行きません。実装の時にサラッと書いた$${\gamma=15}$$のガウスカーネルってどこからきたんですかね？
まず、$${\gamma}$$について。$${\gamma}$$の値を$${1, 2, 4, 8, 16, 32, 64, 128, 256}$$に振ってみるとこうなります。$${\gamma}$$が小さいと上手く分けれていないです。
$${\gamma}$$のチューニング必要性がもっとわかりやすい例としてこのデータを使ってみます。
ガウスカーネルでカーネル主成分分析するとこうなります。$${\gamma}$$が2から16くらいじゃないとうまく機能しなそうです。
なぜこうなるかというと、主成分の取得や分散の計算を、射影された高次元でおこなっているからで、つまりカーネル関数（もしくは$${\phi(x)}$$）によって結果がぶれるということになります。
だから、「データにあわせてカーネル関数を選ぶ」必要があるそうです。$${/gamma}$$を変えるだけでなく、ガウスカーネル以外のカーネル関数まで考えて。その具体例がこの後書いてある内容だそうです。
ということで、今回でカーネル主成分分析の概要ゆるふわ理解は終わり。次はデータにあわせてカーネル関数を選ぶ例を見ていきますー。
ではでは。
この記事が気に入ったらサポートをしてみませんか？