主成分分析（PCA; Principal Component Analysis）を丁寧に解く

この間書いた記事で，codecogsのOnline Latex Equation Editorを使えば，はてなダイアリーでもキレイな数式を書けることが分かったので，昔書いたPCAの導出を載せておこうと思う．
PCAは固有値分解で主成分を抽出する写像を求めることができるんだけど，多くの教科書などでは，「PCAの最適化問題が固有値問題になるよー」と書いて終わりなので，「どうして固有値問題で解くことができるのか」というところを丁寧に導出したい．
考えた結果，すごい回りくどい方法しか思いつかなかったので，他のシンプルな解説もあるのかもしれない．
特に直感的な説明などなくて式を丁寧に見ていっているだけです．たぶん間違ってたり抜けているところがあるので，見つけた人は優しくご指摘下さい．（言い訳がすごい）

主成分分析（PCA）の目的と定式化

次の観測されたN個の観測データから成るデータ集合が与えられたとする．

$\mathbf{x}_1, \cdots, \mathbf{x}_N$ : 観測値
$D$ : 各 $\mathbf{x}_i$ の次元
$K$ : $\{\mathbf{x}_1, \cdots, \mathbf{x}_N\}$ 内の一次独立であるベクトル数（別に意味は分からくてもいい）． $K \leq N$ が成り立つはず．

この観測値の次元 $D$ がすごく大きい場合，すごく扱いづらい．
なので，観測値 $\{\mathbf{x}_n\}_{n=1}^N$ の情報を出来るだけ維持しながら，次元を小さくしたい．
この次元削減を次のような式で実現する．

$\mathbf{t}_n = \mathbf{U}^T \mathbf{x}_n, \quad n = 1, \ldots, N.$

ここで，

$\mathbf{U} \in \mathbb{R}^{D \times R}$ : 次元 $D$ のベクトルを次元 $R$ に小さくする行列．
$R \leq K \leq D$ （ $R$ はもちろん $D$ より小さくする． $R=D$ の場合は次元削減の点で言えば全然意味ないんだけど，そういうのを求めたい変態も少なからずいる．また $\{\mathbf{x}_1, \cdots, \mathbf{x}_N\}$ 内の一次独立であるベクトル数 $K$ よりも小さくしたほうがいい）
$\mathbf{t}_n$ : 行列 $\mathbf{U}$ によって $\mathbf{x}_n$ を次元削減した結果．
$\mathbf{U}^T\mathbf{U} = \mathbf{I}_R$ （ $\mathbf{I}_R$ は $R \times R$ の単位行列． $\mathbf{U}$ の列ベクトルは互いに直交する．理由はよく分からないけど，PCAはこういう写像を求めたいんだから仕方ない）

$\mathbf{U}$ は，「観測値 $\{\mathbf{x}\}_{i=1}^N$ の情報を出来るだけ維持しながら，次元を小さくする」写像である．
この「情報を出来るだけ維持する」ことを評価するために，ここでは，平均二乗誤差（MSE: Mean Squared Error）を小さくすることを考える．
この誤差は， $\mathbf{x}_n$ と観測値の空間 $\mathbb{R}^D$ で $\mathbf{t}_n$ を見たときの誤差のことである．
わかりづらいので，2次元の観測値を1次元に次元削減する問題

${t}_n = \mathbf{u}^T \mathbf{x}_n, \quad n = 1, \ldots, N.$

を例にとって説明してみる．
下のクソみたいな図は， $\mathbf{u}$ による次元削減を元の空間（観測値がある空間）で表したものである．

何が言いたいかというと，

元の空間上でみると，次元削減後の観測値 $t_n$ がある空間というのは，基底 $\mathbf{u}$ で張られた空間（図では直線で表される1次元空間）である．
$t_n$ は，その $\mathbf{u}$ を基底としたときの座標である．つまり， $t_n$ というのは，元の空間でみると，ベクトル $\mathbf{u}t_n$ の長さのこと．
したがって誤差ベクトルというのは次のように定義できる．

$\mathbf{x}_n - \mathbf{u}t_n$

全観測値の誤差ベクトルの長さ（の二乗）を平均を取ると，

$E(\mathbf{U}) = \frac{1}{N} \sum_{n=1}^{N} \| \mathbf{x}_{n} - \mathbf{u} t_n \|^{2}$

このような感じで，PCAの問題は最適化問題として定式化される．

$\min_{\mathbf{U} \in \mathbb{R}^{D \times R}} & \quad E(\mathbf{U}) = \frac{1}{N} \sum_{n=1}^{N} \| \mathbf{x}_{n} - \mathbf{U} \mathbf{t}_{n} \|^{2} \\ \text{subject to} & \quad \mathbf{U}^T \mathbf{U} = \mathbf{I}_R.$

最適化問題を解く

上の $E(\mathbf{U})$ を最大にする $\mathbf{U}$ を求めることを考える．
$E(\mathbf{U})$ は次のように変形できる．

