第3章主成分分析と次元削減（機械学習の数理）

この章の目標

主成分分析を分散の最大化として定式化し、主成分が標本共分散行列の固有ベクトルであることをレイリー商から証明できる

中心化が必要な理由を説明し、主成分分析を特異値分解で計算できる

エッカート–ヤングの定理を証明し、寄与率・白色化・推薦の低ランク近似（欠損値の注意を含む）を理解する

$k$ 平均法の目的関数が単調に減少して有限回で止まることを証明し、局所解に止まる例を挙げられる

ジョンソン–リンデンシュトラウスの補題を正確に述べ、主成分分析との違いを説明できる

前提：02-linear-algebra 第7章（直交射影・実対称行列の対角化）、02-linear-algebra 第8章（特異値分解・レイリー商）。3.8 節の証明では 22-statistics 第1章の正規分布・モーメント母関数・マルコフの不等式を使う。3.3・3.4 節では第1章のデータ漏洩と交差検証を、3.6 節では第2章のリッジ回帰を参照する。

通販サイトが顧客ごとに 200 の商品カテゴリ別の購入額を記録しているとする。顧客は $\mathbb{R}^{200}$ の点だが、眺めることもできず、似た顧客を探す計算も重い。一方で購入額どうしは強く相関しており（子育て中の世帯はおむつも粉ミルクも買う）、少数の要因で大部分が説明できることが多い。情報をなるべく失わずに低い次元で表すことを次元削減 (dimensionality reduction) という。

本章の中心は主成分分析 (principal component analysis, PCA) で、02-linear-algebra 第8章の線形代数がそのまま使える。同じ線形代数は推薦のための低ランク近似の土台でもあるが、欠損の多い実データでは特異値分解をそのまま使うと誤る。後半では、必ず止まるが最適とは限らない $k$ 平均法と、ランダムな線形写像で距離をほぼ保って次元を下げるジョンソン–リンデンシュトラウスの補題を扱う。

データ $x_1, \dots, x_n \in \mathbb{R}^d$ に対し、 $x_i^{\top}$ を第 $i$ 行とする $n \times d$ 行列 $X$ をデータ行列という。ベクトルは縦ベクトルで、転置を $A^{\top}$ と書く（02-linear-algebra の ${}^tA$ 、実行列の $A^{\ast}$ と同じ）。 $\mathbf{1}$ は成分がすべて $1$ のベクトル、 $\lVert A \rVert$ と $\lVert A \rVert_F$ は作用素ノルムとフロベニウスノルムである。

3.1 データ行列と中心化

定義 3.1（標本平均・中心化・標本共分散行列） $\bar{x} = \frac{1}{n}\sum_{i=1}^n x_i$ を標本平均、 $\tilde{x}_i = x_i - \bar{x}$ を中心化 (centering) したデータ、 $\tilde{x}_i^{\top}$ を第 $i$ 行とする $X_c = X - \mathbf{1}\bar{x}^{\top}$ を中心化データ行列といい、
$\hat{\Sigma} = \frac{1}{n}\sum_{i=1}^n \tilde{x}_i\tilde{x}_i^{\top} = \frac{1}{n}X_c^{\top}X_c$
を標本共分散行列 (sample covariance matrix) という。

本章では $n$ で割るが、 $n - 1$ で割る流儀もある（NumPy の np.cov の既定）。定数倍なので固有ベクトルや寄与率は変わらない。 $w^{\top}\hat{\Sigma}w = \frac{1}{n}\sum_i (w^{\top}\tilde{x}_i)^2 \geq 0$ なので $\hat{\Sigma}$ は半正定値の実対称行列であり、固有値 $\lambda_1 \geq \cdots \geq \lambda_d \geq 0$ と正規直交固有ベクトル $v_1, \dots, v_d$ をもつ（02-linear-algebra 第7章定理 7.32）。本章ではこの記号を通して使う。

補題 3.2 任意の $c, w \in \mathbb{R}^d$ について
$\frac{1}{n}\sum_{i=1}^n \lVert x_i - c \rVert^2 = \frac{1}{n}\sum_{i=1}^n \lVert \tilde{x}_i \rVert^2 + \lVert \bar{x} - c \rVert^2, \qquad \frac{1}{n}\sum_{i=1}^n (w^{\top}x_i - w^{\top}\bar{x})^2 = w^{\top}\hat{\Sigma}w$
である。特に $\sum_i \lVert x_i - c \rVert^2$ を最小にする $c$ は $\bar{x}$ だけで、 $\frac{1}{n}\sum_i \lVert \tilde{x}_i \rVert^2 = \operatorname{tr} \hat{\Sigma}$ （全分散）である。

証明. $\sum_i \tilde{x}_i = 0$ なので、 $x_i - c = \tilde{x}_i + (\bar{x} - c)$ を展開すると交差項が消えて第 1 式を得る。第 2 式は $(w^{\top}\tilde{x}_i)^2 = w^{\top}\tilde{x}_i\tilde{x}_i^{\top}w$ を平均すればよく、最後は $\operatorname{tr}(\tilde{x}_i\tilde{x}_i^{\top}) = \lVert \tilde{x}_i \rVert^2$ による。 $\square$

3.2 分散最大化としての主成分分析

データを平均を通る方向 $w$ （単位ベクトル）の直線に射影すると、各点は数 $w^{\top}\tilde{x}_i$ で表される。データの違いをなるべく残すため、その分散 $w^{\top}\hat{\Sigma}w$ （補題 3.2）が大きい方向を選ぶ。

定義 3.3（主成分） $\lVert w \rVert = 1$ のもとで $w^{\top}\hat{\Sigma}w$ を最大にする $w_1$ を第 1 主成分方向、 $w_1, \dots, w_{j-1}$ を定めたあと、さらに $w \perp w_1, \dots, w_{j-1}$ のもとで最大にする $w_j$ を第 $j$ 主成分方向という。 $z_{ij} = w_j^{\top}\tilde{x}_i$ を第 $j$ 主成分得点 (principal component score) という。

定理 3.4（主成分と固有ベクトル）各 $j$ について、 $\lVert w \rVert = 1$ かつ $w \perp v_1, \dots, v_{j-1}$ のもとでの $w^{\top}\hat{\Sigma}w$ の最大値は $\lambda_j$ で、 $w = v_j$ で達成される。したがって $v_1, \dots, v_d$ は主成分方向の列であり、第 $j$ 主成分得点の分散は $\lambda_j$ 、異なる主成分の得点の標本共分散は $0$ である。固有値が相異なれば、主成分方向は符号を除いて一意に $w_j = \pm v_j$ である。

