主成分分析(PCA)

主成分分析(Principal Component Analysis, PCA)は高次元データの分散を最大化する直交方向(主成分)を求めることにより、データの低次元表現を構成する次元削減・探索的データ解析の手法である。行列の固有値分解・特異値分解(SVD)を中核とし、最尤推定・正則化・情報量基準・カーネル法と深く結びつく。

主成分分析(PCA)の概念を最初に提唱したのは、現代統計学の父の一人であるカール・ピアソン(Karl Pearson; 1857年3月27日 – 1936年4月27日)。彼は「散布図上のデータ点に対して、それらとの距離の自乗和が最小になるような直線や平面を引く」という幾何学的なアプローチからこの手法を考案。これは現在の「直交最小二乗法」に近い考え方と言える。

その後、約30年の時を経て、経済学者・統計学者のハロルド・ホテリング(Harold Hotelling;1895年9月29日 – 1973年12月26日)が、心理学や教育学における「複数のテスト結果をどうまとめるか」という課題を解決するために、この手法を洗練させた。彼は「主成分(Principal Components)」という名称を正式に使い、行列の固有値分解を用いた現代的な計算手法を確立した。

設定

$n$ 個の観測ベクトル $\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n \in \mathbb{R}^p$ が与えられているとする。データ行列を $X \in \mathbb{R}^{n \times p}$(各行が $\boldsymbol{x}_i^\top$)、標本平均を

\[\bar{\boldsymbol{x}} = \frac{1}{n}\sum_{i=1}^n \boldsymbol{x}_i\]

とする。以降では中心化済みデータ$\tilde{X} = X - \boldsymbol{1}_n \bar{\boldsymbol{x}}^\top \in \mathbb{R}^{n \times p}$を用いる($\boldsymbol{1}_n = (1,\ldots,1)^\top$)。標本共分散行列は

\[S = \frac{1}{n}\tilde{X}^\top \tilde{X}= \frac{1}{n}\sum_{i=1}^n (\boldsymbol{x}_i - \bar{\boldsymbol{x}})(\boldsymbol{x}_i - \bar{\boldsymbol{x}})^\top\in \mathbb{S}^p_+\]

であり、$S$ は対称半正定値行列($\mathbb{S}^p_+$)である。不偏推定量として $S_u = \frac{n}{n-1}S$ を用いる場合もあるが、PCA の固有ベクトルはいずれも一致する。

PCA の定式化:分散最大化

第一主成分

単位ベクトル $\boldsymbol{w} \in \mathbb{R}^p$($\|\boldsymbol{w}\|_2 = 1$)への射影$\boldsymbol{w}^\top \boldsymbol{x}_i$ の標本分散を最大化する方向を第一主成分方向(First Principal Component Direction)と呼ぶ:

\[\boldsymbol{w}_1= \arg\max_{\|\boldsymbol{w}\|_2 = 1} \boldsymbol{w}^\top S \boldsymbol{w}\]

これはラグランジュ乗数法により$S\boldsymbol{w} = \lambda \boldsymbol{w}$ の固有値問題に帰着する。$\boldsymbol{w}^\top S \boldsymbol{w} = \lambda \boldsymbol{w}^\top \boldsymbol{w} = \lambda$ であるから、分散を最大化する方向は $S$ の最大固有値 $\lambda_1$ に対応する固有ベクトル $\boldsymbol{w}_1$ である。このとき第一主成分スコアの標本分散は $\lambda_1$ に等しい。

第 $k$ 主成分の逐次定義

第 $k$ 主成分方向 $\boldsymbol{w}_k$($k = 2, \ldots, p$)は、先の $k-1$ 個の主成分方向と直交するという制約のもとで分散を最大化する:

\[\boldsymbol{w}_k= \arg\max_{\|\boldsymbol{w}\|_2 = 1,\; \boldsymbol{w} \perp \boldsymbol{w}_j\, (j < k)}\boldsymbol{w}^\top S \boldsymbol{w}\]

