主成分分析：分散最大化はなぜ固有値問題になるのか

0. この記事の要点

身長・体重・胸囲・座高のように、たくさんの数値で表されたデータを、できるだけ情報を残したまま 2 つか 3 つの数値にまとめられないか、という素朴な疑問から始めます。ここでは「情報」を「個体どうしの散らばり」と決めると、いちばん散らばって見える向きを探す問題が、行列の固有値を求める問題にぴったり重なることを順に確かめます。読み終えると、データを何本の軸にまとめればどれだけ残るのかを、固有値から読めるようになります。

• 次元削減は「変数が多すぎて見えない・推定が不安定になる・計算が重い」という三つの困りごとへの対処です。主成分分析（principal component analysis, PCA）はそのうち線形かつ教師なしの方法です。
• PCA の設計原理はただ一つ、「射影したときの分散が最大になる方向を選ぶ」です。これを式に直すと、単位ベクトル $\boldsymbol{u}$ に対する二次形式（成分について 2 次の式）$\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u}$（$S$ は標本共分散行列）の最大化になります。
• 対称行列の二次形式を単位球面上で最大化する答えは、最大固有値と、それに属する固有ベクトルです。したがって第 1 主成分方向は $S$ の最大固有値の固有ベクトルになります（系 3.3）。
• 上位 $k$ 本の主成分が張る部分空間は、データ点からの垂直距離の二乗和を最小にする $k$ 次元アフィン部分空間（原点を通るとは限らない、平らな集合）でもあります。分散最大化と当てはめ誤差最小化は同じ問題の裏表です（定理 5.1）。
• 寄与率 $\lambda_k / \sum_j \lambda_j$ はその主成分が説明する全分散の割合で、$k$ 次元まで採ったときに捨てられる誤差はちょうど $\lambda_{k+1} + \cdots + \lambda_p$ です。
• PCA は変数の単位の取り方に依存し、クラスの区別が保たれる保証もありません。使う前に標準化と目的の確認が必要です（例 7.1、例 7.2）。

この記事で使う事前知識

• 固有値・固有ベクトル: 行列 $A$ をかけても向きが変わらないベクトルと、その伸び縮みの倍率のこと。→ 固有値と固有ベクトル
• スペクトル定理: 実対称行列は、互いに直交する固有ベクトルの基底をもつ、という定理。→ スペクトル定理
• 直交射影と正規直交基底: ベクトルを部分空間へ垂直に落とした影と、長さ 1 で互いに直交する基底のこと。→ 内積空間とグラム・シュミット直交化
• 対角化: 適切な基底をとると行列が対角行列で表せること。→ 対角化とジョルダン標準形
• 分散・共分散: データの散らばりの大きさと、2 つの変数が一緒に動く度合いのこと。→ 確率変数と期待値

1. 動機：変数が多いことの何が困るのか

機械学習で扱うデータは、たいてい 1 個体あたり多数の数値の組として与えられます。$28 \times 28$ 画素の手書き数字画像なら 784 個の数値、遺伝子発現データなら数万個、アンケートなら設問数だけの数値です。この「1 個体を表す数値の個数」を次元と呼び、以下では $p$ と書きます。

次元が高いと、少なくとも三つの困りごとが起きます。

第一に、見えません。人間が直接目で読める散布図は 2 次元、頑張って 3 次元までです。784 次元のデータについて「似ているものどうしが固まっているか」を確かめる方法が、そのままでは存在しません。

第二に、統計的に不利になります。標本数 $n$ に対して次元 $p$ が大きいと、推定すべき量の個数が標本数に対して過大になり、モデルは訓練データの偶然の凹凸まで覚えてしまいます。共分散行列の推定はその典型で、$p \times p$ の対称行列がもつ $p(p+1)/2$ 個の成分を $n$ 個の標本から決めることになります。後で見るように、$n - 1 < p$ のときは標本共分散行列の階数が $p$ に届かず、必ず退化します（注意 6.4）。

第三に、計算量です。多くのアルゴリズムは次元に対して線形以上の計算量をもちます。共分散行列を作るだけでも $O(np^2)$ 回の乗算が要ります。

一方で、実データの多くは「見かけの次元ほど自由度がない」という性質をもちます。身長・体重・胸囲・座高を測れば 4 次元のデータですが、これらは互いに強く相関しており、実質的には「体格の大きさ」というほぼ 1 本の軸で説明できてしまいます。手書き数字の画像でも、隣り合う画素の値はほとんど同じですから、784 個の数値が自由に動くわけではありません。見かけの次元と本質的な次元のこの差が、次元削減の付け入る隙です。

では「情報を落とさずに次元を減らす」とは何でしょうか。何を情報と呼ぶかを決めない限り、この問いは意味をもちません。PCA の答えは明快です — 情報とは**散らばり（分散）**である。

この立場を採る理由は、極端な場合を考えるとわかります。ある方向にデータを射影したとき、その値がすべての個体で同じ（分散が $0$）だったとしましょう。この座標は定数ですから、記録する必要がありません。捨てても個体を区別する能力は一切失われません。逆に、射影した値が大きくばらつく方向は、個体どうしを最もよく引き離す方向です。ならば「射影後の分散が大きい順に軸を採り、小さいものから捨てる」のが自然だろう、というのが PCA の設計原理です。

歴史的には二つの入口があります。Karl Pearson は 1901 年の論文で「空間内の点の集まりに最もよく当てはまる直線・平面」を、点から直線への垂直距離の二乗和を最小にする問題として論じました。Harold Hotelling は 1933 年に、心理測定の文脈で分散最大化の定式化を与え、principal component の名を与えています。当てはめ誤差の最小化と分散の最大化 — 出発点は違いますが、答えは同じ固有値問題になります。この一致は偶然ではなく、ピタゴラスの定理から導かれる必然です（定理 5.1）。

前章までで扱った線形回帰や勾配降下法は、目的変数 $y$ という「正解」をもつ教師あり学習でした。PCA は $y$ を使いません。入力データ $\boldsymbol{x}$ 自身の構造だけを見る教師なし学習です。この違いは後で効いてきます（例 7.2）。

2. 準備：中心化と標本共分散行列

定義 2.1（中心化データ行列と標本共分散行列）↗

$\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n \in \mathbb{R}^p$（$n \ge 2$）を観測データとする。標本平均（重心）を

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

とおき、中心化ベクトルを $\boldsymbol{y}_i = \boldsymbol{x}_i - \bar{\boldsymbol{x}}$ と定める。$\boldsymbol{y}_i^{\mathsf{T}}$ を第 $i$ 行にもつ行列 $X \in \mathbb{R}^{n \times p}$ を中心化データ行列と呼ぶ。さらに

$$S = \frac{1}{n} X^{\mathsf{T}} X = \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{y}_i \boldsymbol{y}_i^{\mathsf{T}} \in \mathbb{R}^{p \times p}$$

を標本共分散行列と呼ぶ。その $(j, l)$ 成分は $s_{jl} = \frac{1}{n}\sum_{i=1}^n y_{ij} y_{il}$ であり、$s_{jj}$ は第 $j$ 変数の標本分散、$s_{jl}$ は第 $j$ 変数と第 $l$ 変数の標本共分散である。

定義の中で $\frac{1}{n} X^{\mathsf{T}} X$ と $\frac{1}{n}\sum_i \boldsymbol{y}_i \boldsymbol{y}_i^{\mathsf{T}}$ を同一視しました。これは行列積の定義から出ます。$X^{\mathsf{T}}$ の第 $i$ 列が $\boldsymbol{y}_i$ なので、$(X^{\mathsf{T}} X)_{jl} = \sum_{i=1}^{n} X_{ij} X_{il} = \sum_{i=1}^{n} y_{ij} y_{il} = \left(\sum_i \boldsymbol{y}_i \boldsymbol{y}_i^{\mathsf{T}}\right)_{jl}$ となります。

注意 2.2

不偏分散を使う流儀では $\frac{1}{n-1}$ で割ります。この行列を $S'$ と書くと $S' = \frac{n}{n-1} S$ のように両者は正の定数倍しか違わないので、固有ベクトル（＝主成分方向）は完全に一致し、固有値は同じ定数倍だけずれます。寄与率は固有値の比なので、これも一致します。以下では表記が軽くなる $\frac{1}{n}$ を採ります。

PCA が見るのはただ一つの量、「ある方向へ射影したときの分散」です。これを定義します。

定義 2.3（方向への射影と射影分散）↗

$\boldsymbol{u} \in \mathbb{R}^p$ を $\|\boldsymbol{u}\| = 1$ を満たすベクトルとする。第 $i$ 個体の $\boldsymbol{u}$ 方向の得点を

$$z_i(\boldsymbol{u}) = \langle \boldsymbol{u}, \boldsymbol{y}_i \rangle = \boldsymbol{u}^{\mathsf{T}}(\boldsymbol{x}_i - \bar{\boldsymbol{x}})$$

と定め、その標本分散を $\boldsymbol{u}$ 方向の射影分散と呼ぶ。