証明. 条件のもとで $w = \sum_{i \geq j} c_iv_i$ , $\sum_i c_i^2 = 1$ と書け（02-linear-algebra 第7章命題 7.8）、 $w^{\top}\hat{\Sigma}w = \sum_{i \geq j} \lambda_ic_i^2 \leq \lambda_j$ で、等号は $w = v_j$ で成り立つ（02-linear-algebra 第8章 8.10 節のレイリー商の議論。 $j = 1$ が命題 8.29）。固有値が相異なれば、等号には $i > j$ で $c_i = 0$ が必要なので $w = \pm v_j$ で、帰納法で $w_j = \pm v_j$ 。得点の分散は $v_j^{\top}\hat{\Sigma}v_j = \lambda_j$ （補題 3.2）、共分散は $v_j^{\top}\hat{\Sigma}v_l = \lambda_lv_j^{\top}v_l = 0$ （ $j \neq l$ ）。 $\square$

$k$ 個の主成分を同時に考えると、分散の最大化と再構成誤差の最小化が一致する。

定理 3.5（主成分部分空間の最適性） $M$ を $d$ 次実対称行列とし、固有値を $\lambda_1 \geq \cdots \geq \lambda_d$ 、正規直交固有ベクトルを $v_1, \dots, v_d$ 、 $V_k = (v_1 \ \cdots \ v_k)$ とする。 $W^{\top}W = I_k$ を満たす任意の $d \times k$ 行列 $W$ について
$\operatorname{tr}(W^{\top}MW) \leq \lambda_1 + \cdots + \lambda_k = \operatorname{tr}(V_k^{\top}MV_k)$
である。 $M = \hat{\Sigma}$ （標本共分散行列）のとき、 $W$ の列空間への直交射影 $WW^{\top}$ について
$\frac{1}{n}\sum_{i=1}^n \lVert \tilde{x}_i - WW^{\top}\tilde{x}_i \rVert^2 = \operatorname{tr} \hat{\Sigma} - \operatorname{tr}(W^{\top}\hat{\Sigma}W) \geq \lambda_{k+1} + \cdots + \lambda_d$
である。つまり第 1〜第 $k$ 主成分方向が張る部分空間は、射影の分散の和を最大にし、同時に平均二乗再構成誤差を最小にする。

証明. $c_l = \lVert W^{\top}v_l \rVert^2$ とおく。 $V = (v_1 \ \cdots \ v_d)$ は直交行列なので $\sum_l c_l = \operatorname{tr}(W^{\top}VV^{\top}W) = \operatorname{tr}(W^{\top}W) = k$ 。 $WW^{\top}$ は直交射影（02-linear-algebra 第7章定理 7.15）で $\lVert WW^{\top}v \rVert^2 = v^{\top}WW^{\top}v = \lVert W^{\top}v \rVert^2$ だから、 $c_l = \lVert WW^{\top}v_l \rVert^2 \leq 1$ 。 $M = \sum_l \lambda_lv_lv_l^{\top}$ より $\operatorname{tr}(W^{\top}MW) = \sum_l \lambda_lc_l$ で、

\sum_{l=1}^d \lambda_lc_l - \sum_{l=1}^k \lambda_l = \sum_{l \leq k} \lambda_l(c_l - 1) + \sum_{l > k} \lambda_lc_l \leq \lambda_k\Bigl(\sum_{l \leq k} (c_l - 1) + \sum_{l > k} c_l\Bigr) = \lambda_k\Bigl(\sum_{l=1}^d c_l - k\Bigr) = 0

（ $l \leq k$ では $c_l - 1 \leq 0$ , $\lambda_l \geq \lambda_k$ 、 $l > k$ では $c_l \geq 0$ , $\lambda_l \leq \lambda_k$ ）。 $W = V_k$ では $c_l$ が $1$ （ $l \leq k$ ）か $0$ で等号が成り立つ。後半：ピタゴラスの定理より $\lVert \tilde{x}_i \rVert^2 = \lVert W^{\top}\tilde{x}_i \rVert^2 + \lVert \tilde{x}_i - WW^{\top}\tilde{x}_i \rVert^2$ で、 $W$ の列 $w_l$ について補題 3.2 より $\frac{1}{n}\sum_i \lVert W^{\top}\tilde{x}_i \rVert^2 = \sum_l w_l^{\top}\hat{\Sigma}w_l = \operatorname{tr}(W^{\top}\hat{\Sigma}W)$ 。 $\square$

例 3.6 5 点 $(1, 4), (2, 3), (3, 5), (4, 7), (5, 6)$ の平均は $\bar{x} = (3, 5)$ 、中心化した点は $(-2, -1), (-1, -2), (0, 0), (1, 2), (2, 1)$ で

\hat{\Sigma} = \begin{pmatrix} 2 & 8/5 \\ 8/5 & 2 \end{pmatrix}, \qquad \lambda_1 = \frac{18}{5},\ v_1 = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ 1 \end{pmatrix}, \qquad \lambda_2 = \frac{2}{5},\ v_2 = \frac{1}{\sqrt{2}}\begin{pmatrix} 1 \\ -1 \end{pmatrix}

である。第 1 主成分得点は $(-3, -3, 0, 3, 3)/\sqrt{2}$ （分散 $3.6$ ）、第 2 主成分得点は $(-1, 1, 0, -1, 1)/\sqrt{2}$ （分散 $0.4$ ）で、両者の共分散は $0$ 。第 1 主成分だけで再構成した点 $\bar{x} + z_{i1}v_1$ （例えば $(1, 4)$ は $(1.5, 3.5)$ に写る）の平均二乗再構成誤差は $\lambda_2 = 0.4$ である。

射影先が平均を通るのは偶然ではない。

命題 3.7（中心化の必要性） $P$ を部分空間 $L$ への直交射影、 $Q = I - P$ とする。データからアフィン部分空間 $m + L$ までの距離の 2 乗の平均は
$E(m) = \frac{1}{n}\sum_{i=1}^n \lVert Q(x_i - m) \rVert^2 = \frac{1}{n}\sum_{i=1}^n \lVert Q\tilde{x}_i \rVert^2 + \lVert Q(\bar{x} - m) \rVert^2$
である。よって $E$ は $m = \bar{x}$ で最小になり、 $k$ 次元アフィン部分空間とデータとの距離の 2 乗の平均の最小値は $\lambda_{k+1} + \cdots + \lambda_d$ で、 $\bar{x} + \operatorname{span}(v_1, \dots, v_k)$ で達成される。