これは $S$ の固有値分解 $S = W \Lambda W^\top$($\Lambda = \mathrm{diag}(\lambda_1, \ldots, \lambda_p)$、$\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_p \geq 0$、$W = [\boldsymbol{w}_1, \ldots, \boldsymbol{w}_p]$:直交行列)の固有ベクトルを順に取ることと等価である。主成分方向 $\boldsymbol{w}_k$ を主成分ローディング(Loading)、対応する固有値 $\lambda_k$ を主成分分散と呼ぶ。

同時最適化としての定式化

逐次定義と等価な定式化として、$q \leq p$ 次元への射影行列 $W_q = [\boldsymbol{w}_1, \ldots, \boldsymbol{w}_q] \in \mathbb{R}^{p \times q}$を同時に求める問題

\[W_q= \arg\max_{W \in \mathbb{R}^{p \times q},\; W^\top W = I_q}\mathrm{tr}(W^\top S W)\]

が考えられる。目的関数 $\mathrm{tr}(W^\top S W) = \sum_{k=1}^q \boldsymbol{w}_k^\top S \boldsymbol{w}_k$は $q$ 個の主成分方向への射影分散の総和であり、解は $S$ の上位 $q$ 個の固有ベクトルからなる行列 $W_q$ となる。

再構成誤差最小化との等価性

PCA の分散最大化定式化は、データの再構成誤差最小化と等価である。アフィン部分空間への射影$\hat{\boldsymbol{x}}_i = \bar{\boldsymbol{x}} + W_q W_q^\top (\boldsymbol{x}_i - \bar{\boldsymbol{x}})$による再構成誤差の総和を最小化する問題

\[\min_{W_q \in \mathbb{R}^{p \times q},\; W_q^\top W_q = I_q}\frac{1}{n}\sum_{i=1}^n \|\boldsymbol{x}_i - \hat{\boldsymbol{x}}_i\|^2= \min_{W_q} \frac{1}{n}\|\tilde{X} - \tilde{X}W_q W_q^\top\|_F^2\]

の解もまた $S$ の上位 $q$ 個の固有ベクトル $W_q$ であることが示せる。分散最大化と再構成誤差最小化が同一の解を持つことは、フロベニウスノルムの直交不変性と$\|\tilde{X}\|_F^2 = \mathrm{tr}(\tilde{X}^\top \tilde{X}) = n\,\mathrm{tr}(S)$(定数)から直ちに導かれる:

\[\frac{1}{n}\|\tilde{X} - \tilde{X}W_q W_q^\top\|_F^2= \mathrm{tr}(S) - \mathrm{tr}(W_q^\top S W_q)\]

すなわち再構成誤差の最小化は $\mathrm{tr}(W_q^\top S W_q)$(総分散)の最大化と等価である。

特異値分解(SVD)による実装

中心化データ行列 $\tilde{X} \in \mathbb{R}^{n \times p}$ の(薄い)SVD を

\[\tilde{X} = U D V^\top\]

とする($U \in \mathbb{R}^{n \times r}$:左特異ベクトル、$D = \mathrm{diag}(d_1, \ldots, d_r)$:特異値($d_1 \geq \cdots \geq d_r \geq 0$)、$V \in \mathbb{R}^{p \times r}$:右特異ベクトル、$r = \mathrm{rank}(\tilde{X})$)。標本共分散行列との関係は

\[S = \frac{1}{n}\tilde{X}^\top \tilde{X}= \frac{1}{n}V D^2 V^\top\]

であり、$S$ の固有値は $\lambda_k = d_k^2/n$、固有ベクトル(主成分ローディング)は $\boldsymbol{w}_k = \boldsymbol{v}_k$($V$ の第 $k$ 列)となる。主成分スコア行列

\[Z = \tilde{X}W_q = U_q D_q\in \mathbb{R}^{n \times q}\]

