Skip to content

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

Prerequisite:ニューラルネットワークと逆伝播:連鎖律を計算グラフの上で逆向きに走らせる

Raw

This content is not available in your language yet.

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

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

Section titled “1. 動機:変数が多いことの何が困るのか”

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

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

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

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

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

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

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

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

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

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

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

Section titled “2. 準備:中心化と標本共分散行列”

Definition 2.1中心化データ行列と標本共分散行列

x1,,xnRp\boldsymbol{x}_1, \ldots, \boldsymbol{x}_n \in \mathbb{R}^pn2n \ge 2)を観測データとする。標本平均(重心)を

xˉ=1ni=1nxi\bar{\boldsymbol{x}} = \frac{1}{n}\sum_{i=1}^{n} \boldsymbol{x}_i

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

S=1nXTX=1ni=1nyiyiTRp×pS = \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)(j, l) 成分は sjl=1ni=1nyijyils_{jl} = \frac{1}{n}\sum_{i=1}^n y_{ij} y_{il} であり、sjjs_{jj} は第 jj 変数の標本分散、sjls_{jl} は第 jj 変数と第 ll 変数の標本共分散である。

定義の中で 1nXTX\frac{1}{n} X^{\mathsf{T}} X1niyiyiT\frac{1}{n}\sum_i \boldsymbol{y}_i \boldsymbol{y}_i^{\mathsf{T}} を同一視しました。これは行列積の定義から出ます。XTX^{\mathsf{T}} の第 ii 列が yi\boldsymbol{y}_i なので、(XTX)jl=i=1nXijXil=i=1nyijyil=(iyiyiT)jl(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} となります。

Remark 2.2

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

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

Definition 2.3方向への射影と射影分散

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

zi(u)=u,yi=uT(xixˉ)z_i(\boldsymbol{u}) = \langle \boldsymbol{u}, \boldsymbol{y}_i \rangle = \boldsymbol{u}^{\mathsf{T}}(\boldsymbol{x}_i - \bar{\boldsymbol{x}})

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

得点の平均は 1nizi(u)=uT(1niyi)=uT0=0\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 です(中心化しているので iyi=0\sum_i \boldsymbol{y}_i = \boldsymbol{0})。したがって射影分散は二乗平均そのもので、

1ni=1nzi(u)2=1ni=1nuTyiyiTu=uT(1ni=1nyiyiT)u=uTSu\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}

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

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

Lemma 2.4標本共分散行列は対称かつ半正定値

Definition 2.1SSST=SS^{\mathsf{T}} = S を満たし、さらに任意の uRp\boldsymbol{u} \in \mathbb{R}^p に対して uTSu0\boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} \ge 0 である。

Proof(Lemma 2.4)

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

ST=(1nXTX)T=1nXT(XT)T=1nXTX=S.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 .

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

uTSu=1nuTXTXu=1n(Xu)T(Xu)=1nXu20.\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 .

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

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

Theorem 2.5実対称行列のスペクトル定理

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

Proof(Theorem 2.5)

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

λj=λjvj2=λjvjTvj=vjT(λjvj)=vjTSvj0\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

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

Remark 2.6

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

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

Section titled “3. 第 1 主成分:分散最大化がレイリー商になる”

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

最大化uTSu(条件 u=1 のもとで)\text{最大化} \quad \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} \quad \text{(条件 } \|\boldsymbol{u}\| = 1 \text{ のもとで)}

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

Theorem 3.1レイリー商の最大値

SRp×pS \in \mathbb{R}^{p\times p} を実対称行列とし、Theorem 2.5 の記号を使う。このとき

maxu=1uTSu=λ1\max_{\|\boldsymbol{u}\| = 1} \boldsymbol{u}^{\mathsf{T}} S \boldsymbol{u} = \lambda_1

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

Proof(Theorem 3.1)

ステップ 1(成分表示). Theorem 2.5 より v1,,vp\boldsymbol{v}_1, \ldots, \boldsymbol{v}_pRp\mathbb{R}^p の正規直交基底なので、任意の u\boldsymbol{u}u=j=1pcjvj\boldsymbol{u} = \sum_{j=1}^p c_j \boldsymbol{v}_jcj=vj,uc_j = \langle \boldsymbol{v}_j, \boldsymbol{u}\rangle)と一意に書けます。正規直交性から