証明. $x$ と $m + L$ の距離は $\lVert Q(x - m) \rVert$ である（02-linear-algebra 第7章定理 7.17）。点 $Qx_i$ の平均は $Q\bar{x}$ なので、補題 3.2 を点 $Qx_i$ と $c = Qm$ に適用すればよい。 $m = \bar{x}$ のあとの最小化は定理 3.5 である。 $\square$

例 3.8（中心化を忘れると） $\frac{1}{n}X^{\top}X = \hat{\Sigma} + \bar{x}\bar{x}^{\top}$ なので、中心化せずに固有ベクトルを求めると平均の方向が混ざる。例 3.6 では $\frac{1}{5}X^{\top}X$ （第 1 行 $(11, 83/5)$ 、第 2 行 $(83/5, 27)$ ）の第 1 固有ベクトルは第 1 軸から約 $57.9^\circ$ の方向で、 $v_1$ （ $45^\circ$ ）より平均 $(3, 5)$ の方向（約 $59.0^\circ$ ）に近い。データが原点から遠いほど、これは散らばりではなく位置を表す。

3.3 特異値分解との関係とエッカート–ヤングの定理

命題 3.9（主成分分析と特異値分解） $X_c$ の特異値分解を $X_c = U\Sigma V^{\top} = \sum_{j=1}^r \sigma_ju_jv_j^{\top}$ （ $r = \operatorname{rank} X_c$ 、02-linear-algebra 第8章定理 8.19。 $\Sigma$ は特異値を並べた $n \times d$ 行列で、 $\hat{\Sigma}$ とは別のもの）とする。

右特異ベクトル $v_j$ は $\hat{\Sigma}$ の固有値 $\lambda_j = \sigma_j^2/n$ の固有ベクトルである（ $j > r$ では $\lambda_j = 0$ ）。

第 $j$ 主成分得点を並べたベクトルは $X_cv_j = \sigma_ju_j$ である。

第 $k$ 主成分までで再構成した $X_cV_kV_k^{\top}$ は、特異値分解を $k$ 項で打ち切った $\sum_{j=1}^k \sigma_ju_jv_j^{\top}$ である。

証明. 1 は $X_c^{\top}X_c = V\Sigma^{\top}\Sigma V^{\top}$ から、2 は $X_cv_j = U\Sigma e_j$ から、3 は $v_j^{\top}V_kV_k^{\top}$ が $j \leq k$ なら $v_j^{\top}$ 、 $j > k$ なら $0$ であることから従う。 $\square$

例 3.6 では $\sigma_1 = 3\sqrt{2}$ , $\sigma_2 = \sqrt{2}$ で $\sigma_j^2/5 = \lambda_j$ である。命題 3.9 の 3 が最良の近似であることを次の定理が保証する。

定理 3.10（エッカート–ヤングの定理, Eckart–Young theorem） $A \in \operatorname{M}_{m,n}(\mathbb{R})$ の特異値分解を $A = \sum_{j=1}^r \sigma_ju_jv_j^{\top}$ （ $\sigma_1 \geq \cdots \geq \sigma_r > 0$ ）、 $k < r$ について $A_k = \sum_{j=1}^k \sigma_ju_jv_j^{\top}$ とする。 $\operatorname{rank} B \leq k$ を満たす任意の $B \in \operatorname{M}_{m,n}(\mathbb{R})$ について
$\lVert A - B \rVert_F^2 \geq \sum_{j=k+1}^r \sigma_j^2 = \lVert A - A_k \rVert_F^2, \qquad \lVert A - B \rVert \geq \sigma_{k+1} = \lVert A - A_k \rVert$
が成り立つ。すなわち $A_k$ は、フロベニウスノルムでも作用素ノルムでも、階数 $k$ 以下の行列による最良近似である。

証明. 等号部分： $A - A_k = \sum_{j > k} \sigma_ju_jv_j^{\top}$ は特異値 $\sigma_{k+1}, \dots, \sigma_r$ の特異値分解の形なので、02-linear-algebra 第8章命題 8.21 による。

フロベニウスノルムの不等式： $A, B$ の第 $i$ 行を $a_i^{\top}, b_i^{\top}$ とする。 $B$ の行空間を含む $k$ 次元部分空間の正規直交基底 $w_1, \dots, w_k$ をとり（02-linear-algebra 第7章系 7.10）、 $W = (w_1 \ \cdots \ w_k)$ とする。 $b_i$ は $W$ の列空間にあるので、最良近似定理（同定理 7.17）とピタゴラスの定理より

\lVert A - B \rVert_F^2 = \sum_i \lVert a_i - b_i \rVert^2 \geq \sum_i \lVert a_i - WW^{\top}a_i \rVert^2 = \sum_i \bigl(\lVert a_i \rVert^2 - \lVert W^{\top}a_i \rVert^2\bigr) = \lVert A \rVert_F^2 - \operatorname{tr}(W^{\top}A^{\top}AW)

定理 3.5 の前半を対称行列 $A^{\top}A$ （固有値は $\sigma_1^2 \geq \cdots \geq \sigma_r^2$ と $0$ ）に適用すると $\operatorname{tr}(W^{\top}A^{\top}AW) \leq \sigma_1^2 + \cdots + \sigma_k^2$ 。 $\lVert A \rVert_F^2 = \sum_j \sigma_j^2$ と合わせて結論を得る。 $\square$

作用素ノルムの不等式は主張にとどめる（証明は 02-linear-algebra 第8章定理 8.27。 $\operatorname{Ker} B$ と $\operatorname{span}(v_1, \dots, v_{k+1})$ が交わることを使う）。同定理のフロベニウスノルムの場合の証明（ワイルの不等式による）とここでの証明は別の道筋である。