は $i$ 行 $k$ 列が観測 $\boldsymbol{x}_i$ の第 $k$ 主成分スコア $z_{ik} = \boldsymbol{w}_k^\top(\boldsymbol{x}_i - \bar{\boldsymbol{x}})$ を与える。

数値計算上の実装は以下の二通りが標準的である。

主成分の性質

無相関性

$k \neq \ell$ に対して第 $k$ と第 $\ell$ 主成分スコアの標本共分散は

\[\frac{1}{n}\sum_{i=1}^n z_{ik} z_{i\ell}= \boldsymbol{w}_k^\top S \boldsymbol{w}_\ell= \lambda_\ell \boldsymbol{w}_k^\top \boldsymbol{w}_\ell= 0\]

となり、主成分スコアは互いに無相関である($S$ の固有ベクトルが直交することから)。

分散の分配

全分散(全変動)は固有値の和として表される:

\[\mathrm{tr}(S) = \sum_{k=1}^p \lambda_k= \sum_{k=1}^p \mathrm{Var}(z_k)\]

第 $k$ 主成分の寄与率(Proportion of Variance Explained, PVE)は

\[\mathrm{PVE}_k = \frac{\lambda_k}{\sum_{j=1}^p \lambda_j}\]

であり、上位 $q$ 個の主成分の累積寄与率は$\sum_{k=1}^q \lambda_k / \sum_{j=1}^p \lambda_j$ で表される。慣例として累積寄与率が $80\%$ または $90\%$ に達する $q$ を次元削減の目標次元として選ぶことが多い。

最良低ランク近似(Eckart–Young–Mirsky 定理)

定理(Eckart–Young–Mirsky, 1936):行列 $\tilde{X} = UDV^\top$(SVD)に対して、フロベニウスノルムの意味での最良ランク $q$ 近似は

\[\tilde{X}_q = U_q D_q V_q^\top= \arg\min_{\mathrm{rank}(B) \leq q} \|\tilde{X} - B\|_F\]

で与えられ、近似誤差は

\[\|\tilde{X} - \tilde{X}_q\|_F^2 = \sum_{k=q+1}^r d_k^2 = n\sum_{k=q+1}^p \lambda_k\]

となる。この定理はスペクトルノルム(最大特異値ノルム)に対しても成立する。PCA による $q$ 次元近似は $\tilde{X}$ の最良ランク $q$ 近似であり、再構成誤差最小化との等価性が代数的に確認される。

主成分数の選択

スクリープロットと累積寄与率

固有値を降順にプロットしたスクリープロット(Scree Plot)において、固有値の変化が急峻から緩やかに変わる「肘点」(Elbow)で主成分数 $q$ を選ぶ方法が古典的に用いられる。視覚的判断に依存する難点があるが、Cattell(1966)のスクリーテストとして標準的な探索的解析手法として定着している。

平行分析

平行分析(Parallel Analysis;Horn, 1965)はデータと同一次元のランダム行列(標準正規エントリー)の固有値の分位点と観測固有値を比較し、ランダム構造を超える固有値を持つ成分のみを採用する。スクリープロットより客観的な基準を提供する。

情報量基準による選択

確率的 PCA(後述)の枠組みでは、次元 $q$ をモデルのパラメータ数と見て AIC・BIC による選択が可能である:

\[\mathrm{BIC}(q)= -2\ell(\hat{\boldsymbol{\theta}}_q) + d_q \log n\]

ここで $d_q = q(p - q) + q + 1$ は確率的 PCA のパラメータ数(ローディング行列の自由度 $q(p-q)$、固有値 $q$、ノイズ分散 $1$)であり、$\ell(\hat{\boldsymbol{\theta}}_q)$ は $q$ 次元確率的 PCA モデルの最大対数尤度である。Minka(2001)による BIC の解析的近似が実用的な基準として広く用いられる。

交差検証

主成分数 $q$ の交差検証による選択として、各観測を除いた LOOCV が理論的に適切である。Wold(1978)の手法では各行列要素を欠損値として扱い、PCA による補完誤差を $q$ の関数として最小化する。