$\begin{align} E(\mathbf{U}) &= \frac{1}{N} \sum_{n=1}^{N} \| \mathbf{x}_{n} - \mathbf{U} \mathbf{t}_{n} \|^{2} \nonumber \\ &= \frac{1}{N} \sum_{n=1}^{N} \| \mathbf{x}_{n} - \mathbf{U} \mathbf{U}^{T} \mathbf{x}_{n} \|^{2} \nonumber \\ &= \frac{1}{N} \sum_{n=1}^{N}(\mathbf{x}_{n} - \mathbf{U}\mathbf{U}^{T} \mathbf{x}_{n})^{T}(\mathbf{x}_{n} - \mathbf{U}\mathbf{U}^{T} \mathbf{x}_{n}) \nonumber \\ &= \frac{1}{N} \sum_{n=1}^{N} (- 2\mathbf{x}_{n}^{T} \mathbf{U}\mathbf{U}^{T}\mathbf{x}_{n}+ \mathbf{x}_{n}^{T}\mathbf{U} \mathbf{I}_{R} \mathbf{U}^{T}\mathbf{x}_{n} ) \nonumber \\ &= \frac{1}{N} \sum_{n=1}^{N}{\rm trace}(\mathbf{U}^{T} \mathbf{x}_{n} \mathbf{x}_{n}^{T} \mathbf{U}) \nonumber \\ &= {\rm trace}(\mathbf{U}^{T} \mathbf{S} \mathbf{U}). \nonumber \end{align}$

ここで， $\mathbf{S} \in \mathbb{R}^{D \times D}$ は，相関係数行列で

$\mathbf{S} = \sum_{i=1}^{N} \mathbf{x}_n \mathbf{x}_n^T$

と定義される．

$\mathbf{S}$ を固有ベクトルと固有値によって対角化する．

$\mathbf{S} = \sum_{i = 1}^{D} \lambda_i \mathbf{v} \mathbf{v}^T = \mathbf{V} \mathbf{D} \mathbf{V}^T$

$\mathbf{v}_i$ : 固有値 $\lambda_i$ に対応する固有ベクトル
$\lambda_1 \geq \lambda_2 \geq \cdots \geq \lambda_D \geq 0$ とする．
$\mathbf{V} = [ \mathbf{v}_1, \ldots, \mathbf{v}_D]$
$\mathbf{D} = \rm{diag}([\lambda_1, \ldots, \lambda_D])$
$\mathbf{S}$ は対称行列なので， $\mathbf{V}^T \mathbf{V} = \mathbf{I}_{D}$ となるように取ることができる（固有ベクトルは互いに直交し，ノルムが1）．

$E(\mathbf{U})$ をこの固有値・固有ベクトルで表すと

$\begin{align} \rm{trace}(\mathbf{U}^{T} \mathbf{S} \mathbf{U}) =& \rm{trace}(\mathbf{U}^{T} \mathbf{V} \mathbf{D} \mathbf{V}^{T} \mathbf{U}) \nonumber \\ =& \lambda_{1} |\mathbf{u}_{1}^{T} \mathbf{v}_{1}|^{2} + \cdots + \lambda_{D} |\mathbf{u}_{1}^{T} \mathbf{v}_{D}|^{2} \nonumber \\ & + \lambda_{1} |\mathbf{u}_{2}^{T} \mathbf{v}_{1}|^{2} + \cdots + \lambda_{D} |\mathbf{u}_{2}^{T} \mathbf{v}_{D}|^{2} + \nonumber \\ & \cdots + \lambda_{1} |\mathbf{u}_{R}^{T} \mathbf{v}_{1}|^{2} + \cdots + \lambda_{D} |\mathbf{u}_{R}^{T} \mathbf{v}_{D}|^{2}. \nonumber \end{align}$

$\{\mathbf{v}_i\}_{i=1}^D$ は， $\mathbb{R}^D$ の正規直交基底なので， $\mathbf{U}$ の列ベクトルをこれらの線形結合で表すことにする．

$\mathbf{u}_i = a_{i1} \mathbf{v}_1 + \cdots + a_{iD}\mathbf{v}_D$

$\mathbf{U} = [ \mathbf{u}_1, \ldots, \mathbf{u}_R]$

したがって，

$|\mathbf{u}_{i}^{T} \mathbf{v}_{j}|^{2} = | (a_{i1} \mathbf{v}_{1} + \cdots + a_{iD} \mathbf{v}_{D})\mathbf{v}_{j}|^{2} = a_{ij}^{2}.$

と書けるので， $E(\mathbf{U})$ は，

$\begin{align} E(\mathbf{U}) = E(\mathbf{A}) =& \lambda_1 a_{11}^2 + \cdots + \lambda_d a_{1D}^2 \nonumber \\ & + \lambda_1 a_{21}^2 + \cdots + \lambda_d a_{2D}^2 \nonumber \\ & \cdots + \lambda_1 a_{R1}^2 + \cdots + \lambda_d a_{RD}^2 \nonumber \end{align}$