最小点の一意性は特異値の重なりで決まる。フロベニウスノルムでは、 $\sigma_k > \sigma_{k+1}$ なら最小点は $A_k$ だけである。実際、上の証明で等号が成り立つには、定理 3.5 の証明の $l > k$ の項（ $\sigma_l^2 < \sigma_k^2$ ）から $c_l = \lVert W^{\top}v_l \rVert^2 = 0$ 、すなわち $W$ の列空間が $\operatorname{span}(v_1, \dots, v_k)$ で、さらに最良近似の等号条件から $b_i = WW^{\top}a_i$ 、つまり $B = AV_kV_k^{\top} = A_k$ でなければならない。 $\sigma_k = \sigma_{k+1}$ なら特異ベクトルの選び方で $A_k$ 自体が変わるので一意でない（ $A = I_2$ , $k = 1$ では単位ベクトル $u$ による $uu^{\top}$ がすべて最小点）。作用素ノルムでは $\sigma_k > \sigma_{k+1}$ でも一意とは限らない（ $A = \operatorname{diag}(3, 2, 1)$ , $k = 1$ では $\operatorname{diag}(3 + t, 0, 0)$ （ $\lvert t \rvert \leq 2$ ）の誤差はすべて $\sigma_2 = 2$ ）。

ヒント

実務では (1) 主成分分析は $X_c$ の特異値分解（np.linalg.svd など）で計算する。 $X_c^{\top}X_c$ の固有値分解は条件数を 2 乗にして小さい特異値を失う。 $\pm(1, 1)$ , $\pm(\delta, 0)$ , $\pm(0, \delta)$ の 6 行の $X_c$ （ $\delta = 10^{-8}$ ）では、 $X_c^{\top}X_c$ の対角成分 $2 + 2\delta^2$ が倍精度で $2$ に丸められて第 2 固有値が $0$ になるが、特異値分解は $\sigma_2 = \sqrt{2}\delta$ を正しく返す。(2) 平均と主成分方向は訓練データだけで求めてテストデータに適用する。全データで求めると前処理にテストデータの情報が漏れる（第1章 1.8 節）。(3) 固有ベクトルの符号は任意で、ライブラリや計算法、データのわずかな違いで反転しうる。

3.4 寄与率と次元の選び方

定義 3.11（寄与率）第 $j$ 主成分の寄与率 (proportion of variance explained) を $\lambda_j/(\lambda_1 + \cdots + \lambda_d)$ 、第 $k$ 主成分までの累積寄与率を $(\lambda_1 + \cdots + \lambda_k)/(\lambda_1 + \cdots + \lambda_d)$ と定める。

定理 3.5 より、累積寄与率は最良の $k$ 次元部分空間への射影で残る分散の割合、すなわち $1 -$ （平均二乗再構成誤差）/（全分散）である。例 3.6 の第 1 主成分の寄与率は $0.9$ である。

次元 $k$ は、累積寄与率が 80% や 90% を超える最小の $k$ 、固有値のグラフ（スクリープロット）の「肘」、相関行列の固有値が $1$ 以上の成分の数などで選ばれることが多いが、これらは経験則であって定理ではない。後に予測などの目的があるなら、 $k$ ごとの性能を交差検証（第1章 1.7 節）で比べるのがよい。

注意 3.12（単位への依存）主成分分析は特徴量の単位に依存する。ある特徴量をメートルからセンチメートルに変えると分散は $10^4$ 倍になり、第 1 主成分はその軸に引き寄せられる（問題 3.3）。単位の違う特徴量は、標準偏差で割って標準化してから（相関行列について）主成分分析することが多い。ただし標準化は雑音のような小さい特徴量も同じ重みに引き上げるので、同種のセンサーの値のように大きさを比べられる場合は標準化しないこともある。

注意

寄与率の小さい主成分を「ノイズ」として捨ててよいとは限らない。主成分分析は目的変数を見ないので、予測に効く方向の分散が小さいこともある。2 クラスのデータが $x_1$ 方向に大きく広がり、クラスの違いが $x_2$ の符号（ばらつき $\pm 0.1$ 程度）にだけ現れるなら、第 2 主成分を捨てると分類はできない。

3.5 白色化

主成分得点は無相関だが分散 $\lambda_j$ はまちまちである。さらに標準偏差で割って、すべての方向の分散を $1$ にそろえる前処理を白色化という。

定義 3.13（白色化） $\hat{\Sigma}$ を正定値とする。 $d$ 次正方行列 $W$ による $z_i = W\tilde{x}_i$ の標本共分散行列 $W\hat{\Sigma}W^{\top}$ が $I_d$ になるとき、 $W$ を白色化 (whitening) 行列という。 $\hat{\Sigma} = V\Lambda V^{\top}$ （ $\Lambda = \operatorname{diag}(\lambda_1, \dots, \lambda_d)$ ）として、 $W = \Lambda^{-1/2}V^{\top}$ を PCA 白色化、 $W = \hat{\Sigma}^{-1/2} = V\Lambda^{-1/2}V^{\top}$ を ZCA 白色化という（ $\hat{\Sigma}^{1/2}$ は 02-linear-algebra 第8章命題 8.15 の正の平方根）。

命題 3.14 $\hat{\Sigma}$ を正定値とする。

$W$ が白色化行列であるための必要十分条件は、ある直交行列 $Q$ で $W = Q\hat{\Sigma}^{-1/2}$ と書けることである。

白色化行列 $W$ について $\lVert Wx - Wy \rVert^2 = (x - y)^{\top}\hat{\Sigma}^{-1}(x - y)$ 。右辺の平方根をマハラノビス距離 (Mahalanobis distance) という。

証明. 1： $W = Q\hat{\Sigma}^{-1/2}$ なら $W\hat{\Sigma}W^{\top} = QQ^{\top} = I$ 。逆に $W\hat{\Sigma}W^{\top} = I$ なら $Q = W\hat{\Sigma}^{1/2}$ は $QQ^{\top} = W\hat{\Sigma}W^{\top} = I$ を満たす直交行列で、 $W = Q\hat{\Sigma}^{-1/2}$ 。2： $W^{\top}W = \hat{\Sigma}^{-1/2}Q^{\top}Q\hat{\Sigma}^{-1/2} = \hat{\Sigma}^{-1}$ 。なお PCA 白色化は $Q = V^{\top}$ の場合である。 $\square$