確率的主成分分析(PPCA)

生成モデルの定義

確率的主成分分析(Probabilistic PCA, PPCA;Tipping–Bishop, 1999)はPCA に潜在変数モデルとしての確率的解釈を与える。潜在変数 $\boldsymbol{z} \in \mathbb{R}^q$($q < p$)と観測変数 $\boldsymbol{x} \in \mathbb{R}^p$ の生成過程を

\[\boldsymbol{z} \sim \mathcal{N}(\boldsymbol{0}, I_q),\qquad\boldsymbol{x} \mid \boldsymbol{z} \sim \mathcal{N}(W\boldsymbol{z} + \boldsymbol{\mu},\, \sigma^2 I_p)\]

と定義する($W \in \mathbb{R}^{p \times q}$:ローディング行列、$\sigma^2 > 0$:等方ノイズ)。観測変数の周辺分布は

\[\boldsymbol{x} \sim \mathcal{N}(\boldsymbol{\mu},\, WW^\top + \sigma^2 I_p)\]

であり、$\boldsymbol{\mu} = \bar{\boldsymbol{x}}$ と推定される。モデルパラメータは $\{W, \sigma^2\}$ であり、共分散行列 $C = WW^\top + \sigma^2 I_p$ が低ランク+等方ノイズの構造を持つ。

最尤推定と古典的 PCA との関係

対数尤度を最大化する MLE として

\[\hat{W} = V_q (D_q^2/n - \hat{\sigma}^2 I_q)^{1/2} R,\qquad\hat{\sigma}^2 = \frac{1}{p-q}\sum_{k=q+1}^p \lambda_k\]

が得られる($R \in \mathbb{R}^{q \times q}$:任意の直交行列、$V_q$:$S$ の上位 $q$ 固有ベクトル、$D_q^2/n = \mathrm{diag}(\lambda_1,\ldots,\lambda_q)$)。$\hat{W}$ の列空間は $S$ の上位 $q$ 固有ベクトルが張る空間と一致し、$\sigma^2 \to 0$($p - q$ 次元の残余固有値の平均がゼロ)の極限で古典的 PCA に一致する。$\hat{\sigma}^2$ は残余固有値の平均として、次元削減で捨てられる変動の等方的ノイズレベルを推定する。

EM アルゴリズムによる推定

PPCA の MLE は閉形式で得られるが(上式)、欠損データ・高次元設定での数値安定性からEM アルゴリズム(Expectation-Maximization)が広く用いられる。

E ステップ:現在のパラメータ $\{W^{(t)}, \sigma^{2(t)}\}$ のもとで潜在変数の事後分布を計算する:

\[\mathbb{E}[\boldsymbol{z}_i \mid \boldsymbol{x}_i]= M^{-1}W^{(t)\top}(\boldsymbol{x}_i - \bar{\boldsymbol{x}}),\quad\mathbb{E}[\boldsymbol{z}_i\boldsymbol{z}_i^\top \mid \boldsymbol{x}_i]= \sigma^{2(t)} M^{-1}+ \mathbb{E}[\boldsymbol{z}_i\mid\boldsymbol{x}_i]\mathbb{E}[\boldsymbol{z}_i\mid\boldsymbol{x}_i]^\top\]

ここで $M = W^{(t)\top}W^{(t)} + \sigma^{2(t)}I_q$。

M ステップ:期待完全対数尤度を最大化するパラメータ更新を行う:

\[W^{(t+1)} = \left[\sum_i (\boldsymbol{x}_i - \bar{\boldsymbol{x}})\mathbb{E}[\boldsymbol{z}_i\mid\boldsymbol{x}_i]^\top\right]\left[\sum_i \mathbb{E}[\boldsymbol{z}_i\boldsymbol{z}_i^\top\mid\boldsymbol{x}_i]\right]^{-1}\]\[\sigma^{2(t+1)} = \frac{1}{np}\sum_i \left[\|\boldsymbol{x}_i - \bar{\boldsymbol{x}}\|^2- 2\mathbb{E}[\boldsymbol{z}_i\mid\boldsymbol{x}_i]^\top W^{(t+1)\top}(\boldsymbol{x}_i - \bar{\boldsymbol{x}})+ \mathrm{tr}\!\left(\mathbb{E}[\boldsymbol{z}_i\boldsymbol{z}_i^\top\mid\boldsymbol{x}_i]W^{(t+1)\top}W^{(t+1)}\right)\right]\]