と $\mathbf{A}$ の関数として書くことができる．
ここで， $\mathbf{A}$ は $R \times D$ の行列で，その要素は， $[\mathbf{A}]_{ij} = a_{ij}$ である．

さらに，制約条件 $\mathbf{U}^T \mathbf{U} = \mathbf{I}_R$ は，

$\begin{align} \mathbf{u}_{i}^{T} \mathbf{u}_{i} & = \left( \sum_{j,k=1}^{D} a_{ij} a_{ik} \mathbf{v}_{j}^{T} \mathbf{v}_{k} \right) \nonumber \\ & = a_{i1}^{2} + \ldots + a_{iD}^{2} = 1 \nonumber \end{align}$
$\mathbf{u}_i^T\mathbf{u}_j = a_{i1}a_{j1} + \ldots + a_{iD}a_{jD} = 0$

と書くことができる．
以上より，PCAの最適化問題は $\mathbf{U}$ を求める問題から変形され， $\mathbf{A}$ を求める問題として以下のように定式化される．

$\begin{align} \max_{\mathbf{A}} \quad & \sum_{i=1}^{R} \lambda_{1} a_{i1}^{2} + \cdots + \lambda_{D}a_{iD}^{2} = \sum_{i=1}^{R} J_{i} (a_{i1}, \cdots, a_{iD}), \nonumber \\ \text{subject to} \quad & a_{i1}^{2} + \cdots + a_{iD}^{2} = 1, \quad \text{for } i = 1, \ldots, R, \nonumber \\ & a_{i1}a_{j1} + \cdots + a_{iD}a_{jD} = 0, \nonumber \\ & \quad \quad \text{for } i = 2, \ldots, R, \quad j = i+1, \cdots, R \nonumber \end{align}$

ここで，評価関数も $J_i (a_{i1}, \ldots, a_{iD})$ の和の形で書けるため，それぞれで最大化値を求めれば良い．
また，重複する制約条件を省略すると， $\{a_{11}, \ldots, a_{1d} \}$ は，その他のパラメータに依存していないことが分かる．
したがって， $\{a_{11}, \ldots, a_{1D} \}$ は

$\begin{align} \max_{a_{11}, \ldots, a_{1D}} \quad & \lambda_{1} a_{11}^{2} + \cdots + \lambda_{D}a_{1D}^{2}, \nonumber \\ \text{subject to} \quad & a_{11}^{2} + \cdots + a_{1D}^{2} = 1 \nonumber \end{align}$

となる．
制約条件を $a_{11}^2 = 1 - a_{12}^2 - \cdots - a_{1D}^2$ として，これを評価関数に代入すると，

$\begin{align} & \lambda_{1}(1 - a_{12}^{2} - \cdots - a_{1D}^{2}) + \lambda_{2}a_{12}^{2} + \cdots + \lambda_{D} a_{1D}^{2} \nonumber \\ = & \lambda_{1} + (\lambda_{2} - \lambda_{1})a_{12}^{2} + \cdots + (\lambda_{D} - \lambda_{1})a_{1D}^{2} \nonumber \end{align}$

となる．
$\lambda_1 \geq \lambda_2 \geq \cdots \lambda_D \geq 0$ より $0 \geq (\lambda_2-\lambda_1) \geq \cdots (\lambda_D-\lambda_1)$ ， $a_{11}^2 + \cdots + a_{1D}^2$ より $0 \leq a_{1i}^2 \leq 1, i = 1, \ldots, D$ という結果から，明らかに $(\lambda_i - \lambda_1)a_{1i}^2 \leq 0$ である．
したがって，評価関数の最大値は $\lambda_1$ であり，そのときのパラメータの値は，

$a_{11} = 1, \; a_{1i} = 0, \; i = 2, \ldots, D$

となる．

$\{a_{11}, \ldots, a_{1D} \}$ が求まったので，次に， $\{a_{21}, \ldots, a_{2D} \}$ の最適化問題

$\begin{align} \max_{a_{21}, \ldots, a_{2D}} \quad & \lambda_{1} a_{21}^{2} + \cdots + \lambda_{D}a_{2D}^{2}, \nonumber \\ \text{subject to} \quad & a_{21}^{2} + \cdots + a_{2D}^{2} = 1 \nonumber \\ & 1a_{21} + 0a_{22} + \cdots + 0a_{2D} = a_{21} = 0 \nonumber \end{align}$

を考える．
まず，明らかに， $a_{21} = 0$ である．
さらに，制約条件を $a_{22}^2 = 1 - a_{23}^2 - \cdots - a_{2D}^2$ として，これを評価関数に代入すると，

$\begin{align} & \lambda_{2}(1 - a_{23}^{2} - \cdots - a_{3D}^{2}) + \lambda_{3}a_{23}^{2} + \cdots + \lambda_{D} a_{2D}^{2} \nonumber \\ = & \lambda_{2} + (\lambda_{3} - \lambda_{2})a_{12}^{2} + \cdots + (\lambda_{D} - \lambda_{2})a_{2D}^{2} \nonumber \end{align}$