得点の平均は $\frac{1}{n}\sum_i z_i(\boldsymbol{u}) = \boldsymbol{u}^{\mathsf{T}}\left(\frac{1}{n}\sum_i \boldsymbol{y}_i\right) = \boldsymbol{u}^{\mathsf{T}}\boldsymbol{0} = 0$ です（中心化しているので $\sum_i \boldsymbol{y}_i = \boldsymbol{0}$）。したがって射影分散は二乗平均そのもので、

$$\frac{1}{n}\sum_{i=1}^{n} z_i(\boldsymbol{u})^2 = \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{u}^{\mathsf{T}} \boldsymbol{y}_i \, \boldsymbol{y}_i^{\mathsf{T}} \boldsymbol{u} = \boldsymbol{u}^{\mathsf{T}} \left( \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{y}_i \boldsymbol{y}_i^{\mathsf{T}} \right) \boldsymbol{u} = \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u}$$

と書けます。途中で使ったのは、スカラー $z_i = \boldsymbol{u}^{\mathsf{T}}\boldsymbol{y}_i$ について $z_i^2 = z_i z_i^{\mathsf{T}} = (\boldsymbol{u}^{\mathsf{T}}\boldsymbol{y}_i)(\boldsymbol{y}_i^{\mathsf{T}}\boldsymbol{u})$ という書き換えと、$\boldsymbol{u}$ が $i$ に依存しないので和の外に出せることだけです。この式が記事全体の出発点です。データの散らばりを方向 $\boldsymbol{u}$ ごとに測る関数が、行列 $S$ の二次形式として書けてしまう。これで問題は線形代数に移りました。

$S$ がどんな行列なのかを押さえておきます。

補題 2.4（標本共分散行列は対称かつ半正定値）✓

定義 2.1 の $S$ は $S^{\mathsf{T}} = S$ を満たし、さらに任意の $\boldsymbol{u} \in \mathbb{R}^p$ に対して $\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} \ge 0$ である。

証明（補題 2.4）

対称性: 転置の性質 $(AB)^{\mathsf{T}} = B^{\mathsf{T}}A^{\mathsf{T}}$ と $(A^{\mathsf{T}})^{\mathsf{T}} = A$ から

$$S^{\mathsf{T}} = \left(\frac{1}{n}X^{\mathsf{T}}X\right)^{\mathsf{T}} = \frac{1}{n}X^{\mathsf{T}}(X^{\mathsf{T}})^{\mathsf{T}} = \frac{1}{n}X^{\mathsf{T}}X = S .$$

半正定値性: 任意の $\boldsymbol{u}$ に対し

$$\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \frac{1}{n}\boldsymbol{u}^{\mathsf{T}}X^{\mathsf{T}}X\boldsymbol{u} = \frac{1}{n}(X\boldsymbol{u})^{\mathsf{T}}(X\boldsymbol{u}) = \frac{1}{n}\|X\boldsymbol{u}\|^2 \ge 0 .$$

ノルムの二乗は非負なので結論を得ます。なお $\|\boldsymbol{u}\| = 1$ のときこの量は射影分散に等しく、「分散は負にならない」という当たり前の事実の言い換えになっています。

∎

対称行列については、次の定理が使えます。これが PCA の理論的な土台です。

定理 2.5（実対称行列のスペクトル定理）✓

$S \in \mathbb{R}^{p \times p}$ が $S^{\mathsf{T}} = S$ を満たすとする。このとき $\mathbb{R}^p$ の正規直交基底 $\boldsymbol{v}_1, \ldots, \boldsymbol{v}_p$ と実数 $\lambda_1 \ge \lambda_2 \ge \cdots \ge \lambda_p$ が存在して、すべての $j$ について $S\boldsymbol{v}_j = \lambda_j \boldsymbol{v}_j$ が成り立つ。行列で書けば、$V = [\boldsymbol{v}_1\ \cdots\ \boldsymbol{v}_p]$、$\Lambda = \operatorname{diag}(\lambda_1, \ldots, \lambda_p)$ として $V^{\mathsf{T}}V = I_p$ かつ $S = V \Lambda V^{\mathsf{T}}$ である。さらに $S$ が半正定値ならば $\lambda_p \ge 0$、すなわちすべての固有値は非負である。

証明（定理 2.5）

前半（正規直交固有基底の存在と固有値の実数性）は本記事では認めて使います。所在は 注意 2.6 に示します。ここでは後半だけ示します。$S$ が半正定値のとき、$\boldsymbol{v}_j$ は単位ベクトルなので

$$\lambda_j = \lambda_j \|\boldsymbol{v}_j\|^2 = \lambda_j \boldsymbol{v}_j^{\mathsf{T}}\boldsymbol{v}_j = \boldsymbol{v}_j^{\mathsf{T}}(\lambda_j \boldsymbol{v}_j) = \boldsymbol{v}_j^{\mathsf{T}} S \boldsymbol{v}_j \ge 0$$

となります。最後の不等号が半正定値性です。とくに 補題 2.4 より標本共分散行列の固有値はすべて非負です。

∎

注意 2.6

スペクトル定理の証明は スペクトル定理 にあります。ここで使う実対称行列の場合は 系 4.3［スペクトル定理］ がそのままの形です。固有値・固有ベクトルの定義と特性多項式による計算は 固有値と固有ベクトル（定義 3.1［固有値と固有ベクトル］）、直交射影と正規直交基底の扱いは 内積空間とグラム・シュミット直交化（定理 7.1［内積空間とグラム・シュミット直交化］）を参照してください。以下では $\boldsymbol{v}_1, \ldots, \boldsymbol{v}_p$ と $\lambda_1 \ge \cdots \ge \lambda_p \ge 0$ を、この定理が与える標本共分散行列 $S$ の正規直交固有系と固有値（降順）として固定します。

3. 第 1 主成分：分散最大化がレイリー商になる

PCA が解きたい問題は、§1 の設計原理をそのまま式にしたものです。

$$\text{最大化} \quad \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} \quad \text{（条件 } \|\boldsymbol{u}\| = 1 \text{ のもとで）}$$

なぜ $\|\boldsymbol{u}\| = 1$ という条件を付けるのでしょうか。付けないと問題が壊れるからです。$\boldsymbol{u}$ を $c$ 倍すると得点も $c$ 倍になり、射影分散は $(c\boldsymbol{u})^{\mathsf{T}}S(c\boldsymbol{u}) = c^2 \boldsymbol{u}^{\mathsf{T}}S\boldsymbol{u}$ と $c^2$ 倍になります。$S \ne O$ なら $c \to \infty$ で発散し、最大値が存在しません。私たちが知りたいのは「どの向きにデータが伸びているか」であって、物差しの長さではありません。長さを $1$ に固定すれば、比較しているのは向きだけになります。

定理 3.1（レイリー商の最大値）✓

$S \in \mathbb{R}^{p\times p}$ を実対称行列とし、定理 2.5 の記号を使う。このとき

$$\max_{\|\boldsymbol{u}\| = 1} \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \lambda_1$$

が成り立ち、最大値は $\boldsymbol{u} = \boldsymbol{v}_1$ で達成される。さらに、最大値を達成する単位ベクトルの全体は、固有値 $\lambda_1$ の固有空間 $\ker(S - \lambda_1 I)$ に含まれる単位ベクトルの全体と一致する。

証明の方針は次のとおりです。$S$ の固有ベクトルを新しい座標軸に選ぶと、二次形式は「固有値を、重み（成分の二乗）つきで平均したもの」に書き直せます。重みの合計は 1 なので、平均が最大固有値を超えることはなく、全部の重みを最大固有値の軸に置いたときにちょうど等しくなります。

証明（定理 3.1）

ステップ 1（成分表示）. 定理 2.5 より $\boldsymbol{v}_1, \ldots, \boldsymbol{v}_p$ は $\mathbb{R}^p$ の正規直交基底なので、任意の $\boldsymbol{u}$ は $\boldsymbol{u} = \sum_{j=1}^p c_j \boldsymbol{v}_j$（$c_j = \langle \boldsymbol{v}_j, \boldsymbol{u}\rangle$）と一意に書けます。正規直交性から

$$\|\boldsymbol{u}\|^2 = \left\langle \sum_j c_j \boldsymbol{v}_j, \sum_l c_l \boldsymbol{v}_l \right\rangle = \sum_{j,l} c_j c_l \langle \boldsymbol{v}_j, \boldsymbol{v}_l\rangle = \sum_{j=1}^{p} c_j^2$$

なので、条件 $\|\boldsymbol{u}\| = 1$ は $\sum_j c_j^2 = 1$ と同値です。

ステップ 2（二次形式の計算）. $S\boldsymbol{u} = \sum_j c_j S\boldsymbol{v}_j = \sum_j c_j \lambda_j \boldsymbol{v}_j$ なので、同じく正規直交性から