EM の各ステップは $O(npq)$ であり、$q \ll p$ の場合に共分散行列の固有値分解($O(p^3)$)より効率的かつ数値安定である。欠損データへの自然な拡張も容易である。

正則化と堅牢 PCA

正則化 PCA

高次元小標本設定($p \gg n$)では標本共分散行列 $S$ の固有値推定に深刻なバイアスが生じる(固有値の分散汚染)。Marchenko–Pastur 則によれば、$p/n \to \gamma \in (0,\infty)$ のとき$S$ の固有値は真の固有値から系統的に歪む。対処として以下の正則化が用いられる。

堅牢 PCA(RPCA)

外れ値や大きなスパースノイズに対してロバストな PCA として、堅牢 PCA(Robust PCA, RPCA;Candès–Li–Ma–Wright, 2011)は観測行列 $M \in \mathbb{R}^{n \times p}$ を低ランク行列 $L$ とスパース行列 $S$ の和に分解する:

\[M = L + S,\quad \mathrm{rank}(L) \ll \min(n,p),\quad \|S\|_0 \ll np\]

凸緩和として主成分追跡(Principal Component Pursuit, PCP)問題

\[\min_{L, S} \|L\|_* + \lambda\|S\|_1\quad \text{s.t.} \quad M = L + S\]

を解く($\|\cdot\|_*$:核ノルム=特異値の和、$\|\cdot\|_1$:要素和絶対値)。核ノルムは行列のランクの凸緩和であり、$\ell_1$ ノルムはスパースノイズの凸緩和として機能する。$\lambda = 1/\sqrt{\max(n,p)}$ と選ぶとき、$L$ と $S$ の非コヒーレント条件(インコヒーレンス条件)のもとで、PCP は確率 $1$ で $L^*$ と $S^*$ を正確に回復する(Candès et al., 2011)。凸問題であるため加速近接勾配法(FISTA)や交互方向乗数法(ADMM)による効率的な解法が知られている。

スパース PCA

古典的 PCA のローディングは一般にすべての変数に非ゼロの係数を持つため、高次元設定での解釈が困難である。スパース PCA(Zou–Hastie–Tibshirani, 2006)は各ローディングのスパース性を課した最適化問題

\[\min_{A, B}\|\tilde{X} - \tilde{X}BA^\top\|_F^2+ \lambda\sum_{k=1}^q \|\boldsymbol{b}_k\|_1\quad \text{s.t.} \quad A^\top A = I_q\]

を解く($A \in \mathbb{R}^{p \times q}$:直交ローディング行列、$B \in \mathbb{R}^{p \times q}$:スパースな近似ローディング)。$\ell_1$ ペナルティ $\lambda\sum_k \|\boldsymbol{b}_k\|_1$ が各ローディングのスパース性を誘導する。アルゴリズムは $A$ を固定して $B$ を更新する(Lasso 問題)と$B$ を固定して $A$ を更新する(SVD 問題)を交互に繰り返す。各ローディングが少数の変数のみに非ゼロ係数を持つため、主成分の解釈可能性が大幅に向上する。

カーネル PCA

非線形な次元削減への拡張として、カーネル PCA(Schölkopf–Smola–Müller, 1998)は特徴写像 $\phi: \mathcal{X} \to \mathcal{F}$($\mathcal{F}$:RKHS)による高次元特徴空間での PCA をカーネルトリックにより実装する。

カーネル行列 $K \in \mathbb{R}^{n \times n}$($K_{ij} = k(\boldsymbol{x}_i, \boldsymbol{x}_j)$)を中心化した行列