例 3.15 例 3.6 では $\hat{\Sigma}^{-1}$ は第 1 行 $(5, -4)$ 、第 2 行 $(-4, 5)$ の行列の $5/18$ 倍である。新しい点 $(4, 4)$ と $(5, 7)$ の平均 $(3, 5)$ からのユークリッド距離は $\sqrt{2} \approx 1.41$ と $2\sqrt{2} \approx 2.83$ だが、マハラノビス距離は $\sqrt{5} \approx 2.24$ と $2\sqrt{5}/3 \approx 1.49$ で、大小が逆転する。 $(4, 4)$ はデータがほとんど散らばらない $v_2$ 方向にずれているからで、異常検知でマハラノビス距離が使われるのはこのためである。

白色化は小さい $\lambda_j$ で割るので雑音を大きく増幅する。実際には $(\Lambda + \varepsilon I)^{-1/2}$ （ $\varepsilon > 0$ は小さい定数）を使うか、小さい固有値の成分を捨ててから白色化する。

3.6 低ランク近似と推薦

$m$ 人の利用者が $n$ 個の商品につけた評価の行列 $R = (r_{ij})$ を考える。好みが少数の要因で決まるなら、 $k$ 次元のベクトル $p_i, q_j$ で $r_{ij} \approx p_i^{\top}q_j$ 、すなわち $R \approx PQ^{\top}$ と書けるはずである。 $R$ がすべてわかっていれば、 $\lVert R - PQ^{\top} \rVert_F$ の最小は打ち切った特異値分解 $R_k$ で達成される（定理 3.10）。しかし実際の評価行列は大部分が欠損している。観測された添字の集合を $\Omega$ 、利用者 $i$ が評価した商品の集合を $\Omega_i$ として、解くべき問題は

\min_{P, Q} \sum_{(i, j) \in \Omega} (r_{ij} - p_i^{\top}q_j)^2 + \lambda\left(\lVert P \rVert_F^2 + \lVert Q \rVert_F^2\right)

であり、特異値分解はそのままでは使えない。欠損を $0$ で埋めて特異値分解すると、「評価していない」を「評価 $0$ 」として近似してしまう。

例 3.16 $3 \times 3$ のランク 1 の行列（第 $i$ 行が $i \cdot (1, 2, 3)$ ）の $(1, 3)$ 成分 $r_{13} = 3$ だけが欠損しているとする。ランク 1 では 2 次の小行列式が $0$ なので $r_{12}r_{23} = r_{13}r_{22}$ となり、観測値と整合する値は $r_{13} = 2 \cdot 6/4 = 3$ だけである。ところが $0$ で埋めた行列の最良ランク 1 近似では $(1, 3)$ 成分の予測は約 $1.08$ で、観測済みの $r_{11} = 1$ まで約 $0.37$ に引き下げられる。観測値の平均で埋めても約 $3.77$ である。上の問題を $k = 1$ , $\lambda = 0$ として観測値だけで解けば、誤差 $0$ で $r_{13} = 3$ を得る。

注意 3.17（非凸性と交互最小二乗法）上の目的関数は $(P, Q)$ について凸でない（ $1 \times 1$ の $(r - pq)^2$ （ $r \neq 0$ ）でも、原点は勾配が $0$ でヘッセ行列の固有値が $\pm 2r$ の鞍点である）。よく使われる交互最小二乗法 (alternating least squares, ALS) では、 $Q$ を固定して各 $p_i$ をリッジ回帰（第2章定理 2.5）の解
$p_i = \Bigl(\sum_{j \in \Omega_i} q_jq_j^{\top} + \lambda I\Bigr)^{-1}\sum_{j \in \Omega_i} r_{ij}q_j$
に更新し、次に $P$ を固定して $Q$ を同様に更新する。各段階でその変数について厳密に最小化するので目的関数は増えないが（定理 3.19 と同じ論法）、大域最適解に到達する保証はない。

ヒント

実務では 推薦の欠損はランダムではない（利用者は気に入りそうな商品を選んで評価する）ので、欠損を $0$ とみなしても観測値だけに当てはめても偏りが生じうる。クリックや購入しか記録がなければ「反応なし」を重みの小さい負例として扱うなどの工夫をする。評価は観測値の一部を隠して行い、時間順のデータでは過去で学習して未来で評価する。観測のない新しい利用者・商品は低ランクモデルだけでは予測できない（コールドスタート問題）。

3.7 $k$ 平均法

顧客を購買傾向の似たグループに分けたい。各グループの中心とのずれの 2 乗和を小さくする分け方を探すのが $k$ 平均法である。

定義 3.18（ $k$ 平均法）割り当て $c\colon \lbrace 1, \dots, n \rbrace \to \lbrace 1, \dots, k \rbrace$ と中心 $\mu_1, \dots, \mu_k \in \mathbb{R}^d$ について $J(c, \mu) = \sum_{i=1}^n \lVert x_i - \mu_{c(i)} \rVert^2$ とする。ロイドのアルゴリズム (Lloyd's algorithm) は、初期中心から次の 2 段階を繰り返し、(a) で割り当てが変わらなくなったら止まる。

(a) 割り当て：各 $i$ を最も近い中心に割り当てる（同点なら今の割り当て先を優先し、次に番号の小さい中心を選ぶ）。

(b) 更新： $C_j = \lbrace i \mid c(i) = j \rbrace$ が空でなければ $\mu_j$ を $\lbrace x_i \rbrace_{i \in C_j}$ の平均にする（空なら変えない）。

定理 3.19（ $k$ 平均法の単調性）ロイドのアルゴリズムの各段階で $J$ は増えず、アルゴリズムは有限回で止まる。止まったとき、空でない各クラスタの中心はそのクラスタの平均で、各点は最も近い中心に割り当てられている。

証明. (a)： $\mu$ を固定すると $J$ は項 $\lVert x_i - \mu_{c(i)} \rVert^2$ の和なので項ごとの最小化で増えず、割り当てが 1 つでも変われば（厳密に近い中心へ移るので）厳密に減る。(b)： $c$ を固定すると $J = \sum_j \sum_{i \in C_j} \lVert x_i - \mu_j \rVert^2$ で、補題 3.2 より内側の和は $\mu_j$ が $C_j$ の平均のとき最小である。有限性：(b) の直後の $J$ は、空でないクラスタの平均までの距離の 2 乗和 $F(c)$ に等しく、割り当て $c$ だけで決まる。止まらずに $c$ が $c'$ に変われば、(a) で厳密に減り (b) で増えないので $F(c') < F(c)$ 。よって同じ割り当ては 2 度現れず、割り当ては高々 $k^n$ 通りなので有限回で止まる。最後の主張は、止まる直前の (b) と最後の (a) から従う。 $\square$