$$\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \left\langle \sum_l c_l \boldsymbol{v}_l, \sum_j c_j \lambda_j \boldsymbol{v}_j \right\rangle = \sum_{j=1}^{p} \lambda_j c_j^2 .$$

つまり二次形式は、固有値を重み $c_j^2$ で平均したものです（重みの和はステップ 1 より $1$）。

ステップ 3（上からの評価）. すべての $j$ で $\lambda_j \le \lambda_1$ かつ $c_j^2 \ge 0$ なので $\lambda_j c_j^2 \le \lambda_1 c_j^2$、したがって

$$\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \sum_j \lambda_j c_j^2 \le \lambda_1 \sum_j c_j^2 = \lambda_1 .$$

ステップ 4（達成）. $\boldsymbol{u} = \boldsymbol{v}_1$ は単位ベクトルであり、$\boldsymbol{v}_1^{\mathsf{T}} S \boldsymbol{v}_1 = \boldsymbol{v}_1^{\mathsf{T}}(\lambda_1 \boldsymbol{v}_1) = \lambda_1$。よって上限 $\lambda_1$ は実際に達成され、最大値です。

ステップ 5（等号成立の条件）. ステップ 3 の不等式の差を取ると

$$\lambda_1 - \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \sum_{j=1}^{p} (\lambda_1 - \lambda_j) c_j^2$$

であり、右辺は非負項の和です。したがって等号が成り立つのは、$\lambda_j < \lambda_1$ であるすべての $j$ について $c_j = 0$ となるとき、かつそのときに限ります。これは $\boldsymbol{u}$ が $\lambda_1$ に属する固有ベクトルたちの張る空間 $\ker(S - \lambda_1 I) = \operatorname{span}\{\boldsymbol{v}_j : \lambda_j = \lambda_1\}$ に入ることと同値です。

∎

この定理をデータの言葉に翻訳すれば、第 1 主成分が得られます。

定義 3.2（主成分方向と主成分得点）↗

定理 2.5 の正規直交固有系 $\boldsymbol{v}_1, \ldots, \boldsymbol{v}_p$（固有値の降順）に対し、$\boldsymbol{v}_k$ を第 $k$ 主成分方向と呼ぶ。第 $i$ 個体の第 $k$ 主成分得点を

$$z_{ik} = \langle \boldsymbol{v}_k, \boldsymbol{y}_i \rangle = \boldsymbol{v}_k^{\mathsf{T}}(\boldsymbol{x}_i - \bar{\boldsymbol{x}})$$

と定める。第 $k$ 主成分得点を並べた $n$ 個の数の組を第 $k$ 主成分と呼ぶこともある。

系 3.3（第 1 主成分は最大固有値の固有ベクトル）✓

すべての単位ベクトル $\boldsymbol{u}$ のうち、射影分散 $\frac{1}{n}\sum_i z_i(\boldsymbol{u})^2$ を最大にするのは $S$ の最大固有値 $\lambda_1$ に属する単位固有ベクトルであり、そのときの射影分散の値は $\lambda_1$ に等しい。

証明（系 3.3）

§2 で示したとおり射影分散は $\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u}$ に等しく、補題 2.4 より $S$ は実対称なので 定理 3.1 がそのまま適用できます。最大値は $\lambda_1$、達成するのは $\lambda_1$ の固有空間内の単位ベクトルです。

∎

注意 3.4

条件付き最大化ですから、ラグランジュの未定乗数法でも同じ結論に達します。$L(\boldsymbol{u}, \mu) = \boldsymbol{u}^{\mathsf{T}}S\boldsymbol{u} - \mu(\boldsymbol{u}^{\mathsf{T}}\boldsymbol{u} - 1)$ とおきます。$S$ が対称であることを使うと

$$\frac{\partial}{\partial u_k}\sum_{j,l} S_{jl}u_j u_l = \sum_{l} S_{kl}u_l + \sum_{j} S_{jk}u_j = 2 (S\boldsymbol{u})_k$$

なので、$\nabla_{\boldsymbol{u}} L = 2S\boldsymbol{u} - 2\mu\boldsymbol{u} = \boldsymbol{0}$、すなわち $S\boldsymbol{u} = \mu\boldsymbol{u}$ が停留条件です。制約付き最適化の停留点は、そのまま固有方程式になります。しかも停留点では $\boldsymbol{u}^{\mathsf{T}}S\boldsymbol{u} = \mu\boldsymbol{u}^{\mathsf{T}}\boldsymbol{u} = \mu$ なので、目的関数の値は乗数そのものです。よって最大を与えるのは最大固有値。ただしこの議論だけでは「最大値が存在すること」を別途言う必要があります（単位球面はコンパクト、二次形式は連続なので最大値は存在します）。定理 3.1 の証明はこの点も込みで初等的に片付いています。

注意 3.5

主成分方向には二つの不定性があります。第一に符号です。$\boldsymbol{v}$ が単位固有ベクトルなら $-\boldsymbol{v}$ もそうで、得点の符号がすべて反転するだけですから、両者は同じ主成分を表します。数値ライブラリの出力の符号が実行環境によって変わることがあるのはこのためで、異常ではありません。第二に、固有値が重複すると固有空間が 2 次元以上になり、その中のどの正規直交基底を採っても構わなくなります（例 3.6）。「第 1 主成分方向は $\boldsymbol{u}$ である」と一意に言えるのは $\lambda_1 > \lambda_2$ のときだけです。

例 3.6（等方的なデータでは主成分方向が定まらない）

$S = \sigma^2 I_p$（$\sigma > 0$）となるデータを考えます。どの方向を選んでも射影分散は $\boldsymbol{u}^{\mathsf{T}}(\sigma^2 I)\boldsymbol{u} = \sigma^2\|\boldsymbol{u}\|^2 = \sigma^2$ で同じです。固有値は $\sigma^2$ のみ（重複度 $p$）で、任意の正規直交基底が主成分系になります。$k$ 次元に落とすと、どの部分空間を選んでも $(p-k)\sigma^2$ の分散を失います。次元削減で得をする余地がまったくない状況です。実データでも上位の固有値が接近していると、標本のわずかな揺らぎで主成分方向が大きく回転します。「第 2 主成分は $\lambda_2 \approx \lambda_3$ のとき解釈しない」というのは実務上の鉄則です。

4. 第 k 主成分と主成分得点の性質

第 1 主成分だけでは足りないとき、次の軸をどう選ぶか。分散が大きい方向をもう一度探すと、また $\boldsymbol{v}_1$ の近くが選ばれてしまい、同じ情報を二度記録することになります。そこで「すでに採った方向と直交する」という条件を課します。直交する方向の得点は、後で見るように互いに無相関になり、情報の重複が起きません。

定理 4.1（逐次的な分散最大化と固有値の対応）✓

定理 2.5 の記号のもとで、$1 \le k \le p$ とする。集合

$$W_k = \{\boldsymbol{u} \in \mathbb{R}^p : \|\boldsymbol{u}\| = 1,\ \langle \boldsymbol{u}, \boldsymbol{v}_1\rangle = \cdots = \langle \boldsymbol{u}, \boldsymbol{v}_{k-1}\rangle = 0 \}$$

（$k = 1$ のときは単位球面全体）の上で、$\boldsymbol{u}^{\mathsf{T}}S\boldsymbol{u}$ の最大値は $\lambda_k$ であり、$\boldsymbol{u} = \boldsymbol{v}_k$ で達成される。

証明（定理 4.1）

$\boldsymbol{u} \in W_k$ を取り、定理 3.1 の証明のステップ 1 と同様に $\boldsymbol{u} = \sum_{j=1}^p c_j \boldsymbol{v}_j$、$c_j = \langle \boldsymbol{v}_j, \boldsymbol{u}\rangle$ と展開します。$W_k$ の定義から $c_1 = \cdots = c_{k-1} = 0$ です。したがって $\sum_{j \ge k} c_j^2 = \|\boldsymbol{u}\|^2 = 1$ であり、ステップ 2 の計算から

$$\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \sum_{j = k}^{p} \lambda_j c_j^2 \le \lambda_k \sum_{j=k}^{p} c_j^2 = \lambda_k$$

を得ます。不等号では、$j \ge k$ のとき固有値が降順に並んでいることから $\lambda_j \le \lambda_k$ であることを使いました。

一方 $\boldsymbol{v}_k$ は単位ベクトルで、$j \le k-1$ について $\langle \boldsymbol{v}_k, \boldsymbol{v}_j \rangle = 0$（正規直交性）なので $\boldsymbol{v}_k \in W_k$ であり、$\boldsymbol{v}_k^{\mathsf{T}}S\boldsymbol{v}_k = \lambda_k$。よって上限は達成され、最大値です。

∎

つまり「分散が最大の方向を採る → それに直交する中で分散が最大の方向を採る → …」という素朴な逐次手続きが、固有値を大きい順に並べる操作とぴったり一致します。定義 3.2 で固有値の降順に番号を付けたのは、この一致を保つためです。

命題 4.2（主成分得点の性質）✓