\[\tilde{K}= K - \frac{1}{n}\boldsymbol{1}_n\boldsymbol{1}_n^\top K- K\frac{1}{n}\boldsymbol{1}_n\boldsymbol{1}_n^\top+ \frac{1}{n^2}\boldsymbol{1}_n\boldsymbol{1}_n^\top K \boldsymbol{1}_n\boldsymbol{1}_n^\top= \left(I - \frac{1}{n}\boldsymbol{1}_n\boldsymbol{1}_n^\top\right)K\left(I - \frac{1}{n}\boldsymbol{1}_n\boldsymbol{1}_n^\top\right)\]

の固有値分解 $\tilde{K} = \sum_{k=1}^n \mu_k \boldsymbol{u}_k \boldsymbol{u}_k^\top$を行い、上位 $q$ 個の固有ベクトルを用いて新たな観測 $\boldsymbol{x}$ の第 $k$ 主成分スコアを

\[z_k(\boldsymbol{x})= \sum_{i=1}^n \frac{u_{ki}}{\sqrt{\mu_k}}\tilde{k}(\boldsymbol{x}_i, \boldsymbol{x})\]

として計算する($\tilde{k}$ は中心化カーネル)。線形カーネル $k(\boldsymbol{x},\boldsymbol{x}') = \boldsymbol{x}^\top\boldsymbol{x}'$ は古典的 PCA に一致し、RBF カーネル $k(\boldsymbol{x},\boldsymbol{x}') = \exp(-\|\boldsymbol{x}-\boldsymbol{x}'\|^2/(2h^2))$ は非線形多様体構造を捉える。計算量は $O(n^3)$($\tilde{K}$ の固有値分解)であり、大規模データには Nyström 近似($O(nm^2)$、$m \ll n$)が用いられる。

PCA と他の次元削減手法の比較

手法最適化基準線形/非線形計算量特徴
PCA分散最大化/再構成誤差最小化線形$O(np\min(n,p))$大域的、直交、解釈容易
カーネル PCARKHS での分散最大化非線形$O(n^3)$カーネル選択が重要
LDAクラス間分散/クラス内分散比最大化線形$O(np^2 + p^3)$教師あり、高々 $K-1$ 次元
t-SNE高次元近傍構造の保存(KL 最小化)非線形$O(n^2)$(近似 $O(n\log n)$)可視化特化、新点写像不可
UMAPトポロジー的構造の保存非線形$O(n^{1.14})$(実用的)可視化・汎化両対応
オートエンコーダ再構成誤差最小化非線形アーキテクチャ依存表現力高、解釈困難
ICA統計的独立性の最大化線形$O(np^2)$非ガウス成分の分離

PCA の限界と注意点

まとめ

PCA は分散最大化・再構成誤差最小化・Eckart–Young–Mirsky の最良低ランク近似という三つの等価な定式化を持ち、SVD を通じて効率的に計算される。固有値の寄与率・スクリープロット・情報量基準・交差検証により次元 $q$ が選択される。確率的 PCA(PPCA)は潜在変数モデルとして正規性仮定のもとでMLE が古典的 PCA のローディングと一致することを示し、EM アルゴリズム・欠損データ・モデル選択への理論的基盤を提供する。堅牢 PCA は核ノルムと $\ell_1$ 正則化の凸緩和によりスパースノイズ存在下での低ランク回復を保証し、スパース PCA は $\ell_1$ ペナルティによりローディングの解釈可能性を高める。カーネル PCA は RKHS の表現定理を通じて非線形多様体構造への拡張を実現する。PCA は線形次元削減の基準点として、LDA・カーネル法・オートエンコーダ・多様体学習など現代の表現学習手法すべての理論的出発点に位置する。

Mathematics is the language with which God has written the universe.





















因子分析 判別分析(LDA/QDA) 正準相関分析(CCA) 指数型分布族 ロジスティック回帰