u2=jcjvj,lclvl=j,lcjclvj,vl=j=1pcj2\|\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

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

ステップ 2(二次形式の計算). Su=jcjSvj=jcjλjvjS\boldsymbol{u} = \sum_j c_j S\boldsymbol{v}_j = \sum_j c_j \lambda_j \boldsymbol{v}_j なので、同じく正規直交性から

uTSu=lclvl,jcjλjvj=j=1pλjcj2.\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 .

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

ステップ 3(上からの評価). すべての jjλjλ1\lambda_j \le \lambda_1 かつ cj20c_j^2 \ge 0 なので λjcj2λ1cj2\lambda_j c_j^2 \le \lambda_1 c_j^2、したがって

uTSu=jλjcj2λ1jcj2=λ1.\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(達成). u=v1\boldsymbol{u} = \boldsymbol{v}_1 は単位ベクトルであり、v1TSv1=v1T(λ1v1)=λ1\boldsymbol{v}_1^{\mathsf{T}} S \boldsymbol{v}_1 = \boldsymbol{v}_1^{\mathsf{T}}(\lambda_1 \boldsymbol{v}_1) = \lambda_1。よって上限 λ1\lambda_1 は実際に達成され、最大値です。

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

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

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

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

Definition 3.2主成分方向と主成分得点

Theorem 2.5 の正規直交固有系 v1,,vp\boldsymbol{v}_1, \ldots, \boldsymbol{v}_p(固有値の降順)に対し、vk\boldsymbol{v}_kkk 主成分方向と呼ぶ。第 ii 個体のkk 主成分得点

zik=vk,yi=vkT(xixˉ)z_{ik} = \langle \boldsymbol{v}_k, \boldsymbol{y}_i \rangle = \boldsymbol{v}_k^{\mathsf{T}}(\boldsymbol{x}_i - \bar{\boldsymbol{x}})

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

Corollary 3.3第 1 主成分は最大固有値の固有ベクトル

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

Proof(Corollary 3.3)

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

Remark 3.4

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

ukj,lSjlujul=lSklul+jSjkuj=2(Su)k\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

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

Remark 3.5

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

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

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

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

Section titled “4. 第 k 主成分と主成分得点の性質”

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

Theorem 4.1逐次的な分散最大化と固有値の対応

Theorem 2.5 の記号のもとで、1kp1 \le k \le p とする。集合

Wk={uRp:u=1, u,v1==u,vk1=0}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=1k = 1 のときは単位球面全体)の上で、uTSu\boldsymbol{u}^{\mathsf{T}}S\boldsymbol{u} の最大値は λk\lambda_k であり、u=vk\boldsymbol{u} = \boldsymbol{v}_k で達成される。

Proof(Theorem 4.1)

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

uTSu=j=kpλjcj2λkj=kpcj2=λk\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

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

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

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

Proposition 4.2主成分得点の性質

Definition 3.2 の主成分得点 zikz_{ik} について、次が成り立つ。

  1. kk について 1ni=1nzik=0\frac{1}{n}\sum_{i=1}^{n} z_{ik} = 0
  2. k,lk, l について 1ni=1nzikzil=λkδkl\frac{1}{n}\sum_{i=1}^{n} z_{ik} z_{il} = \lambda_k \delta_{kl}。すなわち第 kk 主成分得点の標本分散は λk\lambda_k であり、異なる主成分得点どうしの標本共分散は 00 である。
  3. k=1pλk=trS=j=1psjj=1ni=1nyi2\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
Proof(Proposition 4.2)

1. 内積の線形性から 1nizik=vk,1niyi=vk,0=0\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。ここで Definition 2.1 の中心化から iyi=0\sum_i \boldsymbol{y}_i = \boldsymbol{0} を使いました。

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

1nizikzil=vkT(1niyiyiT)vl=vkTSvl=λlvkTvl=λlδkl,\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=lk = l なら λk\lambda_kklk \ne l なら 00。1. より得点の平均が 00 なので、この量はそのまま標本共分散です。

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

trS=tr(1niyiyiT)=1niyi2.\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 .

また trS=jsjj\operatorname{tr}S = \sum_j s_{jj} はトレースの定義そのものです。最後に Theorem 2.5S=VΛVTS = V\Lambda V^{\mathsf{T}}tr(AB)=tr(BA)\operatorname{tr}(AB) = \operatorname{tr}(BA)VTV=IV^{\mathsf{T}}V = I から