定義 3.2 の主成分得点 $z_{ik}$ について、次が成り立つ。

1. 各 $k$ について $\frac{1}{n}\sum_{i=1}^{n} z_{ik} = 0$。
2. 各 $k, l$ について $\frac{1}{n}\sum_{i=1}^{n} z_{ik} z_{il} = \lambda_k \delta_{kl}$。すなわち第 $k$ 主成分得点の標本分散は $\lambda_k$ であり、異なる主成分得点どうしの標本共分散は $0$ である。
3. $\displaystyle\sum_{k=1}^{p} \lambda_k = \operatorname{tr} S = \sum_{j=1}^{p} s_{jj} = \frac{1}{n}\sum_{i=1}^{n}\|\boldsymbol{y}_i\|^2$。

証明（命題 4.2）

1. 内積の線形性から $\frac{1}{n}\sum_i z_{ik} = \left\langle \boldsymbol{v}_k, \frac{1}{n}\sum_i \boldsymbol{y}_i \right\rangle = \langle \boldsymbol{v}_k, \boldsymbol{0}\rangle = 0$。ここで 定義 2.1 の中心化から $\sum_i \boldsymbol{y}_i = \boldsymbol{0}$ を使いました。

2. $z_{ik} z_{il} = (\boldsymbol{v}_k^{\mathsf{T}}\boldsymbol{y}_i)(\boldsymbol{y}_i^{\mathsf{T}}\boldsymbol{v}_l)$ と書けるので

$$\frac{1}{n}\sum_i z_{ik} z_{il} = \boldsymbol{v}_k^{\mathsf{T}}\left(\frac{1}{n}\sum_i \boldsymbol{y}_i\boldsymbol{y}_i^{\mathsf{T}}\right)\boldsymbol{v}_l = \boldsymbol{v}_k^{\mathsf{T}} S \boldsymbol{v}_l = \lambda_l \boldsymbol{v}_k^{\mathsf{T}}\boldsymbol{v}_l = \lambda_l \delta_{kl},$$

最後は正規直交性です。$k = l$ なら $\lambda_k$、$k \ne l$ なら $0$。1. より得点の平均が $0$ なので、この量はそのまま標本共分散です。

3. まず $\operatorname{tr}(\boldsymbol{y}\boldsymbol{y}^{\mathsf{T}}) = \sum_j y_j^2 = \|\boldsymbol{y}\|^2$ とトレースの線形性から

$$\operatorname{tr} S = \operatorname{tr}\left(\frac{1}{n}\sum_i \boldsymbol{y}_i\boldsymbol{y}_i^{\mathsf{T}}\right) = \frac{1}{n}\sum_i \|\boldsymbol{y}_i\|^2 .$$

また $\operatorname{tr}S = \sum_j s_{jj}$ はトレースの定義そのものです。最後に 定理 2.5 の $S = V\Lambda V^{\mathsf{T}}$ と $\operatorname{tr}(AB) = \operatorname{tr}(BA)$、$V^{\mathsf{T}}V = I$ から

$$\operatorname{tr}S = \operatorname{tr}(V\Lambda V^{\mathsf{T}}) = \operatorname{tr}(\Lambda V^{\mathsf{T}}V) = \operatorname{tr}\Lambda = \sum_k \lambda_k .$$∎

2. は PCA のもう一つの顔です。もとの変数は互いに相関していますが、主成分得点に移ると相関が消えます。$V^{\mathsf{T}}$ による座標変換は、共分散行列を対角化する変換にほかなりません（対角化とジョルダン標準形、定理 3.2［対角化とジョルダン標準形］）。3. は「全分散はどの座標系で測っても同じで、それが固有値の総和に等しい」という保存則で、次節の寄与率の分母になります。

5. もう一つの顔：当てはめ誤差の最小化

Pearson の出発点は分散ではなく距離でした。ここでは、その定式化が同じ固有ベクトルに行き着くことを示します。$V \subseteq \mathbb{R}^p$ を線形部分空間、$P_V$ をそれへの直交射影とします。

定理 5.1（最良近似アフィン部分空間）✓

$1 \le k \le p$ とする。$k$ 次元線形部分空間 $V \subseteq \mathbb{R}^p$ と点 $\boldsymbol{c} \in \mathbb{R}^p$ が定めるアフィン部分空間 $A = \boldsymbol{c} + V$ に対し、当てはめ誤差を

$$E(\boldsymbol{c}, V) = \frac{1}{n}\sum_{i=1}^{n} \operatorname{dist}(\boldsymbol{x}_i, A)^2 = \frac{1}{n}\sum_{i=1}^{n} \left\| (I - P_V)(\boldsymbol{x}_i - \boldsymbol{c}) \right\|^2$$

と定める。このとき

$$\min_{\boldsymbol{c} \in \mathbb{R}^p,\ \dim V = k} E(\boldsymbol{c}, V) = \lambda_{k+1} + \lambda_{k+2} + \cdots + \lambda_p$$

であり、最小値は $\boldsymbol{c} = \bar{\boldsymbol{x}}$、$V = \operatorname{span}\{\boldsymbol{v}_1, \ldots, \boldsymbol{v}_k\}$ で達成される（$k = p$ のとき右辺は $0$ と読む）。

この証明は 4 段階です。まず、直線や平面を重心に通すのが最良であることを示し（ステップ 1）、次にピタゴラスの定理で「誤差を小さくする」ことを「射影した長さの二乗和を大きくする」ことに言い換えます（ステップ 2。下の図のように、1 点の長さの二乗は、影の長さの二乗と、影から点までの距離の二乗に分かれます）。ステップ 3 では、その和が上位 $k$ 個の固有値の和を超えないことを、各固有値にかかる重み $w_j$ が 0 以上 1 以下で合計が $k$ であることから示します。固有値の大きいほうから重みを詰めるのが最良、という見通しです。

ステップ 2 のピタゴラスの三角形。点 y（原点からの位置ベクトル）を直線 V に垂直に落とした影が P_V y です。長さの二乗は、影の部分と、影から点までの垂直な部分に分かれます。

証明（定理 5.1）

ステップ 0（距離の式）. $A = \boldsymbol{c} + V$ 上の点は $\boldsymbol{c} + \boldsymbol{w}$（$\boldsymbol{w} \in V$）と書けます。$\|\boldsymbol{x} - \boldsymbol{c} - \boldsymbol{w}\|$ を $\boldsymbol{w} \in V$ について最小にするのは直交射影 $\boldsymbol{w} = P_V(\boldsymbol{x} - \boldsymbol{c})$ であり（定理 7.2［内積空間とグラム・シュミット直交化］）、そのときの残差が $(I - P_V)(\boldsymbol{x}-\boldsymbol{c})$ です。よって定義式の二つの表現は一致します。

ステップ 1（重心が最適であること）. $\boldsymbol{x}_i - \boldsymbol{c} = \boldsymbol{y}_i + (\bar{\boldsymbol{x}} - \boldsymbol{c})$ と分解し、$\boldsymbol{d} = (I-P_V)(\bar{\boldsymbol{x}} - \boldsymbol{c})$ とおくと $(I-P_V)(\boldsymbol{x}_i - \boldsymbol{c}) = (I-P_V)\boldsymbol{y}_i + \boldsymbol{d}$。二乗して和を取ると

$$n E(\boldsymbol{c}, V) = \sum_i \|(I-P_V)\boldsymbol{y}_i\|^2 + 2\left\langle \sum_i (I-P_V)\boldsymbol{y}_i,\ \boldsymbol{d} \right\rangle + n\|\boldsymbol{d}\|^2 .$$

中央の項は、$\sum_i \boldsymbol{y}_i = \boldsymbol{0}$ と $I - P_V$ の線形性から $\sum_i (I-P_V)\boldsymbol{y}_i = (I-P_V)\boldsymbol{0} = \boldsymbol{0}$ となって消えます。残りは $\|\boldsymbol{d}\|^2 \ge 0$ なので、$\boldsymbol{c} = \bar{\boldsymbol{x}}$（このとき $\boldsymbol{d} = \boldsymbol{0}$）が最適です。以下 $\boldsymbol{c} = \bar{\boldsymbol{x}}$ とし、$E(V) = \frac{1}{n}\sum_i \|(I-P_V)\boldsymbol{y}_i\|^2$ と書きます。データを中心化してよい理由がここにあります。

ステップ 2（ピタゴラス：誤差最小化＝分散最大化）. 直交射影の性質から $P_V \boldsymbol{y} \perp (I - P_V)\boldsymbol{y}$ なので、

$$\|\boldsymbol{y}_i\|^2 = \|P_V \boldsymbol{y}_i\|^2 + \|(I-P_V)\boldsymbol{y}_i\|^2 .$$

$i$ について平均を取り 命題 4.2 の 3. を使うと

$$E(V) = \operatorname{tr}S - \frac{1}{n}\sum_i \|P_V\boldsymbol{y}_i\|^2 .$$