保証されるのは止まることと $J$ が減ることだけで、止まった点が最小点とは限らない。

例 3.20（局所解に止まる例）4 点 $(0, 0), (0, 1), (4, 0), (4, 1)$ を $k = 2$ で分ける。初期中心を $\mu_1 = (0, 0)$ , $\mu_2 = (0, 1)$ とすると、(a) で $(0, 0), (4, 0)$ が中心 1 に、 $(0, 1), (4, 1)$ が中心 2 に割り当てられ（ $J = 32$ ）、(b) で $\mu_1 = (2, 0)$ , $\mu_2 = (2, 1)$ 、 $J = 16$ となって止まる。しかし左右に分ければ $J = 4 \cdot (1/2)^2 = 1$ で（空でない 2 組への 7 通りの分け方の中で最小）、初期中心を $(0, 0), (4, 0)$ にとればこちらに到達する。

$k$ 平均法の大域最小点を求める問題は一般に NP 困難であることが知られている（本書では証明しない。NP 困難は 23-optimization 第7章で扱う）。実際には、初期中心を互いに離れるように確率的に選ぶ k-means++ などの初期化と、複数回の実行で $J$ が最小のものを採ることが行われる。ユークリッド距離を使うので特徴量の単位に依存し、最適な $J$ は $k$ について増えない（ $k = n$ なら $0$ ）ので $J$ の大小で $k$ は選べない。

3.8 ジョンソン–リンデンシュトラウスの補題

$d = 10^4$ 次元のベクトルが $n = 10^6$ 個あり、近いものを探したいとする。主成分分析は捨てた方向だけで異なる 2 点を同じ点に写すので、個々の距離を保証しない。ここではすべての 2 点間の距離をほぼ保つことを求める。驚くべきことに、データを見ずに選んだランダムな線形写像で、 $d$ によらない次元までそれができる。

定理 3.21（ジョンソン–リンデンシュトラウスの補題, Johnson–Lindenstrauss lemma） $0 < \varepsilon < 1$ 、 $x_1, \dots, x_n \in \mathbb{R}^d$ （ $n \geq 2$ ）とし、正の整数 $k$ が
$k \geq \frac{4\ln n}{\varepsilon^2/2 - \varepsilon^3/3}$
を満たすとする。成分が独立に $N(0, 1/k)$ に従う $k \times d$ のランダム行列 $G$ について、確率 $1/n$ 以上で、すべての $i, j$ で
$(1 - \varepsilon)\lVert x_i - x_j \rVert^2 \leq \lVert Gx_i - Gx_j \rVert^2 \leq (1 + \varepsilon)\lVert x_i - x_j \rVert^2 \tag{1}$
が成り立つ。特に (1) を満たす線形写像 $\mathbb{R}^d \to \mathbb{R}^k$ が存在する。 $\varepsilon^2/2 - \varepsilon^3/3 \geq \varepsilon^2/6$ なので $k \geq 24\varepsilon^{-2}\ln n$ なら十分で、 $k = O(\varepsilon^{-2}\log n)$ でよい。

証明. $a = \varepsilon^2/2 - \varepsilon^3/3$ とおき、 $u = x_i - x_j \neq 0$ を固定する（ $u = 0$ なら (1) は自明）。 $G$ の第 $l$ 行 $g_l^{\top}$ について、 $g_l^{\top}u$ は独立な正規変数の 1 次結合なので $N(0, \lVert u \rVert^2/k)$ に従い、 $G$ の別々の行から作られるので $l$ について独立である（22-statistics 第1章系 1.17 の 1、命題 1.4 の 3）。よって $Y = k\lVert Gu \rVert^2/\lVert u \rVert^2$ は $k$ 個の独立な標準正規変数の 2 乗和で、(1) は $(1 - \varepsilon)k \leq Y \leq (1 + \varepsilon)k$ と同値である。 $Z \sim N(0, 1)$ , $t < 1/2$ について $E[e^{tZ^2}] = \frac{1}{\sqrt{2\pi}}\int e^{-(1 - 2t)z^2/2}\ dz = (1 - 2t)^{-1/2}$ なので、独立な和のモーメント母関数は積であること（22-statistics 第1章 1.4 節）から $E[e^{tY}] = (1 - 2t)^{-k/2}$ 。マルコフの不等式（同定理 1.25）を $e^{tY}$ に使い $t = \varepsilon/(2(1 + \varepsilon))$ とおくと

P(Y \geq (1 + \varepsilon)k) \leq \frac{E[e^{tY}]}{e^{t(1 + \varepsilon)k}} = \left((1 + \varepsilon)e^{-\varepsilon}\right)^{k/2} \leq e^{-ka/2}

（最後は $\ln(1 + \varepsilon) \leq \varepsilon - \varepsilon^2/2 + \varepsilon^3/3$ による。差は $\varepsilon = 0$ で $0$ 、導関数は $\varepsilon^3/(1 + \varepsilon) \geq 0$ ）。同様に $e^{-tY}$ と $t = \varepsilon/(2(1 - \varepsilon))$ から $P(Y \leq (1 - \varepsilon)k) \leq ((1 - \varepsilon)e^{\varepsilon})^{k/2} \leq e^{-k\varepsilon^2/4} \leq e^{-ka/2}$ （ $\ln(1 - \varepsilon) \leq -\varepsilon - \varepsilon^2/2$ と $\varepsilon^2/2 \geq a$ による）。各組で (1) が破れる確率は $2e^{-ka/2}$ 以下で、 $k \geq 4\ln n/a$ より $e^{-ka/2} \leq n^{-2}$ だから、 $n(n - 1)/2$ 組のどこかで破れる確率は $n(n - 1)e^{-ka/2} \leq 1 - 1/n$ 以下であり、すべての組で (1) が成り立つ確率は $1/n$ 以上である。 $\square$

同じ計算で、 $0 < \delta < 1$ について $k \geq (4\ln n + 2\ln(1/\delta))/a$ なら $e^{-ka/2} \leq \delta n^{-2}$ となり、すべての組で (1) が成り立つ確率は $1 - \delta$ 以上になる。

$k$ は $n$ の対数にしか依存せず、元の次元 $d$ にはよらない。写像はデータを見ずに選べるので、後から来た点にも同じ $G$ を使える。