trS=tr(VΛVT)=tr(ΛVTV)=trΛ=kλk.\operatorname{tr}S = \operatorname{tr}(V\Lambda V^{\mathsf{T}}) = \operatorname{tr}(\Lambda V^{\mathsf{T}}V) = \operatorname{tr}\Lambda = \sum_k \lambda_k .
  1. は PCA のもう一つの顔です。もとの変数は互いに相関していますが、主成分得点に移ると相関が消えますVTV^{\mathsf{T}} による座標変換は、共分散行列を対角化する変換にほかなりません(対角化とジョルダン標準形Theorem 3.2[対角化とジョルダン標準形])。3. は「全分散はどの座標系で測っても同じで、それが固有値の総和に等しい」という保存則で、次節の寄与率の分母になります。

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

Section titled “5. もう一つの顔:当てはめ誤差の最小化”

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

Theorem 5.1最良近似アフィン部分空間

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

E(c,V)=1ni=1ndist(xi,A)2=1ni=1n(IPV)(xic)2E(\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

と定める。このとき

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

であり、最小値は c=xˉ\boldsymbol{c} = \bar{\boldsymbol{x}}V=span{v1,,vk}V = \operatorname{span}\{\boldsymbol{v}_1, \ldots, \boldsymbol{v}_k\} で達成される(k=pk = p のとき右辺は 00 と読む)。

Proof(Theorem 5.1)

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

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

nE(c,V)=i(IPV)yi2+2i(IPV)yi, d+nd2.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 .

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

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

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

ii について平均を取り Proposition 4.2 の 3. を使うと

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

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

ステップ 3(最大化). u1,,uk\boldsymbol{u}_1, \ldots, \boldsymbol{u}_kVV の正規直交基底とすると PVy=m=1kum,yumP_V\boldsymbol{y} = \sum_{m=1}^k \langle \boldsymbol{u}_m, \boldsymbol{y}\rangle \boldsymbol{u}_m、したがって PVy2=m=1kum,y2\|P_V\boldsymbol{y}\|^2 = \sum_{m=1}^k \langle \boldsymbol{u}_m,\boldsymbol{y}\rangle^2 です。§2 の射影分散の計算を各 um\boldsymbol{u}_m に適用すると

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

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

  • 0wj10 \le w_j \le 1wj=PVvj2w_j = \|P_V \boldsymbol{v}_j\|^2 であり、PVvjvj=1\|P_V\boldsymbol{v}_j\| \le \|\boldsymbol{v}_j\| = 1ベッセルの不等式(Remark 7.3)[内積空間とグラム・シュミット直交化]、あるいはステップ 2 のピタゴラス)。
  • j=1pwj=k\sum_{j=1}^p w_j = k:和の順序を変えると jwj=m=1kj=1pvj,um2=m=1kum2=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 = kvj\boldsymbol{v}_j が正規直交基底であることによるパーセバルの等式)。

このとき

j=1pλjwjj=1kλj=jkλj(wj1)+j>kλjwjjkλk(wj1)+j>kλkwj=λk(j=1pwjk)=0.\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 .

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

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

minVE(V)=trSj=1kλj=j=1pλjj=1kλj=j=k+1pλj\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

を得ます(Proposition 4.2 の 3. を使いました)。

Remark 5.2

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

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

Section titled “6. 寄与率:どこまで採ればよいか”

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

Definition 6.1寄与率と累積寄与率

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

rk=λkλ1++λpr_k = \frac{\lambda_k}{\lambda_1 + \cdots + \lambda_p}

と定め、上位 kk 本までの累積寄与率

Rk=m=1krm=λ1++λkλ1++λpR_k = \sum_{m=1}^{k} r_m = \frac{\lambda_1 + \cdots + \lambda_k}{\lambda_1 + \cdots + \lambda_p}

と定める。

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

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

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

Example 6.25 点の完全な計算

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

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

y1=(2,2),y2=(1,1),y3=(0,1),y4=(1,0),y5=(2,2).\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).

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

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

S=15(109910)=(21.81.82).S = \frac{1}{5}\begin{pmatrix} 10 & 9 \\ 9 & 10 \end{pmatrix} = \begin{pmatrix} 2 & 1.8 \\ 1.8 & 2 \end{pmatrix}.

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