$\operatorname{tr}S$ は $V$ に依らない定数なので、$E(V)$ の最小化は射影後の分散の和 $\frac{1}{n}\sum_i\|P_V\boldsymbol{y}_i\|^2$ の最大化と同値です。Pearson の問題と Hotelling の問題が一致するのは、このピタゴラスの定理 1 本によります。

ステップ 3（最大化）. $\boldsymbol{u}_1, \ldots, \boldsymbol{u}_k$ を $V$ の正規直交基底とすると $P_V\boldsymbol{y} = \sum_{m=1}^k \langle \boldsymbol{u}_m, \boldsymbol{y}\rangle \boldsymbol{u}_m$、したがって $\|P_V\boldsymbol{y}\|^2 = \sum_{m=1}^k \langle \boldsymbol{u}_m,\boldsymbol{y}\rangle^2$ です。§2 の射影分散の計算を各 $\boldsymbol{u}_m$ に適用すると

$$\frac{1}{n}\sum_i \|P_V\boldsymbol{y}_i\|^2 = \sum_{m=1}^{k} \boldsymbol{u}_m^{\mathsf{T}} S \boldsymbol{u}_m .$$

これを $\lambda_1 + \cdots + \lambda_k$ で上から抑えます。$w_j = \sum_{m=1}^{k}\langle \boldsymbol{v}_j, \boldsymbol{u}_m\rangle^2$ とおくと、定理 3.1 の証明ステップ 2 と同じ展開から $\sum_m \boldsymbol{u}_m^{\mathsf{T}}S\boldsymbol{u}_m = \sum_{j=1}^{p}\lambda_j w_j$ です。重み $w_j$ は二つの制約を満たします。

• $0 \le w_j \le 1$：$w_j = \|P_V \boldsymbol{v}_j\|^2$ であり、$\|P_V\boldsymbol{v}_j\| \le \|\boldsymbol{v}_j\| = 1$（ベッセルの不等式（注意 7.3）［内積空間とグラム・シュミット直交化］、あるいはステップ 2 のピタゴラス）。
• $\sum_{j=1}^p w_j = k$：和の順序を変えると $\sum_j w_j = \sum_{m=1}^{k}\sum_{j=1}^{p}\langle \boldsymbol{v}_j,\boldsymbol{u}_m\rangle^2 = \sum_{m=1}^k \|\boldsymbol{u}_m\|^2 = k$（$\boldsymbol{v}_j$ が正規直交基底であることによるパーセバルの等式）。

このとき

$$\sum_{j=1}^{p}\lambda_j w_j - \sum_{j=1}^{k}\lambda_j = \sum_{j\le k}\lambda_j (w_j - 1) + \sum_{j > k}\lambda_j w_j \le \sum_{j\le k}\lambda_k (w_j - 1) + \sum_{j>k}\lambda_k w_j = \lambda_k\left(\sum_{j=1}^p w_j - k\right) = 0 .$$

不等号の根拠は次のとおりです。$j \le k$ では $w_j - 1 \le 0$ かつ $\lambda_j \ge \lambda_k$ なので $\lambda_j(w_j-1) \le \lambda_k(w_j-1)$。$j > k$ では $w_j \ge 0$ かつ $\lambda_j \le \lambda_k$ なので $\lambda_j w_j \le \lambda_k w_j$。したがって $\sum_m \boldsymbol{u}_m^{\mathsf{T}}S\boldsymbol{u}_m \le \lambda_1 + \cdots + \lambda_k$ です。

ステップ 4（達成と結論）. $V = \operatorname{span}\{\boldsymbol{v}_1,\ldots,\boldsymbol{v}_k\}$、$\boldsymbol{u}_m = \boldsymbol{v}_m$ と取れば $\sum_m \boldsymbol{v}_m^{\mathsf{T}}S\boldsymbol{v}_m = \sum_{m\le k}\lambda_m$ で等号が成立します。ステップ 2 に戻して

$$\min_V E(V) = \operatorname{tr}S - \sum_{j=1}^{k}\lambda_j = \sum_{j=1}^{p}\lambda_j - \sum_{j=1}^{k}\lambda_j = \sum_{j=k+1}^{p}\lambda_j$$

を得ます（命題 4.2 の 3. を使いました）。

∎

注意 5.2

この定理と線形回帰の違いに注意してください。回帰は目的変数の方向（縦方向）に測った残差の二乗和を最小にします（定義 3.1［線形回帰と最小二乗法］）。PCA は直線・平面への垂直距離を最小にします。したがって $y$ を $x$ で回帰した直線と、$x$ を $y$ で回帰した直線と、PCA の第 1 主成分軸は、一般に三本とも別の直線です。どれが正しいかは「何を誤差とみなすか」という問題設定で決まります。説明変数にも測定誤差がある場合に垂直距離を採る立場は、全最小二乗法（total least squares）と呼ばれます。

6. 寄与率：どこまで採ればよいか

定理 5.1 は、$k$ 次元まで採ったときに捨てる誤差が $\lambda_{k+1} + \cdots + \lambda_p$ ちょうどであることを教えてくれます。全分散が $\sum_j \lambda_j$（命題 4.2）ですから、比を取れば「何割を残したか」が固有値だけで読めます。

定義 6.1（寄与率と累積寄与率）↗

$\operatorname{tr}S = \sum_{j=1}^p \lambda_j > 0$ とする。第 $k$ 主成分の寄与率を

$$r_k = \frac{\lambda_k}{\lambda_1 + \cdots + \lambda_p}$$

と定め、上位 $k$ 本までの累積寄与率を

$$R_k = \sum_{m=1}^{k} r_m = \frac{\lambda_1 + \cdots + \lambda_k}{\lambda_1 + \cdots + \lambda_p}$$

と定める。

$0 \le r_k \le 1$ かつ $\sum_k r_k = 1$ は定義と $\lambda_k \ge 0$ から直ちに従います。意味は次の等式に集約されます。定理 5.1 の最小誤差を $E_k^{*}$ と書けば

$$\frac{E_k^{*}}{\operatorname{tr}S} = \frac{\lambda_{k+1}+\cdots+\lambda_p}{\lambda_1+\cdots+\lambda_p} = 1 - R_k .$$

つまり累積寄与率 $R_k$ は、$k$ 次元に落としたときに保たれる分散の割合であり、同時に $1 - R_k$ が相対的な平均二乗再構成誤差です。「累積寄与率 90% で打ち切る」という操作は、「もとの散らばりの 10% ぶんの二乗誤差を許容して圧縮する」と言い換えられます。

注意

「累積寄与率 80% 以上」「固有値の減り方が緩やかになる肘（スクリープロットの折れ目）で切る」「相関行列の PCA では $\lambda > 1$ の主成分だけ採る（カイザー基準）」といった目安がよく使われますが、これらは定理ではなく経験則です。下流のタスク（分類精度、復元品質など）が測れるなら、$k$ は交差検証で選ぶほうが確実だと思います。

例 6.2（5 点の完全な計算）

$\mathbb{R}^2$ の 5 点 $(1,2), (2,3), (3,5), (4,4), (5,6)$ を PCA にかけます。

重心と中心化. $\bar{\boldsymbol{x}} = \frac{1}{5}(1+2+3+4+5,\ 2+3+5+4+6) = (3, 4)$。よって中心化ベクトルは

$$\boldsymbol{y}_1 = (-2,-2),\quad \boldsymbol{y}_2 = (-1,-1),\quad \boldsymbol{y}_3 = (0,1),\quad \boldsymbol{y}_4 = (1,0),\quad \boldsymbol{y}_5 = (2,2).$$

検算として $\sum_i \boldsymbol{y}_i = (-2-1+0+1+2,\ -2-1+1+0+2) = (0,0)$。

共分散行列. $\sum_i y_{i1}^2 = 4+1+0+1+4 = 10$、$\sum_i y_{i2}^2 = 4+1+1+0+4 = 10$、$\sum_i y_{i1}y_{i2} = 4+1+0+0+4 = 9$ なので

$$S = \frac{1}{5}\begin{pmatrix} 10 & 9 \\ 9 & 10 \end{pmatrix} = \begin{pmatrix} 2 & 1.8 \\ 1.8 & 2 \end{pmatrix}.$$

固有値と固有ベクトル. $\det(S - \lambda I) = (2-\lambda)^2 - 1.8^2 = 0$ より $2 - \lambda = \pm 1.8$、すなわち $\lambda_1 = 3.8$、$\lambda_2 = 0.2$。$\lambda_1$ に対しては $(S - 3.8 I)\boldsymbol{u} = \begin{pmatrix} -1.8 & 1.8 \\ 1.8 & -1.8\end{pmatrix}\boldsymbol{u} = \boldsymbol{0}$ から $u_1 = u_2$、よって $\boldsymbol{v}_1 = \frac{1}{\sqrt{2}}(1,1)^{\mathsf{T}}$。同様に $\lambda_2$ に対しては $u_1 = -u_2$ から $\boldsymbol{v}_2 = \frac{1}{\sqrt{2}}(1,-1)^{\mathsf{T}}$。内積は $\frac{1}{2}(1 - 1) = 0$ で直交しており、定理 2.5 と整合します。