import numpy as np

rng = np.random.default_rng(0)
n, d = 200, 10_000
X = rng.standard_normal((n, d))       # 200 点（10,000 次元）
iu = np.triu_indices(n, 1)            # 点の組 i < j

def sqdist(A):                        # すべての組の距離の 2 乗
    sq = (A**2).sum(axis=1)
    return (sq[:, None] + sq[None, :] - 2 * A @ A.T)[iu]

D = sqdist(X)
for k in [50, 200, 1000]:
    G = rng.standard_normal((k, d)) / np.sqrt(k)   # 成分は独立に N(0, 1/k)
    r = sqdist(X @ G.T) / D
    print(f"k={k:4d}: 比の最小 {r.min():.3f}, 最大 {r.max():.3f}")

k=  50: 比の最小 0.406, 最大 1.872
k= 200: 比の最小 0.638, 最大 1.369
k=1000: 比の最小 0.831, 最大 1.176

$n = 200$ で定理 3.21 が要求する次元は $\varepsilon = 0.5$ なら $255$ 、 $\varepsilon = 0.2$ なら $1223$ で、実験では $k = 1000$ で比が $[0.83, 1.18]$ に収まった。

注意 3.22（主成分分析との比較）主成分分析はデータに合わせて部分空間を選び平均二乗誤差の意味で最適だが、個々の距離は保証しない。ジョンソン–リンデンシュトラウスの補題はデータによらない写像で全組の距離を保証する代わりに、必要な次元が大きくなりうる（ $n = 10^6$ , $\varepsilon = 0.1$ では定理 3.21 の要求は $k \geq 11842$ で、3.8 節の冒頭の $d = 10^4$ を超える。定理 3.21 は十分条件なので、実際にはもっと小さい $k$ で足りることも多い）。データが低次元の部分空間の近くにあるなら主成分分析が、そうでない高次元データの近傍探索ではランダム射影が向く。

まとめ

主成分方向は標本共分散行列の固有ベクトル、分散は固有値である（レイリー商）。第 $k$ 主成分までの部分空間は射影の分散の和を最大にし、再構成誤差を最小にする。最良の近似アフィン部分空間は平均を通るので、中心化が必要である。
主成分分析は中心化データ行列の特異値分解そのもので、分散は $\sigma_j^2/n$ である。 $X_c^{\top}X_c$ は作らずに計算する。
エッカート–ヤングの定理：打ち切った特異値分解は、フロベニウスノルムでも作用素ノルムでも最良のランク $k$ 近似である。
累積寄与率は最良の $k$ 次元部分空間に残る分散の割合である。次元の選び方の基準は経験則で、主成分分析は単位に依存する。
白色化行列は $Q\hat{\Sigma}^{-1/2}$ （ $Q$ は直交行列）に限られ、白色化後の距離はマハラノビス距離である。
欠損の多い評価行列に特異値分解はそのまま使えない。観測値への当てはめは非凸である。
$k$ 平均法は目的関数を単調に減らして有限回で止まるが、局所解に止まりうる。
ジョンソン–リンデンシュトラウスの補題：ランダムな線形写像で、 $n$ 点の距離の 2 乗の比をすべて $1 \pm \varepsilon$ の範囲に保ったまま $O(\varepsilon^{-2}\log n)$ 次元に写せる。

演習問題

問題 3.1 ★ 4 点 $(0, 0), (2, 2), (4, 4), (6, 2)$ について、標本共分散行列（ $n$ で割る）、主成分方向、各主成分の分散と寄与率、第 1 主成分得点を求めよ。

解答

平均は $(3, 2)$ 、中心化した点は $(-3, -2), (-1, 0), (1, 2), (3, 0)$ で、 $\hat{\Sigma}$ は第 1 行 $(5, 2)$ 、第 2 行 $(2, 2)$ の行列である。固有多項式は $t^2 - 7t + 6 = (t - 6)(t - 1)$ なので $\lambda_1 = 6$ , $\lambda_2 = 1$ 、主成分方向は $v_1 = (2, 1)/\sqrt{5}$ , $v_2 = (1, -2)/\sqrt{5}$ （符号は任意）、寄与率は $6/7 \approx 0.857$ と $1/7$ 。第 1 主成分得点は $(-8, -2, 4, 6)/\sqrt{5}$ で、2 乗の平均は $(64 + 4 + 16 + 36)/20 = 6 = \lambda_1$ に一致する。

問題 3.2 ★ 第 1 行 $(3, 0)$ 、第 2 行 $(4, 5)$ の行列 $A$ の特異値分解（02-linear-algebra 第8章例 8.20： $\sigma_1 = 3\sqrt{5}$ , $u_1 = (1, 3)/\sqrt{10}$ , $v_1 = (1, 1)/\sqrt{2}$ , $\sigma_2 = \sqrt{5}$ ）から最良のランク 1 近似 $A_1$ を求め、 $\lVert A - A_1 \rVert_F^2 = \sigma_2^2$ を確かめよ。第 1 行を $0$ にした行列 $B$ の誤差と比べよ。

解答

$A_1 = \sigma_1u_1v_1^{\top}$ は $\frac{3\sqrt{5}}{\sqrt{20}} = \frac{3}{2}$ より第 1 行 $(1.5, 1.5)$ 、第 2 行 $(4.5, 4.5)$ の行列。 $A - A_1$ の成分は $1.5, -1.5, -0.5, 0.5$ で 2 乗和は $5 = \sigma_2^2$ 。一方 $\lVert A - B \rVert_F^2 = 9 > 5$ で、定理 3.10 のとおり $A_1$ のほうが近い。

問題 3.3 ★★ 例 3.6 のデータで第 2 特徴量だけを 10 倍したとき、標本共分散行列、第 1 主成分の寄与率、第 1 主成分方向が第 1 軸となす角を求めよ。元の方向 $(1, 1)$ をこの座標で表したものと比べ、各特徴量を標準化した場合とも比べよ。

解答