主成分得点. zi1=yi1+yi22z_{i1} = \frac{y_{i1}+y_{i2}}{\sqrt{2}}zi2=yi1yi22z_{i2} = \frac{y_{i1}-y_{i2}}{\sqrt{2}} なので

(zi1)i=15=12(4,2,1,1,4),(zi2)i=15=12(0,0,1,1,0).(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).

分散を確かめます。15izi12=1516+4+1+1+162=3810=3.8=λ1\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_115izi22=150+0+1+1+02=210=0.2=λ2\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。共分散は 1512(0+0+(1)(1)+(1)(1)+0)=0\frac{1}{5}\cdot\frac{1}{2}\big(0 + 0 + (1)(-1) + (1)(1) + 0\big) = 0Proposition 4.2 の 2. のとおりです。

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

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

第 1 主成分(λ₁ = 3.8)第 2 主成分(λ₂ = 0.2)y₁y₂
ex-two-dim のデータ(中心化済み)と 2 本の主成分軸。太い軸が第 1 主成分方向、破線が第 2 主成分方向で、細い点線は第 1 主成分軸への垂線です。点 (0,1) と (1,0) は第 1 主成分軸上の同じ点に落ちます。この 2 点を区別しているのは第 2 主成分だけです。

Example 6.3データが直線上に乗るとき

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

S=1ni(tia)(tia)T=(1niti2)aaT=st2aaTS = \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}}

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

Remark 6.4

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

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

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

flowchart TD
A["生データ x_1, ..., x_n"] --> B["重心を引いて中心化"]
B --> C["標本共分散行列 S を作る"]
C --> D["固有値分解: λ_1 ≥ ... ≥ λ_p ≥ 0"]
D --> E["累積寄与率などで次元 k を決める"]
E --> F["上位 k 本の固有ベクトルに射影"]
F --> G["k 次元表現 z_i を得る"]
D --> H["捨てる誤差 = λ_(k+1) + ... + λ_p"]
PCA の手順。固有値分解までは決定的な計算で、判断が入るのは k の決定だけです。

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

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

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

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

Example 7.1単位を変えると主成分が変わる

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

S~=DSD=(42020400).\tilde{S} = D S D = \begin{pmatrix} 4 & 20 \\ 20 & 400 \end{pmatrix}.

trS~=404\operatorname{tr}\tilde S = 404detS~=1600400=1200\det \tilde S = 1600 - 400 = 1200 なので、固有値は λ=404±4042412002=404±1584162\lambda = \frac{404 \pm \sqrt{404^2 - 4\cdot 1200}}{2} = \frac{404\pm\sqrt{158416}}{2}、数値では λ1401.01\lambda_1 \approx 401.01λ22.99\lambda_2 \approx 2.99。第 1 主成分方向は (4401.01)u1+20u2=0(4 - 401.01)u_1 + 20u_2 = 0 から u219.85u1u_2 \approx 19.85\,u_1、すなわち第 2 座標軸から約 2.92.9^\circ しか離れていません。もとの 4545^\circ とはまったく別の方向です。寄与率も 0.750.75 から 0.9930.993 へ跳ね上がります。

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

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

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

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

S(100001+0.01)S \approx \begin{pmatrix} 100 & 0 \\ 0 & 1 + 0.01 \end{pmatrix}

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

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

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

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

Exercise 8.1標準

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

Solution

(1) 特性多項式は

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

なので λ1=6\lambda_1 = 6λ2=1\lambda_2 = 1。検算:λ1+λ2=7=trS\lambda_1 + \lambda_2 = 7 = \operatorname{tr}Sλ1λ2=6=detS\lambda_1\lambda_2 = 6 = \det S

(2) λ1=6\lambda_1 = 6 について S6I=(1224)S - 6I = \begin{pmatrix} -1 & 2 \\ 2 & -4 \end{pmatrix}。第 1 行から u1+2u2=0-u_1 + 2u_2 = 0、すなわち u1=2u2u_1 = 2u_2(第 2 行 2u14u2=02u_1 - 4u_2 = 0 は同じ条件です)。長さを 11 にして v1=15(2,1)T\boldsymbol{v}_1 = \frac{1}{\sqrt5}(2,1)^{\mathsf{T}}

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

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