主成分得点. $z_{i1} = \frac{y_{i1}+y_{i2}}{\sqrt{2}}$、$z_{i2} = \frac{y_{i1}-y_{i2}}{\sqrt{2}}$ なので

$$(z_{i1})_{i=1}^{5} = \tfrac{1}{\sqrt2}(-4,\, -2,\, 1,\, 1,\, 4), \qquad (z_{i2})_{i=1}^{5} = \tfrac{1}{\sqrt2}(0,\, 0,\, -1,\, 1,\, 0).$$

分散を確かめます。$\frac{1}{5}\sum_i z_{i1}^2 = \frac{1}{5}\cdot\frac{16+4+1+1+16}{2} = \frac{38}{10} = 3.8 = \lambda_1$、$\frac{1}{5}\sum_i z_{i2}^2 = \frac{1}{5}\cdot\frac{0+0+1+1+0}{2} = \frac{2}{10} = 0.2 = \lambda_2$。共分散は $\frac{1}{5}\cdot\frac{1}{2}\big(0 + 0 + (1)(-1) + (1)(1) + 0\big) = 0$。命題 4.2 の 2. のとおりです。

寄与率. $\operatorname{tr}S = 4 = 3.8 + 0.2$ なので $r_1 = 3.8/4 = 0.95$、$r_2 = 0.05$。第 1 主成分だけで全分散の 95% を説明します。

再構成誤差. $k = 1$ に落とすと、$\boldsymbol{y}_i$ の残差は $z_{i2}\boldsymbol{v}_2$ で、その二乗ノルムは $z_{i2}^2 = 0, 0, 0.5, 0.5, 0$。平均は $\frac{1}{5}(0+0+0.5+0.5+0) = 0.2 = \lambda_2$ で、定理 5.1 と一致します。復元値は $\hat{\boldsymbol{x}}_i = \bar{\boldsymbol{x}} + z_{i1}\boldsymbol{v}_1$、たとえば $i=3$ なら $z_{31} = 1/\sqrt2$ より $\hat{\boldsymbol{x}}_3 = (3,4) + \frac{1}{\sqrt2}\cdot\frac{1}{\sqrt2}(1,1) = (3.5, 4.5)$ です。

5 点の完全な計算の例のデータ（中心化済み）と 2 本の主成分軸。太い軸が第 1 主成分方向、破線が第 2 主成分方向で、細い点線は第 1 主成分軸への垂線です。点 (0,1) と (1,0) は第 1 主成分軸上の同じ点に落ちます。この 2 点を区別しているのは第 2 主成分だけです。

例 6.3（データが直線上に乗るとき）

すべての点が 1 本の直線上にある場合、すなわち単位ベクトル $\boldsymbol{a}$ と実数 $t_1, \ldots, t_n$（$\sum_i t_i = 0$ としてよい）によって $\boldsymbol{y}_i = t_i \boldsymbol{a}$ と書ける場合を考えます。このとき

$$S = \frac{1}{n}\sum_i (t_i\boldsymbol{a})(t_i\boldsymbol{a})^{\mathsf{T}} = \left(\frac{1}{n}\sum_i t_i^2\right)\boldsymbol{a}\boldsymbol{a}^{\mathsf{T}} = s_t^2\, \boldsymbol{a}\boldsymbol{a}^{\mathsf{T}}$$

です。この行列は $\boldsymbol{a}$ を $s_t^2\boldsymbol{a}$ に写し（$\boldsymbol{a}^{\mathsf{T}}\boldsymbol{a} = 1$ より）、$\boldsymbol{a}$ に直交するベクトル $\boldsymbol{w}$ を $s_t^2 \boldsymbol{a}(\boldsymbol{a}^{\mathsf{T}}\boldsymbol{w}) = \boldsymbol{0}$ に写します。したがって固有値は $\lambda_1 = s_t^2$（固有ベクトル $\boldsymbol{a}$）と $0$（重複度 $p-1$）。寄与率は $r_1 = 1$、$k=1$ での再構成誤差は $\lambda_2 + \cdots + \lambda_p = 0$。PCA は、データが実際に乗っている直線をそのまま見つけ出します。逆に言えば、固有値 $0$ はデータに厳密な線形従属関係があることの検出器です。

注意 6.4

$S\boldsymbol{u} = \frac{1}{n}\sum_i \boldsymbol{y}_i \langle \boldsymbol{y}_i, \boldsymbol{u}\rangle$ は常に $\operatorname{span}\{\boldsymbol{y}_1,\ldots,\boldsymbol{y}_n\}$ に属するので $\operatorname{rank}S \le \dim\operatorname{span}\{\boldsymbol{y}_i\}$ です。さらに $\sum_i \boldsymbol{y}_i = \boldsymbol{0}$ という線形関係があるため、$n$ 本のベクトル $\boldsymbol{y}_i$ が張る空間の次元は高々 $n-1$。したがって

$$\operatorname{rank}S \le \min(p,\ n-1)$$

です。標本数が次元以下（$n \le p$）のデータでは必ず $\lambda_p = 0$ となり、共分散行列は正則になりません。遺伝子発現データのように $p \gg n$ の場面では、非零の固有値は高々 $n-1$ 個しかなく、PCA が返せる主成分の本数も $n-1$ 本までです。

7. 使う前に確認すべきこと

PCA の手順。固有値分解までは決定的な計算で、判断が入るのは k の決定だけです。

中心化を忘れない. 平均を引かずに $\frac{1}{n}\sum_i \boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}}$ を固有値分解すると何が起きるでしょうか。$\boldsymbol{x}_i = \boldsymbol{y}_i + \bar{\boldsymbol{x}}$ を代入して展開すると、交差項が $\sum_i \boldsymbol{y}_i = \boldsymbol{0}$ で消えるので

$$\frac{1}{n}\sum_i \boldsymbol{x}_i\boldsymbol{x}_i^{\mathsf{T}} = S + \bar{\boldsymbol{x}}\,\bar{\boldsymbol{x}}^{\mathsf{T}}$$

となります。余計に足された $\bar{\boldsymbol{x}}\bar{\boldsymbol{x}}^{\mathsf{T}}$ は $\bar{\boldsymbol{x}}$ 方向に固有値 $\|\bar{\boldsymbol{x}}\|^2$ をもつ階数 1 の行列です。原点から遠いデータでは $\|\bar{\boldsymbol{x}}\|^2$ が $\lambda_1$ より大きくなりがちで、そのとき「第 1 主成分」は散らばりの方向ではなく重心の方向を指してしまいます。定理 5.1 のステップ 1 が示すとおり、重心を通らせるのが当てはめの意味でも最適です。

スケールに依存する. 変数ごとに単位を取り替えると、共分散行列は $\tilde{S} = D S D$（$D$ は正の対角行列）に変わりますが、固有ベクトルは $D\boldsymbol{v}_j$ にはなりません。

例 7.1（単位を変えると主成分が変わる）

$S = \begin{pmatrix} 4 & 2 \\ 2 & 4\end{pmatrix}$ とします。固有値は $4 \pm 2$、すなわち $\lambda_1 = 6$（固有ベクトル $\frac{1}{\sqrt2}(1,1)$）と $\lambda_2 = 2$（$\frac{1}{\sqrt2}(1,-1)$）で、第 1 主成分軸は $45^\circ$ の方向です。ここで第 2 変数の単位を $\frac{1}{10}$ に細かくする（たとえば m を dm にする）と、$D = \operatorname{diag}(1,10)$ として

$$\tilde{S} = D S D = \begin{pmatrix} 4 & 20 \\ 20 & 400 \end{pmatrix}.$$

$\operatorname{tr}\tilde S = 404$、$\det \tilde S = 1600 - 400 = 1200$ なので、固有値は $\lambda = \frac{404 \pm \sqrt{404^2 - 4\cdot 1200}}{2} = \frac{404\pm\sqrt{158416}}{2}$、数値では $\lambda_1 \approx 401.01$、$\lambda_2 \approx 2.99$。第 1 主成分方向は $(4 - 401.01)u_1 + 20u_2 = 0$ から $u_2 \approx 19.85\,u_1$、すなわち第 2 座標軸から約 $2.9^\circ$ しか離れていません。もとの $45^\circ$ とはまったく別の方向です。寄与率も $0.75$ から $0.993$ へ跳ね上がります。

これは PCA の欠陥ではなく、分散という量が単位に依存する（次元をもつ）ことの当然の帰結です。身長 [cm] と年収 [円] のように単位が比較不能な変数が混ざるときは、各変数を標準偏差で割ってから PCA にかけます。これは相関行列の固有値問題を解くことと同じで、相関行列 PCA と呼ばれます。

教師なしであることを忘れない. PCA はラベルを見ません。分散が大きい方向が、知りたい区別に対応しているとは限らないのです。

例 7.2（分散最大の方向が判別に役立たない例）

