PCA、NIPALS
主成分分析は教師なし機械学習の一つで、データをより小さい次元に圧縮（要約）することで観察をより容易にすることを目的としています。主成分とはデータのより意味のあると思われる方向のことを指します。機械学習なので、それをある指標を用いて算出しますが、ここでは意味のある方向＝より分散の大きい方向のことを指します。
X-Y次元平面上に散らばる点を考えてみます。例えば、Y軸方向への分散が非常に小さければ、そのデータを語るうえで重要となるのはX軸のみでよいでしょう。Y=X軸上に点が集中していれば、今度はどちらの軸も重要に重視しなければなりません。
注意したいのが、PCAを行った際に主成分軸に意味をつけるのは人間であるということです。PCAは指標に従い変換しますが、第１主成分軸方向はどのような成分が含まれているのかを確認しつつ、意味付けを行います。軸の方向もあまり意味がなく、第１主成分軸に例えば年収の成分が多く含まれているからと言って、正の方向にそのまま年収の高いサンプルがより存在しているかというとそうではないこともあります。
データの観察という利用方法以外でも、データの次元数が多く、計算資源的に厳しい場合にPCAの上からN個分の主成分のみを用いたりもします。例えば画像を深層学習でよく用いられているモデルを用いて特徴量化（ベクトルへの変換）を行うことがあります。その時深層学習モデルからの出力は数百次元など非常に高次元になってしまいます。これを主成分３２個のみ抽出したりということです。
PCA
では、主成分をどのように求めるかについてみてみます。
アルゴリズム
入力 $X$ (NサンプルxD次元行列)を用意します。
- 入力 $X$ の各次元の平均値 $\vec{\mu}$ を計算し、それを $X$ から差し引いたものを $\tilde{X}$ とします。
- $\tilde{X}$ の共分散行列 $\Sigma$ を求める
- $\Sigma$ の固有値と固有ベクトルを求め、固有値の絶対値が大きい順に固有ベクトルを並べた行列 $U=[v_1, ,,, , v_D]$ を求める
- $U$ と $X$ の積 $XU$ を求める
以上で計算することができます。Scikit-learnの結果と比較したのですが、符号が反転している箇所がありました。これは実装による違いだと思いますが、一貫した分析ができるのであれば特に問題はないと思っているため（＋整合性を合わせるための調査をしていないため）修正は行っていません。どうやらScikit-learnでは符号を強制的に変更しているようです。（そもそもSVDを用いて求めているようですが。PCAとSVDの関連についての記事もご覧いただいた方が良いと思います）
NIPALS
上記のアルゴリズムでは、すべての主成分を求めていました。しかし、最初に述べたようにそのような状況はあまりなく、主成分を上から数個でいいような状況がほとんどです。そのような場合には、ほかのアルゴリズムを用いたほうが良いでしょう。Scikit-learnでは内部で分岐が行われています。n_componentsつまりいくつの主成分を残すかによって異なります。
- max(N,D)<=500 または、最尤推定で必要な数を求める場合
- 全主成分を求めてから指定した一部の主成分のみに絞り込みます
- 前者の条件はサンプル数も特徴量数も少ない場合です。コメントにこの条件はSmall problemである書いてあります。
- 1 <= n_components < 0.8 * min(N,D)
- サンプル数、特徴量数のいずれかの８割より少ない場合です。この場合はRandomizedPCAを用いているようです。
- それ以外、例えばn_componentsにDと同じ程度の数値を指定している場合は全主成分を求めてから絞り込みます
上記はsvd_solverをautoにしてした時ですが、svd_solver="arpack"ではscipy内の関数を利用しています。
NIPALS(Nonlinear Iterative Partial Least Squares)も第１主成分から順に求めるアルゴリズムです。こちらのページを参考にしました。各主成分ごとにいくつかのイテレーションを回していきます。PCA同様に入力 $X$ (NサンプルxD次元行列)を用意します。平均を差し引いたものとします。
- $X$の第１列目(主成分ではない) $t_1$ を取得。$t_1$ の要素数はNと一致する
- $p_1 = X^T t_1 / t_1^T t_1$ を計算。$p_1$ の要素数はDと一致する
- $p_1$ を正規化
- $t_1 = X p_1 / p_1^T p_1$ で更新
- 収束するまで２～４を繰り返す。おそらく $t_1$ の前イテレーションからの差分が少なくなったタイミング
- 第２主成分以降を求める場合は、$X = X - t_1 p_1^T$ で入力を更新してから再度１～５（添え字は適宜更新）を繰り返す
import numpy as np
from sklearn.decomposition import PCA
from sklearn.datasets import load_iris
X = load_iris()["data"]
class NIPALS:
def __init__(self, n_components, max_iter=100, tol=1e-6):
self.n_components = n_components
self.max_iter = max_iter
self.tol = tol
def fit_transform(self, X):
self.mean_X = np.mean(X, axis=0)
self.X_shift = X - self.mean_X
out = np.zeros((X.shape[0], self.n_components))
for n in range(self.n_components):
th_old = self.X_shift[:,n]
for itr in range(self.max_iter):
ph = np.dot(self.X_shift.T, th_old) / np.sum(th_old**2)
ph = ph / np.sqrt(np.sum(ph**2))
th = np.dot(self.X_shift, ph) / np.sum(ph**2)
if np.linalg.norm(th_old - th)<=self.tol:
break
th_old = th
out[:,n] = th
th_x_ph = np.matmul(th.reshape(X.shape[0], 1), ph.reshape(1, X.shape[1]))
self.X_shift -= th_x_ph
return out
nipals = NIPALS(n_components=2, tol=1e-6)
print("*"*100)
print("*"*100)
print(" NIPALS")
print(nipals.fit_transform(X)[0:10, :])
print("*"*100)
print("*"*100)
print(" PCA: FULL")
pca = PCA(n_components=2)
print(pca.fit_transform(X)[0:10,:])
実装しては見たのですが、第１主成分のみもとめる場合でも全主成分を求めるものに比べておそかったです。原因としては、実装の仕方云々もあるとは思うのですが、そもそもの入力データが適していなかったのではないかと思っています。二つに手法では計算量がおおよそ以下のようになると思います
- 共分散分析で求める方法
- 共分散分析の計算量：行列の１要素あたりNを、約$D^2$分求めます。$O(N D^2)$
- Jacobi法の計算量：DxDの正方行列なので、イテレーションk回とすると$O(k D^2)$
- NIPALSで求める方法(各主成分ごと)
- ステップ２と４での行列ベクトル積の計算量：$O(ND)$ を数イテレーション繰り返す
上記のJacobi法の計算量があまり問題にならない領域（～５０次元）のためだと思います。もっと高次元で行えば処理時間が逆転すると思います。
その他の手法について
おととしのアドベントカレンダーや、そこで引用されていた論文を参考にしてください。
以上