$\hat{\Sigma}' = \operatorname{diag}(1, 10)\ \hat{\Sigma}\ \operatorname{diag}(1, 10)$ は第 1 行 $(2, 16)$ 、第 2 行 $(16, 200)$ で、固有値は $101 \pm \sqrt{10057}$ （約 $201.28$ と $0.715$ ）、寄与率は約 $0.9965$ （元は $0.9$ ）。第 1 固有ベクトルは $(16, \lambda_1 - 2) \approx (16, 199.28)$ の方向で、第 1 軸と約 $85.4^\circ$ をなす。元の方向 $(1, 1)$ はこの座標では $(1, 10)$ の方向（約 $84.3^\circ$ ）なので、主成分は同じ方向の座標変換にはならず、分散の大きい特徴量に引き寄せられ、寄与率も見かけ上上がる。標準化すると、どちらの単位でも相関行列は第 1 行 $(1, 0.8)$ 、第 2 行 $(0.8, 1)$ （固有値 $1.8, 0.2$ ）で単位によらない。

問題 3.4 ★★ 1 次元のデータ $0, 1, 4, 5, 8, 9$ を $k = 3$ で分ける。(1) 初期中心を $4, 8, 9$ としてロイドのアルゴリズムを実行し、止まったときの $J$ を求めよ。(2) $J$ の最小値が $1.5$ であることを示せ。

解答

(1) (a) で $0, 1, 4, 5$ が中心 $4$ に（ $5$ から $4$ までは $1$ 、 $8$ までは $3$ ）、 $8, 9$ はそれぞれ自分の中心に割り当てられ、(b) で中心は $2.5, 8, 9$ 。次の (a) でも $5$ は $2.5$ のほうが近く（ $2.5 < 3$ ）割り当ては変わらないので止まり、 $J = 2.5^2 + 1.5^2 + 1.5^2 + 2.5^2 = 17$ 。

(2) $\lbrace 0, 1 \rbrace, \lbrace 4, 5 \rbrace, \lbrace 8, 9 \rbrace$ で $J = 6 \cdot 0.5^2 = 1.5$ 。下界：最小 $a$ 、最大 $b$ のクラスタの寄与は、平均 $m$ について $(a - m)^2 + (b - m)^2 \geq (b - a)^2/2$ 以上である。空のクラスタがあるときは、2 点以上のクラスタの 1 点を移しても $J$ は増えないので、大きさが $(2, 2, 2)$ , $(3, 2, 1)$ , $(4, 1, 1)$ の場合を調べればよい。2 点の幅は $1$ 以上、3 点の幅は $4$ 以上、4 点の幅は $5$ 以上なので、それぞれ $J \geq 1.5$ , $J \geq 8$ , $J \geq 12.5$ 。よって最小値は $1.5$ で、(1) は局所解である。

問題 3.5 ★★ 顧客の解約予測（特徴量 300 個）について次の報告があった。「全データを標準化して主成分分析し、累積寄与率 90% となる 25 成分を残した。その後データを訓練用とテスト用に分けてロジスティック回帰を学習し、テストの正解率は 91% だった。残りの 275 成分は寄与率が小さいのでノイズである。」問題点を指摘せよ。

解答

(1) 標準化と主成分分析をテストデータを含む全データで行っており、データ漏洩である（第1章 1.8 節）。訓練データだけで求めた変換をテストデータに適用すべきである（目的変数を使わないので影響は小さいことも多いが、手順として誤り）。(2) 累積寄与率 90% は経験則で、成分数は訓練データ内の交差検証で選ぶべきである。(3) 分散の小さい方向が解約の予測に効くこともある（3.4 節の WARNING）。「ノイズ」と断定する根拠はなく、捨てた成分を含めたモデルと比べる必要がある。

問題 3.6 ★★ 評価行列 $R$ の第 1 行が $(2, 4)$ 、第 2 行が $(3, ?)$ で、 $(2, 2)$ 成分が欠損している。(1) ランク 1 の行列として観測値と整合する $r_{22}$ を求めよ。(2) 欠損を $0$ で埋めた最良ランク 1 近似の $(2, 2)$ 成分は約 $1.47$ である（計算機で確かめよ）。なぜ (1) から大きくずれるのか。(3) $k = 1$ , $\lambda = 0$ の交互最小二乗法で $q = (q_1, q_2)$ を固定したときの $p_1, p_2$ の更新式を書き、目的関数が増えない理由を述べよ。

解答

(1) ランク 1 なら $r_{11}r_{22} = r_{12}r_{21}$ なので $r_{22} = 4 \cdot 3/2 = 6$ 。(2) 埋めた $0$ を含む全成分との 2 乗誤差を最小にするので、近似は $(2, 2)$ 成分を $0$ に近づけようとし、ほかの成分もゆがむ（ $(1, 1)$ 成分は約 $2.78$ になる）。「未評価」を「評価 0」と扱った偏りである。(3) 観測は $(1, 1), (1, 2), (2, 1)$ なので $p_1 = (2q_1 + 4q_2)/(q_1^2 + q_2^2)$ , $p_2 = 3/q_1$ （ $q_1 \neq 0$ ）。これは $q$ を固定したときの $\sum_{(i, j) \in \Omega}(r_{ij} - p_iq_j)^2$ の各 $p_i$ についての厳密な最小化（1 変数の最小二乗）なので、値は増えない。 $q$ の更新も同様である。 $q = (1, 2)$ なら $p = (2, 3)$ で誤差 $0$ 、予測は $p_2q_2 = 6$ となる。

問題 3.7 ★★ (1) $n = 10^4$ , $\varepsilon = 0.2$ のとき定理 3.21 が要求する次元 $k$ を求めよ。(2) 標準基底 $e_1, \dots, e_n \in \mathbb{R}^d$ （ $n \leq d$ ）を、 $d$ 個の座標から $k$ 個を選んで残し定数倍する写像 $f$ で写す。 $k \leq n - 2$ なら、座標の選び方と定数によらず定理 3.21 の (1) 式が破れる組があることを示せ。

解答

(1) $a = 0.02 - 0.008/3 \approx 0.017333$ 、 $4\ln 10^4/a \approx 36.841/0.017333 \approx 2125.5$ なので $k = 2126$ 。(2) 残す座標に含まれない $e_i$ は $0$ に写る。 $k \leq n - 2$ なら少なくとも 2 つの $e_i, e_j$ が $0$ に写り、 $\lVert f(e_i) - f(e_j) \rVert^2 = 0 < 2(1 - \varepsilon) = (1 - \varepsilon)\lVert e_i - e_j \rVert^2$ 。座標を選ぶ方法は疎なデータに弱いが、成分が正規分布に従う $G$ ならデータによらず定理 3.21 が成り立つ。