$\mathbb{R}^2$ に 2 つのクラスがあり、クラス A は $(0, 1)$ を中心に、クラス B は $(0,-1)$ を中心に、どちらも第 1 座標方向に標準偏差 $10$、第 2 座標方向に標準偏差 $0.1$ でばらついているとします。全データをまとめた共分散行列は、対称性から近似的に

$$S \approx \begin{pmatrix} 100 & 0 \\ 0 & 1 + 0.01 \end{pmatrix}$$

です（第 2 成分の分散は、クラス中心の差による $1$ とクラス内のばらつき $0.01$ の和）。よって $\lambda_1 = 100$、第 1 主成分方向は第 1 座標軸で、寄与率は $100/101.01 \approx 0.99$。ところが 2 つのクラスは第 1 座標では完全に重なっており、区別しているのは第 2 座標だけです。第 1 主成分だけを残すと、寄与率 99% を確保しながらクラスの情報を 100% 失います。

ラベルが使える場合は、クラス間分散とクラス内分散の比を最大化する線形判別分析（LDA）のような教師あり手法が適切です。PCA を前処理に使うなら、「分散の大きさ」と「タスクにとっての有用性」は別物だと意識してください。

線形であることを忘れない. PCA が探すのは部分空間、すなわち平らな集合です。データが曲がった曲面（渦巻き状の 2 次元曲面など）に沿って分布している場合、それを平面で近似することになり、うまく展開できません。非線形に拡張する道具としてはカーネル PCA や、ニューラルネットワークを使ったオートエンコーダがあります。逆に、隠れ層の活性化が恒等写像で損失が二乗誤差の場合、オートエンコーダの大域最適解が張る部分空間は PCA のそれと一致することが知られています（Baldi–Hornik, 1989）。線形の世界では PCA が最適解であり、非線形性を入れて初めてその先へ進めるわけです。

確率モデルとして見ることもできる. ここまでは標本共分散行列という記述統計の量だけを扱いましたが、$\boldsymbol{x}$ が期待値 $\boldsymbol{\mu}$、共分散行列 $\Sigma$ の確率ベクトルであるとして、$\operatorname{Var}(\boldsymbol{u}^{\mathsf{T}}\boldsymbol{x}) = \boldsymbol{u}^{\mathsf{T}}\Sigma\boldsymbol{u}$ を最大化する問題としても同じ理論が成り立ちます（確率変数と期待値、定義 5.1［確率変数と期待値］）。この視点を進めると、潜在変数モデルとしての確率的 PCA が得られます（確率論とベイズ統計の役割）。

8. 演習

演習 8.1標準

2 変数データの標本共分散行列が $S = \begin{pmatrix} 5 & 2 \\ 2 & 2\end{pmatrix}$ であるとする。(1) 固有値を求めよ。(2) 正規化した主成分方向を求め、直交していることを確かめよ。(3) 第 1 主成分の寄与率を求めよ。(4) 第 1 主成分だけを残したときの平均二乗再構成誤差を求めよ。

解答

(1) 特性多項式は

$$\det(S - \lambda I) = (5-\lambda)(2-\lambda) - 2\cdot 2 = \lambda^2 - 7\lambda + 6 = (\lambda - 6)(\lambda - 1)$$

なので $\lambda_1 = 6$、$\lambda_2 = 1$。検算：$\lambda_1 + \lambda_2 = 7 = \operatorname{tr}S$、$\lambda_1\lambda_2 = 6 = \det S$。

(2) $\lambda_1 = 6$ について $S - 6I = \begin{pmatrix} -1 & 2 \\ 2 & -4 \end{pmatrix}$。第 1 行から $-u_1 + 2u_2 = 0$、すなわち $u_1 = 2u_2$（第 2 行 $2u_1 - 4u_2 = 0$ は同じ条件です）。長さを $1$ にして $\boldsymbol{v}_1 = \frac{1}{\sqrt5}(2,1)^{\mathsf{T}}$。

$\lambda_2 = 1$ について $S - I = \begin{pmatrix} 4 & 2 \\ 2 & 1\end{pmatrix}$。第 1 行から $4u_1 + 2u_2 = 0$、すなわち $u_2 = -2u_1$。よって $\boldsymbol{v}_2 = \frac{1}{\sqrt5}(1,-2)^{\mathsf{T}}$。

内積は $\frac{1}{5}(2\cdot 1 + 1\cdot(-2)) = 0$ で直交しており、定理 2.5 と整合します。

(3) $\operatorname{tr}S = 7$ なので $r_1 = 6/7 \approx 0.857$。

(4) 定理 5.1 より、捨てられる誤差は残りの固有値の和 $\lambda_2 = 1$。全分散 $7$ の $1/7 \approx 14.3\%$ で、$1 - r_1$ に一致します。

演習 8.2標準

2 変数を標準化した後の相関行列 $R = \begin{pmatrix} 1 & \rho \\ \rho & 1 \end{pmatrix}$（$-1 < \rho < 1$）を考える。固有値と正規化した固有ベクトルを $\rho$ の式で求め、第 1 主成分の寄与率を書け。さらに $\rho \to 0$ と $\rho \to 1$ の極限で何が起きるかを述べよ。

解答

$\det(R - \lambda I) = (1-\lambda)^2 - \rho^2 = 0$ から $1 - \lambda = \pm\rho$、よって固有値は $1+\rho$ と $1-\rho$。

$\rho \ne 0$ とすると、固有値 $1+\rho$ については $R - (1+\rho)I = \begin{pmatrix} -\rho & \rho \\ \rho & -\rho\end{pmatrix}$ より $u_1 = u_2$、固有ベクトルは $\frac{1}{\sqrt2}(1,1)^{\mathsf{T}}$。固有値 $1-\rho$ については同様に $u_1 = -u_2$ で $\frac{1}{\sqrt2}(1,-1)^{\mathsf{T}}$。

大小関係は $\rho$ の符号で決まります。$\rho > 0$ なら $\lambda_1 = 1+\rho$（方向 $\frac{1}{\sqrt2}(1,1)$）、$\rho<0$ なら $\lambda_1 = 1-\rho$（方向 $\frac{1}{\sqrt2}(1,-1)$）。いずれの場合も $\lambda_1 = 1 + |\rho|$、$\operatorname{tr}R = 2$ なので

$$r_1 = \frac{1 + |\rho|}{2}.$$

$\rho \to 0$：$r_1 \to 1/2$。実際 $\rho = 0$ では $R = I_2$ となり 例 3.6 の状況で、主成分方向が定まりません。2 変数が無相関なら、1 次元に落とすと必ず分散の半分を失います。次元削減の利得がない、というのが正しい結論です。

$\rho \to 1$：$r_1 \to 1$。2 変数がほぼ同じ情報をもち、$\frac{1}{\sqrt2}(1,1)$ 方向の 1 次元でほぼ完全に表せます。捨てる誤差は $\lambda_2 = 1 - \rho \to 0$。

なお、相関行列は各変数の分散を $1$ に揃えた後の共分散行列なので、$\operatorname{tr}R = p$（ここでは $2$）が常に成り立ちます。カイザー基準「$\lambda > 1$ の主成分だけ採る」は、「平均的な 1 変数分（$\operatorname{tr}R/p = 1$）より多く説明する主成分を残す」という意味です。

演習 8.3難

$Q \in \mathbb{R}^{p\times p}$ を直交行列（$Q^{\mathsf{T}}Q = I$）、$\boldsymbol{b} \in \mathbb{R}^p$ を定ベクトルとし、$\tilde{\boldsymbol{x}}_i = Q\boldsymbol{x}_i + \boldsymbol{b}$ とおく。(1) 新しいデータの標本共分散行列が $\tilde{S} = QSQ^{\mathsf{T}}$ となることを示せ。(2) $\tilde S$ の固有値は $S$ の固有値と一致し、対応する正規直交固有ベクトルが $Q\boldsymbol{v}_j$ で与えられることを示せ。(3) PCA の寄与率と再構成誤差が回転・鏡映・平行移動で不変であることを結論せよ。(4) 一方 $D = \operatorname{diag}(1,10)$、$\tilde{\boldsymbol{x}}_i = D\boldsymbol{x}_i$ の場合には、$S = \begin{pmatrix}4&2\\2&4\end{pmatrix}$ に対して $D\boldsymbol{v}_1$ が $\tilde S = DSD$ の固有ベクトルにならないことを確かめよ。

解答

(1) 新しい重心は平均の線形性から $\frac{1}{n}\sum_i (Q\boldsymbol{x}_i + \boldsymbol{b}) = Q\bar{\boldsymbol{x}} + \boldsymbol{b}$。よって中心化ベクトルは

$$\tilde{\boldsymbol{y}}_i = (Q\boldsymbol{x}_i + \boldsymbol{b}) - (Q\bar{\boldsymbol{x}} + \boldsymbol{b}) = Q(\boldsymbol{x}_i - \bar{\boldsymbol{x}}) = Q\boldsymbol{y}_i$$