(3) trS=7\operatorname{tr}S = 7 なので r1=6/70.857r_1 = 6/7 \approx 0.857

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

Exercise 8.2標準

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

Solution

det(RλI)=(1λ)2ρ2=0\det(R - \lambda I) = (1-\lambda)^2 - \rho^2 = 0 から 1λ=±ρ1 - \lambda = \pm\rho、よって固有値は 1+ρ1+\rho1ρ1-\rho

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

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

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

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

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

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

Exercise 8.3

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

Solution

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

y~i=(Qxi+b)(Qxˉ+b)=Q(xixˉ)=Qyi\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

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

S~=1ni(Qyi)(Qyi)T=Q(1niyiyiT)QT=QSQT.\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) Svj=λjvjS\boldsymbol{v}_j = \lambda_j\boldsymbol{v}_j の両辺に左から QQ を掛け、途中に I=QTQI = Q^{\mathsf{T}}Q を挿入します。

S~(Qvj)=QSQTQvj=QSvj=λj(Qvj).\tilde S (Q\boldsymbol{v}_j) = QSQ^{\mathsf{T}}Q\boldsymbol{v}_j = QS\boldsymbol{v}_j = \lambda_j (Q\boldsymbol{v}_j).

また Qvj,Qvl=vjTQTQvl=vjTvl=δjl\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} なので {Qvj}\{Q\boldsymbol{v}_j\} は正規直交系、とくに pp 個の一次独立なベクトルです。よってこれが S~\tilde S の正規直交固有基底を与え、固有値の組(重複度込み)は SS のそれと一致します。

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

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

S~(110)=(4+20020+4000)=(2044020).\tilde S \begin{pmatrix} 1 \\ 10\end{pmatrix} = \begin{pmatrix} 4 + 200 \\ 20 + 4000 \end{pmatrix} = \begin{pmatrix} 204 \\ 4020 \end{pmatrix}.

もし (1,10)T(1,10)^{\mathsf{T}} が固有ベクトルなら、像も第 2 成分が第 1 成分の 1010 倍、すなわち (204,2040)T(204, 2040)^{\mathsf{T}} でなければなりません。402020404020 \ne 2040 なので固有ベクトルではありません。Example 7.1 で見たとおり、実際の第 1 主成分方向は u2/u119.85u_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 の部分空間と一致することを示した論文。

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

Proposition 8.4SVD と PCA の対応

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

X=UΣVT,URn×r,  Σ=diag(σ1σr>0),  VRp×rX = 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}

UTU=VTV=IrU^{\mathsf{T}}U = V^{\mathsf{T}}V = I_r)とする。このとき VV の第 jj 列は S=1nXTXS = \frac{1}{n}X^{\mathsf{T}}X の固有ベクトルであり、対応する固有値は λj=σj2/n\lambda_j = \sigma_j^2/n1jr1 \le j \le r)である。残る prp - r 個の固有値は 00 である。さらに第 jj 主成分得点は zij=σjuijz_{ij} = \sigma_j u_{ij} で与えられる。

Proof(Proposition 8.4)

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

S=1nXTX=1nVΣUTUΣVT=1nVΣ2VT=j=1rσj2nvjvjTS = \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}}

vj\boldsymbol{v}_jVV の第 jj 列)。両辺に vl\boldsymbol{v}_l を掛けると、VV の列の正規直交性から Svl=σl2nvlS\boldsymbol{v}_l = \frac{\sigma_l^2}{n}\boldsymbol{v}_l を得ます。一方 wRp\boldsymbol{w} \in \mathbb{R}^pVV の全列と直交すれば Sw=0S\boldsymbol{w} = \boldsymbol{0} なので、VV の列直交補空間(次元 prp-r)はすべて固有値 00 の固有空間に含まれます。以上で pp 個の固有値が出そろいました。

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

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

NumPy での確認. Example 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 が固有値を昇順で返すので並べ替えが要ります。出力される固有ベクトルの符号は環境によって反転しうる点は Remark 3.5 のとおりで、得点の符号がまとめて変わるだけなので寄与率も再構成誤差も影響を受けません。

Report an error in this article ・Operated by: Mugen Giken LLCPricingTermsLegal notice

© 2026 夢現技研合同会社 ・Feeding the text to an LLM is welcome. Code samples are MIT licensed.