で、平行移動 $\boldsymbol{b}$ は消えます。したがって

$$\tilde S = \frac{1}{n}\sum_i (Q\boldsymbol{y}_i)(Q\boldsymbol{y}_i)^{\mathsf{T}} = Q\left(\frac{1}{n}\sum_i \boldsymbol{y}_i\boldsymbol{y}_i^{\mathsf{T}}\right)Q^{\mathsf{T}} = QSQ^{\mathsf{T}}.$$

(2) $S\boldsymbol{v}_j = \lambda_j\boldsymbol{v}_j$ の両辺に左から $Q$ を掛け、途中に $I = Q^{\mathsf{T}}Q$ を挿入します。

$$\tilde S (Q\boldsymbol{v}_j) = QSQ^{\mathsf{T}}Q\boldsymbol{v}_j = QS\boldsymbol{v}_j = \lambda_j (Q\boldsymbol{v}_j).$$

また $\langle Q\boldsymbol{v}_j, Q\boldsymbol{v}_l\rangle = \boldsymbol{v}_j^{\mathsf{T}}Q^{\mathsf{T}}Q\boldsymbol{v}_l = \boldsymbol{v}_j^{\mathsf{T}}\boldsymbol{v}_l = \delta_{jl}$ なので $\{Q\boldsymbol{v}_j\}$ は正規直交系、とくに $p$ 個の一次独立なベクトルです。よってこれが $\tilde S$ の正規直交固有基底を与え、固有値の組（重複度込み）は $S$ のそれと一致します。

(3) 寄与率は固有値の比だけで決まるので不変です。再構成誤差も 定理 5.1 より $\lambda_{k+1}+\cdots+\lambda_p$ で、これも不変。主成分方向だけがデータと一緒に $Q$ で回ります。幾何的には、PCA が座標系ではなく点の配置そのものを見ていることを意味します。

(4) $S = \begin{pmatrix}4&2\\2&4\end{pmatrix}$ の固有ベクトルは $\boldsymbol{v}_1 = \frac{1}{\sqrt2}(1,1)^{\mathsf{T}}$（$\lambda_1 = 6$）です。$D\boldsymbol{v}_1$ は $(1,10)^{\mathsf{T}}$ の定数倍。一方 $\tilde S = DSD = \begin{pmatrix} 4 & 20 \\ 20 & 400\end{pmatrix}$ なので

$$\tilde S \begin{pmatrix} 1 \\ 10\end{pmatrix} = \begin{pmatrix} 4 + 200 \\ 20 + 4000 \end{pmatrix} = \begin{pmatrix} 204 \\ 4020 \end{pmatrix}.$$

もし $(1,10)^{\mathsf{T}}$ が固有ベクトルなら、像も第 2 成分が第 1 成分の $10$ 倍、すなわち $(204, 2040)^{\mathsf{T}}$ でなければなりません。$4020 \ne 2040$ なので固有ベクトルではありません。例 7.1 で見たとおり、実際の第 1 主成分方向は $u_2/u_1 \approx 19.85$ の方向です。

結論として、PCA は直交変換（等長変換）に対しては同変ですが、一般の可逆線形変換に対してはそうではありません。だからこそ変数のスケーリングが本質的な前処理になります。

参考文献

• K. Pearson, “On Lines and Planes of Closest Fit to Systems of Points in Space”, Philosophical Magazine 2 (1901), 559–572. PCA を「点集合への最良当てはめ直線・平面」として導入した原論文。
• H. Hotelling, “Analysis of a Complex of Statistical Variables into Principal Components”, Journal of Educational Psychology 24 (1933), 417–441, 498–520. 分散最大化による定式化と principal component の命名。
• I. T. Jolliffe, Principal Component Analysis, 2nd ed., Springer, 2002. PCA の標準的なモノグラフ。第 1〜3 章に本記事の内容と、寄与率・主成分の解釈に関する詳細な議論があります。
• C. M. Bishop, Pattern Recognition and Machine Learning, Springer, 2006 — 第 12 章「Continuous Latent Variables」。分散最大化と誤差最小化の両方の定式化、および確率的 PCA。
• G. H. Golub and C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013 — 特異値分解とその数値計算を扱う章。
• P. Baldi and K. Hornik, “Neural networks and principal component analysis: Learning from examples without local minima”, Neural Networks 2 (1989), 53–58. 線形オートエンコーダの大域最適解が PCA の部分空間と一致することを示した論文。

Appendix: 特異値分解で計算する

共分散行列を作らずに済ませる. 実務の PCA は、$S$ を作ってから固有値分解するのではなく、中心化データ行列 $X$ の特異値分解（SVD）を直接計算するのが標準です。理由は次の命題と、その後の数値的な注意にあります。

命題 8.4（SVD と PCA の対応）✓

$X \in \mathbb{R}^{n\times p}$ を 定義 2.1 の中心化データ行列、$r = \operatorname{rank}X$ とし、その簡約特異値分解を

$$X = U \Sigma V^{\mathsf{T}}, \qquad U \in \mathbb{R}^{n\times r},\ \ \Sigma = \operatorname{diag}(\sigma_1 \ge \cdots \ge \sigma_r > 0),\ \ V \in \mathbb{R}^{p \times r}$$

（$U^{\mathsf{T}}U = V^{\mathsf{T}}V = I_r$）とする。このとき $V$ の第 $j$ 列は $S = \frac{1}{n}X^{\mathsf{T}}X$ の固有ベクトルであり、対応する固有値は $\lambda_j = \sigma_j^2/n$（$1 \le j \le r$）である。残る $p - r$ 個の固有値は $0$ である。さらに第 $j$ 主成分得点は $z_{ij} = \sigma_j u_{ij}$ で与えられる。

証明（命題 8.4）

$U^{\mathsf{T}}U = I_r$ を使って

$$S = \frac{1}{n}X^{\mathsf{T}}X = \frac{1}{n}V\Sigma U^{\mathsf{T}} U \Sigma V^{\mathsf{T}} = \frac{1}{n} V\Sigma^2 V^{\mathsf{T}} = \sum_{j=1}^{r}\frac{\sigma_j^2}{n}\, \boldsymbol{v}_j\boldsymbol{v}_j^{\mathsf{T}}$$

（$\boldsymbol{v}_j$ は $V$ の第 $j$ 列）。両辺に $\boldsymbol{v}_l$ を掛けると、$V$ の列の正規直交性から $S\boldsymbol{v}_l = \frac{\sigma_l^2}{n}\boldsymbol{v}_l$ を得ます。一方 $\boldsymbol{w} \in \mathbb{R}^p$ が $V$ の全列と直交すれば $S\boldsymbol{w} = \boldsymbol{0}$ なので、$V$ の列直交補空間（次元 $p-r$）はすべて固有値 $0$ の固有空間に含まれます。以上で $p$ 個の固有値が出そろいました。

得点については $z_{ij} = \langle \boldsymbol{v}_j, \boldsymbol{y}_i\rangle$ が行列 $XV$ の $(i,j)$ 成分であり、$XV = U\Sigma V^{\mathsf{T}}V = U\Sigma$ なので $z_{ij} = \sigma_j u_{ij}$ です。

∎

なぜ SVD のほうがよいのか. $S = \frac{1}{n}X^{\mathsf{T}}X$ を明示的に作ると、$X$ の特異値が二乗されるため条件数も二乗されます。$\sigma_r/\sigma_1$ が計算機の丸め誤差の水準に近いとき、$X^{\mathsf{T}}X$ の段階で小さい固有値の情報が失われます。$X$ に直接 SVD をかければこの二乗を経由しません。計算量の面でも、$p \gg n$ なら $p\times p$ 行列を作らずに済みます。

NumPy での確認. 例 6.2 の数値を両方の方法で再現します。

import numpy as np

X_raw = np.array([[1., 2.], [2., 3.], [3., 5.], [4., 4.], [5., 6.]])
X = X_raw - X_raw.mean(axis=0)      # 中心化
n = X.shape[0]

# 方法 1: 共分散行列の固有値分解（eigh は実対称行列用、固有値は昇順）
S = X.T @ X / n
lam, V = np.linalg.eigh(S)
lam, V = lam[::-1], V[:, ::-1]      # 降順に並べ替える

# 方法 2: 中心化データ行列の特異値分解
U, sigma, Vt = np.linalg.svd(X, full_matrices=False)
lam_svd = sigma ** 2 / n

print(lam)               # [3.8 0.2]
print(lam_svd)           # [3.8 0.2]
print(lam / lam.sum())   # 寄与率 [0.95 0.05]
print(V[:, 0])           # 第 1 主成分方向（符号は不定）
print(X @ V[:, :1])      # 第 1 主成分得点

eigh が固有値を昇順で返すので並べ替えが要ります。出力される固有ベクトルの符号は環境によって反転しうる点は 注意 3.5 のとおりで、得点の符号がまとめて変わるだけなので寄与率も再構成誤差も影響を受けません。