#はじめに
次元削減の主要な手法である主成分分析の理論を理解するために、かなり細かいところまで追いかけて理解したので、皆さんの理解の助けになれば良いなと思って投稿します。
因みに私は機械学習の勉強を始めて3ヶ月ちょっとで、数学の知識も高校レベル程度なので、そういった前知識が無いような方でも理解できると思います。一方で前知識が不足している分、理論が間違っているところもあるかと思います。。
#事前準備
主成分分析の理論を理解する上で、必要な知識や性質がかなり多いので、そういった事前準備をここで行います。事前準備はかなり長いので、問題ない方は飛ばして頂いて大丈夫です。また、証明は読まずに、そういった性質が成り立つことを踏まえて頂くだけでも良いですし、まずは飛ばして、定義や性質を使った時に事前準備の該当箇所に戻ってきても良いかと思います。
ここでは、$X$ を $n$ 次正方行列、$Y$ を $m\times n$ 行列とし、以下のように表すとしておきます。
X =\begin{pmatrix}
x^{(1)}_1 & x^{(2)}_1 & \cdots & x^{(n)}_1\\
x^{(1)}_2 & x^{(2)}_2 & \cdots & x^{(n)}_2\\
\vdots & \vdots & \ddots & \vdots\\
x^{(1)}_n & x^{(2)}_n & \cdots & x^{(n)}_n\\
\end{pmatrix},
Y =\begin{pmatrix}
y^{(1)}_1 & y^{(2)}_1 & \cdots & y^{(n)}_1\\
y^{(1)}_2 & y^{(2)}_2 & \cdots & y^{(n)}_2\\
\vdots & \vdots & \ddots & \vdots\\
y^{(1)}_m & y^{(2)}_m & \cdots & y^{(n)}_m\\
\end{pmatrix}
##固有値と固有ベクトル
ここでは固有値と固有ベクトルについて紹介します。
###定義
$n$ 次正方行列 $X$ とベクトル $\boldsymbol{u}$ に対して、
X\boldsymbol{u} = \lambda \boldsymbol{u}
を満たす $\boldsymbol{0}$ でない $\boldsymbol{u}$ とスカラー $\lambda$ が存在するとき、$\lambda$ を $X$ の固有値、$\boldsymbol{u}$ を $X$ の固有ベクトルといいます。
一般的に行列をベクトルにかけるとベクトルは原点を中心に回転し、大きさも変わります。例えば、$X$ が2次正方行列の場合、
X =\begin{pmatrix}
1 & -2 \\
-3 & 4 \\
\end{pmatrix},
\boldsymbol{x} =
\begin{pmatrix}
1 \\
2 \\
\end{pmatrix}
とすると、
\begin{eqnarray}
X \boldsymbol{x} =
\begin{pmatrix}
1 & -2 \\
-3 & 4 \\
\end{pmatrix}
\begin{pmatrix}
1 \\
2 \\
\end{pmatrix}
=
\begin{pmatrix}
-3 \\
5 \\
\end{pmatrix}
\end{eqnarray}
となります。つまり全く異なるベクトルになってしまいます。簡単に図で表すとこんな感じです。
因みに
X =\begin{pmatrix}
\cos\theta & -\sin\theta \\
\sin\theta & \cos\theta \\
\end{pmatrix}
とすれば、ベクトルの大きさは変わらず、原点を中心に $\theta$ だけ回転させたベクトルになります。
少し話が逸れましたが、簡単に解釈すると、固有ベクトルは、ある行列をかけても大きさが変わるだけで、方向は変わらないベクトルということができます。改めて固有値と固有ベクトルの定義を見てみるとそのことが分かると思います。
X\boldsymbol{u} = \lambda \boldsymbol{u}
###求め方
\bigm\vert
X - \lambda E
\bigm\vert
= 0
を解くことで、固有値を求めることができます。
この $\vert X - \lambda E \vert$ は行列式を表し、この方程式を固有方程式といいます。
なぜこの方程式で固有値が求められるかは、固有値と固有ベクトルの定義から導かれる連立方程式を計算していけば分かります。
###性質1
$n$ 次正方行列は固有値を $n$ 個持つ。
固有方程式を具体的に考えると、左辺には
\bigl(x^{(1)}_1 - \lambda\bigr)\bigl(x^{(2)}_2 - \lambda\bigr)\cdots \bigl(x^{(n)}_n - \lambda\bigr)
という項があります。これは $\lambda$ の $n$ 字式なので、重解があった場合の重複度も考慮すれば、必ず $n$ 個の固有値が見つかります。つまり、$n$ 次正方行列は固有値を $n$ 個持つと言えます。
###性質2
行列 $X$ が対称行列であれば、$X$ の異なる固有値に対応する固有ベクトル同士は直行する。
この性質を証明するために必要な内積の性質があるので、その準備をします。
####内積の性質
\boldsymbol{a}=
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_m
\end{pmatrix},
\boldsymbol{b}=
\begin{pmatrix}
b_1 \\
b_2 \\
\vdots\\
b_n
\end{pmatrix}
とする。このとき、
\boldsymbol{a}\cdot(Y\boldsymbol{b})= (Y^T\boldsymbol{a})\cdot\boldsymbol{b}
が成り立つ。
これを示します。まず、左辺は
\begin{eqnarray}
\boldsymbol{a}\cdot(Y\boldsymbol{b})&=&
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_m
\end{pmatrix}\cdot
\left({
\begin{pmatrix}
y^{(1)}_1 & y^{(2)}_1 & \cdots & y^{(n)}_1\\
y^{(1)}_2 & y^{(2)}_2 & \cdots & y^{(n)}_2\\
\vdots & \vdots & \ddots & \vdots\\
y^{(1)}_m & y^{(2)}_m & \cdots & y^{(n)}_m\\
\end{pmatrix}}
\begin{pmatrix}
b_1 \\
b_2 \\
\vdots\\
b_n
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_m
\end{pmatrix}\cdot
\begin{pmatrix}
y^{(1)}_1b_1 + y^{(2)}_1b_2 + \cdots + y^{(n)}_1b_n\\
y^{(1)}_2b_1 + y^{(2)}_2b_2 + \cdots + y^{(n)}_2b_n\\
\vdots \\
y^{(1)}_mb_1 + y^{(2)}_mb_2 + \cdots + y^{(n)}_mb_n\\
\end{pmatrix}\\
&=&
a_1\Bigl(y^{(1)}_1b_1 + y^{(2)}_1b_2 + \cdots + y^{(n)}_1b_n\Bigr)+
a_2\Bigl(y^{(1)}_2b_1 + y^{(2)}_2b_2 + \cdots + y^{(n)}_2b_n\Bigr)+
\cdots+
a_m\Bigl(y^{(1)}_mb_1 + y^{(2)}_mb_2 + \cdots + y^{(n)}_mb_n\Bigr)
\end{eqnarray}
となり、次に右辺は、
\begin{eqnarray}
(Y^T\boldsymbol{a})\cdot\boldsymbol{b} &=&
\left({
\begin{pmatrix}
y^{(1)}_1 & y^{(1)}_2 & \cdots & y^{(1)}_m\\
y^{(2)}_1 & y^{(2)}_2 & \cdots & y^{(2)}_m\\
\vdots & \vdots & \ddots & \vdots\\
y^{(n)}_1 & y^{(n)}_2 & \cdots & y^{(n)}_m\\
\end{pmatrix}}
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_m
\end{pmatrix}
\right)\cdot
\begin{pmatrix}
b_1 \\
b_2 \\
\vdots\\
b_n
\end{pmatrix}\\
&=&
\begin{pmatrix}
y^{(1)}_1a_1 + y^{(1)}_2a_2 + \cdots + y^{(1)}_ma_m\\
y^{(2)}_1a_1 + y^{(2)}_2a_2 + \cdots + y^{(2)}_ma_m\\
\vdots\\
y^{(n)}_1a_1 + y^{(n)}_2a_2 + \cdots + y^{(n)}_ma_m\\
\end{pmatrix}
\begin{pmatrix}
b_1 \\
b_2 \\
\vdots\\
b_n
\end{pmatrix}\\
&=&
b_1\Bigl(y^{(1)}_1a_1 + y^{(1)}_2a_2 + \cdots + y^{(1)}_ma_m\Bigr)+
b_2\Bigl(y^{(2)}_1a_1 + y^{(2)}_2a_2 + \cdots + y^{(2)}_ma_m\Bigr)+
\cdots+
b_n\Bigl(y^{(n)}_1a_1 + y^{(n)}_2a_2 + \cdots + y^{(n)}_ma_m\Bigr)\\
&=&
a_1\Bigl(y^{(1)}_1b_1 + y^{(2)}_1b_2 + \cdots + y^{(n)}_1b_n\Bigr)+
a_2\Bigl(y^{(1)}_2b_1 + y^{(2)}_2b_2 + \cdots + y^{(n)}_2b_n\Bigr)+
\cdots+
a_m\Bigl(y^{(1)}_mb_1 + y^{(2)}_mb_2 + \cdots + y^{(n)}_mb_n\Bigr)
\end{eqnarray}
となります。以上から、
\boldsymbol{a}\cdot(Y\boldsymbol{b})= (Y^T\boldsymbol{a})\cdot\boldsymbol{b}
が成り立つと言えます。
では、この内積の性質を踏まえて、性質2を証明していきましょう。
$X$ を対称行列とし、$X$ の異なる2つの固有値 $\lambda_i, \lambda_j$ を任意にとります。また、それぞれに対応する固有ベクトルを $\boldsymbol{u_i}, \boldsymbol{u_j}$ とします。
このとき、固有値、固有ベクトルの定義から以下のように表せます。
X\boldsymbol{u_i} = \lambda_i\boldsymbol{u_i}\\
X\boldsymbol{u_j} = \lambda_j\boldsymbol{u_j}
ここで、
\boldsymbol{u_i}\cdot (X\boldsymbol{u_j})
を考え、2通りで表してみます。先ほど証明した内積の性質と、$X$ が対称行列であることから
\begin{eqnarray}
\boldsymbol{u_i}\cdot (X\boldsymbol{u_j})&=&
(X^T\boldsymbol{u_i})\cdot\boldsymbol{u_j}\\
&=&
(X\boldsymbol{u_i})\cdot\boldsymbol{u_j}\\
&=&
(\lambda_i\boldsymbol{u_i})\cdot\boldsymbol{u_j}\\
&=&
\lambda_i\boldsymbol{u_i}\boldsymbol{u_j}
\end{eqnarray}
と表せ、
\begin{eqnarray}
\boldsymbol{u_i}\cdot(X\boldsymbol{u_j})
&=&
\boldsymbol{u_i}\cdot(\lambda_j\boldsymbol{u_j})\\
&=&
\lambda_j\boldsymbol{u_i}\boldsymbol{u_j}
\end{eqnarray}
のようにも表せます。勿論これら2つは等しい内積を表すので、
\lambda_i\boldsymbol{u_i}\boldsymbol{u_j}
=
\lambda_j\boldsymbol{u_i}\boldsymbol{u_j}\\
\leftrightarrows(\lambda_i-\lambda_j)\boldsymbol{u_i}\boldsymbol{u_j}=0
となります。ここで、$\lambda_i, \lambda_j$ は異なるので
\lambda_i-\lambda_j \neq 0
よって、
\boldsymbol{u_i}\boldsymbol{u_j}=0
$\boldsymbol{u_i}, \boldsymbol{u_j}$ の内積が $0$ なので、$\boldsymbol{u_i}$と $\boldsymbol{u_j}$ は直交すると言えます。
##スカラーをベクトルで微分
ここでは、スカラーをベクトルで微分するとどうなるかについて紹介します。
スカラーをベクトルで微分すると聞くと、「どういうこと？」と思いますが、定義に則って計算するだけなので簡単です。
###定義
\boldsymbol{a}=
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_n
\end{pmatrix}
とし、 $\mathcal{L}(\boldsymbol{a})$ を $\boldsymbol{a}$ の関数とする。このとき、$\mathcal{L}(\boldsymbol{a})$ の $\boldsymbol{a}$ による微分は次のように表される。
\frac{\partial\mathcal{L}(\boldsymbol{a})}{\partial \boldsymbol{a}} =
\begin{pmatrix}
\displaystyle\frac{\partial\mathcal{L}(\boldsymbol{a})}{\partial a_1} \\
\displaystyle\frac{\partial\mathcal{L}(\boldsymbol{a})}{\partial a_2} \\
\vdots\\
\displaystyle\frac{\partial\mathcal{L}(\boldsymbol{a})}{\partial a_n}
\end{pmatrix}
###性質
$X$ が $n$ 次正方行列のとき、
\frac{\partial}{\partial \boldsymbol{a}}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) = \bigl(X + X^T\bigr)\boldsymbol{a}
が成り立つ。
定義から、
\frac{\partial}{\partial \boldsymbol{a}}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) =
\begin{pmatrix}
\displaystyle\frac{\partial}{\partial a_1}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) \\
\displaystyle\frac{\partial}{\partial a_2}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) \\
\vdots\\
\displaystyle\frac{\partial}{\partial a_n}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr)
\end{pmatrix}
となります。ここで、
\begin{eqnarray}
\boldsymbol{a}^TX\boldsymbol{a}
&=&
\begin{pmatrix}
a_1&
a_2&
\cdots&
a_n
\end{pmatrix}
\begin{pmatrix}
x^{(1)}_1 & x^{(1)}_2 & \cdots & x^{(1)}_n\\
x^{(2)}_1 & x^{(2)}_2 & \cdots & x^{(2)}_n\\
\vdots & \vdots & \ddots & \vdots\\
x^{(n)}_1 & x^{(n)}_2 & \cdots & x^{(n)}_n\\
\end{pmatrix}
\begin{pmatrix}
a_1 \\
a_2 \\
\vdots\\
a_n
\end{pmatrix}\\
&=&
a_1\Bigl(x^{(1)}_1a_1 + x^{(2)}_1a_2 + \cdots + x^{(n)}_1a_n\Bigr)+
a_2\Bigl(x^{(1)}_2a_1 + x^{(2)}_2a_2 + \cdots + x^{(n)}_2a_n\Bigr)+
\cdots+
a_n\Bigl(x^{(1)}_na_1 + x^{(2)}_na_2 + \cdots + x^{(n)}_na_n\Bigr)
\end{eqnarray}
であるので、$\boldsymbol{a}$ の任意の成分 $a_k$ について、
\begin{eqnarray}
\frac{\partial}{\partial a_k}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr)
&=&
a_1x^{(k)}_1 +
a_2x^{(k)}_2 +
\cdots +
\bigl(a_1x^{(1)}_k + a_2x^{(2)}_k + \cdots + 2a_kx^{(k)}_k + \cdots + a_nx^{(n)}_k\bigr) +
\cdots +
a_nx^{(k)}_n\\
&=&
\bigl(x^{(k)}_1 + x^{(1)}_k\bigr)a_1 +
\bigl(x^{(k)}_2 + x^{(2)}_k\bigr)a_2 +
\cdots +
2x^{(k)}_ka_k +
\cdots +
\bigl(x^{(k)}_n + x^{(n)}_k\bigr)a_n
\end{eqnarray}
となります。したがって、
\begin{eqnarray}
\frac{\partial}{\partial \boldsymbol{a}}(\boldsymbol{a}^TX\boldsymbol{a}) &=&
\begin{pmatrix}
\displaystyle\frac{\partial}{\partial a_1}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) \\
\displaystyle\frac{\partial}{\partial a_2}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr) \\
\vdots\\
\displaystyle\frac{\partial}{\partial a_n}\bigl(\boldsymbol{a}^TX\boldsymbol{a}\bigr)
\end{pmatrix}\\
&=&
\begin{pmatrix}
2x^{(1)}_1a_1 +
\bigl(x^{(1)}_2 + x^{(2)}_1\bigr)a_2 +
\cdots +
\bigl(x^{(1)}_n + x^{(n)}_1\bigr)a_n\\
\bigl(x^{(2)}_1 + x^{(1)}_2\bigr)a_1 +
2x^{(2)}_2a_2 +
\cdots +
\bigl(x^{(2)}_n + x^{(n)}_2\bigr)a_n \\
\vdots\\
\bigl(x^{(n)}_1 + x^{(1)}_n\bigr)a_1 +
\bigl(x^{(n)}_2 + x^{(2)}_n\bigr)a_2 +
\cdots +
2x^{(n)}_na_n
\end{pmatrix}\\
&=&
\bigl(X + X^T\bigr)\boldsymbol{a}
\end{eqnarray}
が成り立つと言えます。
##行列の平均
次の分散共分散行列を紹介するためにも、ここで平均について紹介します。
\boldsymbol{y^{(k)}}=
\begin{pmatrix}
y^{(k)}_1 \\
y^{(k)}_2 \\
\vdots\\
y^{(k)}_m
\end{pmatrix},
\boldsymbol{y_k}=
\begin{pmatrix}
y^{(1)}_k \\
y^{(2)}_k \\
\vdots\\
y^{(n)}_k
\end{pmatrix}
とし、$\boldsymbol{y^{(k)}}$ の平均 $\overline{y^{(k)}}$ を次のように表すとします。
\begin{eqnarray}
\overline{y^{(k)}} &=&
\frac{1}{m}\bigl(y^{(k)}_1+y^{(k)}_2+\cdots+y^{(k)}_m\bigr)\\
&=&
\frac{1}{m}\sum\limits_{l=1}^my^{(k)}_l
\end{eqnarray}
このとき、$Y$ の平均 $\overline{\boldsymbol{Y}}$ を次のように定めます。
\begin{eqnarray}
\overline{\boldsymbol{Y}} &=&
\frac{1}{m}\sum\limits_{k=1}^my_k\\
&=&
\frac{1}{m}\bigl(y_1+y_2+\cdots+y_m\bigr)\\
&=&
\frac{1}{m}\hspace{3pt}
\left({
\begin{pmatrix}
y^{(1)}_1 \\ y^{(2)}_1 \\ \vdots \\ y^{(n)}_1
\end{pmatrix}}+
\begin{pmatrix}
y^{(1)}_2 \\ y^{(2)}_2 \\ \vdots \\ y^{(n)}_2
\end{pmatrix}+
\cdots +
\begin{pmatrix}
y^{(1)}_m \\ y^{(2)}_m \\ \vdots \\ y^{(n)}_m
\end{pmatrix}
\right)\\
&=&
\begin{pmatrix}
\displaystyle\frac{1}{m} \bigl(y^{(1)}_1 + y^{(1)}_2 + \cdots + y^{(1)}_m \bigr)\\
\displaystyle\frac{1}{m} \bigl(y^{(2)}_1 + y^{(2)}_2 + \cdots + y^{(2)}_m \bigr)\\
\vdots \\
\displaystyle\frac{1}{m} \bigl(y^{(n)}_1 + y^{(n)}_2 + \cdots + y^{(n)}_m \bigr)\\
\end{pmatrix}\\
&=&
\begin{pmatrix}
\overline{y^{(1)}}\\
\overline{y^{(2)}}\\
\vdots \\
\overline{y^{(n)}}
\end{pmatrix}
\end{eqnarray}
##分散共分散行列
ここでは、分散共分散行列の表し方について紹介します。
必要ないかもしれませんが、自分で手を動かして確認したので、一応紹介します。
V =
\frac{1}{m}\sum\limits_{k=1}^m\Bigl(y_k-\overline{\boldsymbol{Y}}\Bigr)\Bigl(y_k-\overline{\boldsymbol{Y}}\Bigr)^T
とすると、$V$ は $Y$ の分散共分散行列を表します。この確認をしてみましょう。
\begin{eqnarray}
V &=&
\displaystyle\frac{1}{m}\sum\limits_{k=1}^m\bigl(y_k-\overline{\boldsymbol{Y}}\bigr)\bigl(y_k-\overline{\boldsymbol{Y}}\bigr)^T\\
&=&
\displaystyle\frac{1}{m}\sum\limits_{k=1}^m
\begin{pmatrix}
y^{(1)}_k - \overline{y^{(1)}}\\
y^{(2)}_k - \overline{y^{(2)}}\\
\vdots \\
y^{(n)}_k - \overline{y^{(n)}}
\end{pmatrix}
\begin{pmatrix}
y^{(1)}_k - \overline{y^{(1)}}&
y^{(2)}_k - \overline{y^{(2)}}&
\cdots &
y^{(n)}_k - \overline{y^{(n)}}
\end{pmatrix}\\
&=&
\displaystyle\frac{1}{m}\sum\limits_{k=1}^m
\begin{pmatrix}
\bigl(y^{(1)}_k - \overline{y^{(1)}}\bigr)^2 &
\bigl(y^{(1)}_k - \overline{y^{(1)}}\bigr)\bigl(y^{(2)}_k-\overline{y^{(2)}}\bigr) &
\cdots &
\bigl(y^{(1)}_k - \overline{y^{(1)}}\bigr)\bigl(y^{(n)}_k-\overline{y^{(n)}}\bigr)\\
\bigl(y^{(2)}_k - \overline{y^{(2)}}\bigr)\bigl(y^{(1)}_k-\overline{y^{(1)}}\bigr) &
\bigl(y^{(2)}_k - \overline{y^{(2)}}\bigr)^2 &
\cdots &
\bigl(y^{(2)}_k - \overline{y^{(2)}}\bigr)\bigl(y^{(n)}_k-\overline{y^{(n)}}\bigr)\\
\vdots & \vdots & \ddots & \vdots \\
\bigl(y^{(n)}_k - \overline{y^{(n)}}\bigr)\bigl(y^{(1)}_k - \overline{y^{(1)}}\bigr) &
\bigl(y^{(n)}_k - \overline{y^{(n)}}\bigr)\bigl(y^{(2)}_k - \overline{y^{(2)}}\bigr) &
\cdots &
\bigl(y^{(n)}_k - \overline{y^{(n)}}\bigr)^2 \\
\end{pmatrix}\\
&=&
\begin{pmatrix}
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(1)}_1 - \overline{y^{(1)}}\biggr)^2 +
\biggl(y^{(1)}_2 - \overline{y^{(1)}}\biggr)^2 +
\cdots +
\biggl(y^{(1)}_m - \overline{y^{(1)}}\biggr)^2\Biggr) &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(1)}_1 - \overline{y^{(1)}}\biggr)\biggl(y^{(2)}_1 - \overline{y^{(2)}}\biggr) +
\biggl(y^{(1)}_2 - \overline{y^{(1)}}\biggr)\biggl(y^{(2)}_2 - \overline{y^{(2)}}\biggr) +
\cdots +
\biggl(y^{(1)}_m - \overline{y^{(1)}}\biggr)\biggl(y^{(2)}_m - \overline{y^{(2)}}\biggr)\Biggr) &
\cdots &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(1)}_1 - \overline{y^{(1)}}\biggr)\biggl(y^{(n)}_1 - \overline{y^{(n)}}\biggr) +
\biggl(y^{(1)}_2 - \overline{y^{(1)}}\biggr)\biggl(y^{(n)}_2 - \overline{y^{(n)}}\biggr) +
\cdots +
\biggl(y^{(1)}_m - \overline{y^{(1)}}\biggr)\biggl(y^{(n)}_m - \overline{y^{(n)}}\biggr)\Biggr)\\
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(2)}_1 - \overline{y^{(2)}}\biggr)\biggl(y^{(1)}_1 - \overline{y^{(1)}}\bigr) +
\biggl(y^{(2)}_2 - \overline{y^{(2)}}\biggr)\biggl(y^{(1)}_2 - \overline{y^{()}}\biggr) +
\cdots +
\biggl(y^{(2)}_m - \overline{y^{(2)}}\biggr)\biggl(y^{(1)}_m - \overline{y^{()}}\biggr)\Biggr) &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(2)}_1 - \overline{y^{(2)}}\biggr)^2 +
\biggl(y^{(2)}_2 - \overline{y^{(2)}}\biggr)^2 +
\cdots +
\biggl(y^{(2)}_m - \overline{y^{(2)}}\biggr)^2\Biggr) &
\cdots &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(2)}_1 - \overline{y^{(2)}}\biggr)\biggl(y^{(n)}_1 - \overline{y^{(n)}}\biggr) +
\biggl(y^{(2)}_2 - \overline{y^{(2)}}\biggr)\biggl(y^{(n)}_2 - \overline{y^{(n)}}\biggr) +
\cdots +
\biggl(y^{(2)}_m - \overline{y^{(2)}}\biggr)\biggl(y^{(n)}_m - \overline{y^{(n)}}\biggr)\Biggr)\\
\vdots & \vdots & \ddots & \vdots \\
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(n)}_1 - \overline{y^{(n)}}\biggr)\biggl(y^{(1)}_1 - \overline{y^{(1)}}\biggr) +
\biggl(y^{(n)}_2 - \overline{y^{(n)}}\biggr)\biggl(y^{(1)}_2 - \overline{y^{(1)}}\biggr) +
\cdots +
\biggl(y^{(n)}_m - \overline{y^{(n)}}\biggr)\biggl(y^{(1)}_m - \overline{y^{(1)}}\biggr)\Biggr) &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(n)}_1 - \overline{y^{(n)}}\biggr)\biggl(y^{(2)}_1 - \overline{y^{(2)}}\biggr) +
\biggl(y^{(n)}_2 - \overline{y^{(n)}}\biggr)\biggl(y^{(2)}_2 - \overline{y^{(2)}}\biggr) +
\cdots +
\biggl(y^{(n)}_m - \overline{y^{(n)}}\biggr)\biggl(y^{(2)}_m - \overline{y^{(2)}}\biggr)\Biggr) &
\cdots &
\displaystyle\frac{1}{m}\Biggl(
\biggl(y^{(n)}_1 - \overline{y^{(n)}}\biggr)^2 +
\biggl(y^{(n)}_2 - \overline{y^{(n)}}\biggr)^2 +
\cdots +
\biggl(y^{(n)}_m - \overline{y^{(n)}}\biggr)^2\Biggr)
\end{pmatrix}
\end{eqnarray}
となります。ここで、$V(y^{(i)})$ が $y^{(i)}$ の分散、$Cov(y^{(i)}, y^{(j)})$ が $y^{(i)}$ と $y^{(j)}$ の共分散を表すとし、
V^{(i,\hspace{2pt}j)} =
\begin{cases}
V(y^{(i)}) & (i = j)\\
Cov(y^{(i)}, y^{(j)}) & (i \neq j)
\end{cases}
とおけば、$V$ は以下のように表せます。
V =
\begin{pmatrix}
V^{(1,1)} & V^{(1,2)} & \cdots & V^{(1,n)}\\
V^{(2,1)} & V^{(2,2)} & \cdots & V^{(2,n)}\\
\vdots & \vdots & \ddots & \vdots\\
V^{(n,1)} & V^{(n,2)} & \cdots & V^{(n,n)}\\
\end{pmatrix}
これは、明らかに $Y$ の分散共分散行列を表しているので、
\frac{1}{n}\sum\limits_{k=1}^n\bigl(y_k-\overline{\boldsymbol{Y}}\bigr)\bigl(y_k-\overline{\boldsymbol{Y}}\bigr)^T
が $Y$ の分散共分散行列を表していると言えます。
##ラグランジュの未定乗数法
いよいよ事前準備もここで最後です。ただし、いちばん理解するのが難しいかもしれないです。
ここでは、ラグランジュの未定乗数法という最適化の手法を紹介します。
###目的
\boldsymbol{x}=
\begin{pmatrix}
x_1 \\
x_2 \\
\vdots\\
x_n
\end{pmatrix}
とする。このとき、制約条件 $g(\boldsymbol{x})=0$ の下で、目的関数 $f(\boldsymbol{x})$ を最大化または最小化する $\boldsymbol{x}$ を求める。
###解法
ある定数 $\lambda (\neq 0)$ を導入し、
\mathcal{L}(\boldsymbol{x}, \lambda) = f(\boldsymbol{x}) + \lambda g(\boldsymbol{x})
とおいたとき、
\begin{cases}
\displaystyle\frac{\partial}{\partial \boldsymbol{x}}\mathcal{L}(\boldsymbol{x}, \lambda) =0\hspace{1pt}\\
\displaystyle\frac{\partial}{\partial \lambda}\mathcal{L}(\boldsymbol{x}, \lambda) =0
\end{cases}
を解くことで、$f(\boldsymbol{x})$ を最大化、最小化させ得る $\boldsymbol{x}$ を求めることができる。
このように最適化問題を解くことをラグランジュの未定乗数法といい、 $\mathcal{L}(\boldsymbol{x}, \lambda)$ をラグランジュ関数、$\lambda$ をラグランジュ乗数といいます。
###解法が成り立つことの証明
では、なぜこの解法が成り立つかを証明していきますが、この証明をするためには必要な性質があるので、まずはその準備をします。
###勾配
必要な性質は勾配に関するものなので、ここでは勾配について紹介します。
####定義
関数 $f(\boldsymbol{x})$ の勾配 $\boldsymbol{\nabla}f$ は以下のように表される。
\boldsymbol{\nabla}f =
\begin{pmatrix}
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_1} \\
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_2} \\
\vdots\\
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_n}
\end{pmatrix}
####意味
この勾配がどのような意味を持つのかを考えていきましょう。
\boldsymbol{t}=
\begin{pmatrix}
x_1 \\
x_2 \\
\vdots\\
x_n
\end{pmatrix}
というベクトルを考え、そこから
\Delta \boldsymbol{t}=
\begin{pmatrix}
\Delta x_1 \\
\Delta x_2 \\
\vdots\\
\Delta x_n
\end{pmatrix}
だけ微小変化した場合を考えます。このとき、関数 $f(\boldsymbol{x})$ の変化 $\Delta f$ は以下のように表せます。
\begin{eqnarray}
\Delta f &=&
f(\boldsymbol{t} + \Delta\boldsymbol{t}) - f(\boldsymbol{t})\\
&=&
f
\left({
\begin{pmatrix}
x_1 + \Delta x_1\\
x_2 + \Delta x_2\\
\vdots\\
x_n + \Delta x_n
\end{pmatrix}
}\right)
-f
\left({
\begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}
}\right)
\end{eqnarray}
ここで、テイラー展開を用いて一次近似すれば、
f
\left({
\begin{pmatrix}
x_1 + \Delta x_1\\
x_2 + \Delta x_2\\
\vdots\\
x_n + \Delta x_n
\end{pmatrix}
}\right) \fallingdotseq
f
\left({
\begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}
}\right)+
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_1}\Delta x_1 +
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_2}\Delta x_2+
\cdots+
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_n}\Delta x_n
と表すことができるので、以下のように $\Delta f$ は $\boldsymbol{\nabla}f, \Delta\boldsymbol{t}$ の内積として表現できます。
\begin{eqnarray}
\Delta f &=&
f
\left({
\begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}
}\right)+
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_1}\Delta x_1 +
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_2}\Delta x_2+
\cdots+
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_n}\Delta x_n -
f
\left({
\begin{pmatrix}
x_1\\
x_2\\
\vdots\\
x_n
\end{pmatrix}
}\right)\\
&=&
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_1}\Delta x_1 +
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_2}\Delta x_2+
\cdots+
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_n}\Delta x_n\\
&=&
\begin{pmatrix}
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_1} &
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_2} &
\cdots&
\displaystyle\frac{\partial f(\boldsymbol{x})}{\partial x_n}
\end{pmatrix}
\begin{pmatrix}
\Delta x_1\\
\Delta x_2\\
\vdots\\
\Delta x_n
\end{pmatrix}\\
&=&
\boldsymbol{\nabla}f\Delta\boldsymbol{t}
\end{eqnarray}
さらに、$\boldsymbol{\nabla}f$ と $\Delta\boldsymbol{t}$ の成す角を $\theta$ とすれば、次のように $\Delta f$ を表すこともできます。
\Delta f =
\bigm\vert
\boldsymbol{\nabla}f
\bigm\vert
\bigm\vert
\Delta\boldsymbol{t}
\bigm\vert
\cos \theta
ここで、今の状況を簡単に図示してみると次のようになります。
したがって、$\Delta\boldsymbol{t}$ を一定とすれば、$\theta = 0$ のとき、$\Delta f$ が最大になります。$\theta = 0$ とは、$\Delta\boldsymbol{t}$ が $\boldsymbol{\nabla}f$ と同じ方向を向いているということです。つまり、$\boldsymbol{\nabla}f$ の方向に動いたとき $f(\boldsymbol{x})$ が最も大きく増加します。
以上から、勾配とは関数値が最も急激に増大する方向を表している、と言えます。
####勾配と等高線の性質
勾配は等高線と直交する。
ラクランジュの未定乗数法を理解する上で、この性質がかなり重要です。数式だけだと理解するのが難しいので、図を用いて幾何的に証明していきます。
3次元までしか図示できないので、ここでは3次元空間で考えていきます。
曲面 $z = f(x, y)$ と点 $(a, b)$ における接平面を考えます。図示すると次のようになります。
ここで、$z = f(a,b)$ における $f(x, y)$ と接平面の断面を考えると、$f(x, y)$ の断面は等高線に、接平面の断面は接線になります。これを図示すると次のようになります。
ここで、
\displaystyle\frac{\partial f(x,y)}{\partial x} = f_X,
\displaystyle\frac{\partial f(x,y)}{\partial y} = f_Y
とおけば、この接平面の方程式は以下のように表すことができます。
z = f(a,b) + f_X(a,b)(x-a) + f_Y(a,b)(y-b)
この接平面と $z = f(a,b)$ の交線が接線になるので、
f(a,b) = f(a,b) + f_X(a,b)(x-a) + f_Y(a,b)(y-b)\\
\leftrightarrows
f_X(a,b)(x-a) + f_Y(a,b)(y-b) = 0
が接線の方程式になります。この接線の方程式は以下のように変形することができます。
f_X(a,b)(x-a) + f_Y(a,b)(y-b) = 0\\
\leftrightarrows
\begin{pmatrix}
f_X(a,b)\\
f_Y(a,b)
\end{pmatrix}
\begin{pmatrix}
x-a\\
y-b
\end{pmatrix}=0\\
\leftrightarrows
\begin{pmatrix}
f_X(a,b)\\
f_Y(a,b)
\end{pmatrix}
\left({
\begin{pmatrix}
x\\
y
\end{pmatrix}}
-
\begin{pmatrix}
a\\
b
\end{pmatrix}
\right)=0
この
\begin{pmatrix}
x\\
y
\end{pmatrix}
-
\begin{pmatrix}
a\\
b
\end{pmatrix}
を図示すると次のようになります。$(x,y)$ は接線上の任意の点を表しているので、これは接線の傾きを表しています。
したがって、
\begin{pmatrix}
f_X(a,b)\\
f_Y(a,b)
\end{pmatrix}
\begin{pmatrix}
x-a\\
y-b
\end{pmatrix}=0
は、点 $(a,b)$ における勾配と等高線の接線の方向ベクトルの内積が $0$、すなわち、勾配と等高線が垂直に交わることを意味しています。
別解として、ここまで複雑に考えなくても、この
\Delta f =
\boldsymbol{\nabla}f\Delta\boldsymbol{t}
という式から、簡単に考えることもできます。
等高線というのは、関数値の値が変化しない点の集合なので、$\Delta f = 0$ になります。よって、
0 =
\boldsymbol{\nabla}f\Delta\boldsymbol{t}
となり、内積が $0$ であることから、$\boldsymbol{\nabla}f$ と $\Delta\boldsymbol{t}$ は直交することが分かります。等高線ということから、 $\Delta\boldsymbol{t}$ は関数値が変わらない方向へのベクトルを表しています。したがって、勾配と等高線は直交すると言うことができます。このように考えることもできます。
これで準備完了なので、なぜあの解法で最適解を見つけることができるか証明していきます。
ここでも図示したいので、3次元空間で考えていきます。
目的関数 $f(x,y)$ を制約条件 $g(x,y) = 0$ の下最大化する場合を考えます。
そこで、
\alpha>\beta>\gamma
を満たす $f(x,y)$ の等高線
\alpha = f(x,y), \beta = f(x,y), \gamma = f(x,y)
と $g(x,y) = 0$ を下に図示します。ただし、$\beta = f(x,y)$ と $g(x,y)=0$ は接するとします。
$f(x,y)$ を最大化したいので、等高線の中心の方へ行きたいですが、制約条件として $g(x,y)=0$ を満たさなければなりません。$g(x,y)=0$ を満たすとは、すなわち、$g(x,y)=0$ と $f(x,y)$ が交点を持つことです。等高線を3つ準備したので、それぞれについて考えてみます。
i) $\alpha = f(x,y)$
$f(x,y)$ の値としては大きいが、$g(x,y)=0$ を満たしません。$(g(x,y)=0$ と $f(x,y)$ が交点を持たない。$)$よって、最適解にはなり得ません。
ii) $\gamma = f(x,y)$
$f(x,y)$ が $g(x,y)=0$ と交点を持ちますが、$f(x,y)$ の値は最大ではありません。なぜなら、$g(x,y)=0$ と交点を持ち、$\gamma$ よりも大きい関数値をもつ等高線が存在するからです。
iii) $\beta = f(x,y)$
$f(x,y)$ が $g(x,y)=0$ と交点を持ち、制約条件の下 $\beta$ よりも大きい関数値は存在しません。したがって、最適解になり得ます。
以上から、$f(x,y)$ と $g(x,y) = 0$ が接する場合を考えれば良いです。これら2つの関数は接しているので、接線を共有しています。ここで、勾配と等高線の性質から、$f(x,y)$ の勾配 $\boldsymbol{\nabla}f$ と$g(x,y)$ の勾配 $\boldsymbol{\nabla}g$ は共通接線に対して垂直になります。したがって、$\boldsymbol{\nabla}f$ と $\boldsymbol{\nabla}g$ は平行になります。よって、ある定数 $\lambda$ を用いて
\boldsymbol{\nabla}f=\lambda\boldsymbol{\nabla}g
と表すことができます。$-\lambda$ を $\lambda$ と置き直せば、以下のように表せます。
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0
ここまでを整理すると、以下の条件を満たす $(x,y)$ を求めれば良いということになります。
\begin{cases}
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0\\
g(x,y) = 0
\end{cases}
ここで、$x,y,\lambda$ についての関数 $\mathcal{L}(x,y,\lambda)$ を
\mathcal{L}(x,y,\lambda) = f(x,y) + \lambda g(x,y)
のように定めると、
\begin{cases}
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0\\
g(x,y) = 0
\end{cases}
は
\begin{cases}
\displaystyle\frac{\partial}{\partial x}\mathcal{L}(x,y, \lambda) =0\hspace{1pt}\\
\displaystyle\frac{\partial}{\partial y}\mathcal{L}(x,y,\lambda) =0\hspace{1pt}\\
\displaystyle\frac{\partial}{\partial \lambda}\mathcal{L}(x,y,\lambda) =0
\end{cases}
と同値になります。実際、
\displaystyle\frac{\partial}{\partial x}\mathcal{L}(x,y,\lambda) = 0\\
\leftrightarrows
\displaystyle\frac{\partial}{\partial x}f(x,y)+
\lambda \displaystyle\frac{\partial}{\partial x}g(x,y)=0,\\
\displaystyle\frac{\partial}{\partial y}\mathcal{L}(x,y,\lambda) = 0\\
\leftrightarrows
\displaystyle\frac{\partial}{\partial y}f(x,y)+
\lambda \displaystyle\frac{\partial}{\partial y}g(x,y)=0
であるので、まとめれば、
\boldsymbol{\nabla}f + \lambda\boldsymbol{\nabla}g = 0
となり、また、
\displaystyle\frac{\partial}{\partial \lambda}\mathcal{L}(x,y,\lambda) = 0\\
\leftrightarrows
g(x,y)=0
であるので、同値であることが分かります。
今回は説明のために3次元で話しましたが、これらの性質は一般的にも成り立ちます。
とても長い事前準備でしたが、これで準備万端です。（準備が長すぎるので、後で事前準備の部分は違う記事にまとめてシャープにしようとかと思ってます。。）
#主成分分析
簡潔に言うと、主成分分析とは、データの情報をできるだけ保持したまま、高次元のデータを低次元へ圧縮することです。では、どのようにすれば、情報を保持したまま低次元へ圧縮できるか見ていきましょう。
##理論
理論を紹介します。
ここでは、
\boldsymbol{x} =
\begin{pmatrix}
\boldsymbol{x}_1 \\
\boldsymbol{x}_2 \\
\vdots\\
\boldsymbol{x}_m
\end{pmatrix}
とし、
\boldsymbol{x}_k=
\begin{pmatrix}
x^{(1)}_k \\
x^{(2)}_k \\
\vdots\\
x^{(n)}_k
\end{pmatrix},
\boldsymbol{x}^{(k)} =
\begin{pmatrix}
x^{(k)}_1 \\
x^{(k)}_2 \\
\vdots\\
x^{(k)}_m
\end{pmatrix}
とします。これだと分かりづらいので、特徴量が $n$ 次元で、データ点が $m$ 個あるようなデータをイメージすると分かりやすいと思います。
n 次元から1次元へ次元削減
まず、この $n$ 次元データを1次元データへ次元削減することを考えます。
$n$ 次元のデータとは、$n$ 個の軸から成る空間のデータです。言い換えれば、$n$ 個の直交する方向ベクトルで作られる空間にデータが存在しているということです。したがって、このデータをただ1つのベクトル
\boldsymbol{u}=
\begin{pmatrix}
u_1 \\
u_2 \\
\vdots\\
u_m
\end{pmatrix}
へ射影することにより、1次元へ次元削減をすることができます。
これで、次元削減の方法は分かりました。では、『データの情報をできるだけ保持する』ためにはどうすれば良いでしょうか。
結論から言うと、データの分散ができるだけ大きくなるような方向ベクトルへ射影すれば良いです。ここは、直感的に理解するものなのかもしれません。例えば、$\boldsymbol{u}$ 軸へ射影した結果が以下の2通りだったとします。
このとき、どちらの射影の方がデータの情報を持っていそうでしょうか。
元々のデータの散らばり具合ができるだけ失われていない方がデータの情報を持っていそうですよね。言い換えると、射影後のデータの分散が大きい方が、元データの情報を持っていると考えられる訳です。
よって、射影後のデータ $\boldsymbol{u}^T\boldsymbol{x}$ の分散が最大になるような $\boldsymbol{u}$ を見つけ、その $\boldsymbol{u}$ による射影で次元削減を行えば良さそうです。
では、射影後の分散がどのように表されるか見ていきましょう。
まず、$\boldsymbol{x}$ の平均を $\overline{\boldsymbol{x}}$、$\boldsymbol{u}$ による射影後の平均を $\overline{\boldsymbol{x}}_{\boldsymbol{u}}$ とすれば、
\begin{eqnarray}
\overline{\boldsymbol{x}}_{\boldsymbol{u}}&=&
\displaystyle\frac{1}{m}\sum^m_{k=1}\boldsymbol{u}^T\boldsymbol{x}_k\\
&=&
\displaystyle\frac{1}{m}\Bigl(
\bigl(u_1x^{(1)}_1 + \cdots + u_mx^{(n)}_1\bigr) + \cdots + \bigl(u_1x^{(1)}_m + \cdots + u_mx^{(n)}_m\bigr)\Bigr)\\
&=&
u_1\biggl(\displaystyle\frac{1}{m}\bigl(x^{(1)}_1 + \cdots + x^{(1)}_m\bigr)\biggr) +
\cdots +
u_m\biggl(\displaystyle\frac{1}{m}\bigl(x^{(n)}_1 + \cdots + x^{(n)}_m\bigr)\biggr)\\
&=&
u_1\overline{x^{(1)}} + \cdots + u_m\overline{x^{(m)}}\\
&=&
\boldsymbol{u}^T\overline{\boldsymbol{x}}
\end{eqnarray}
と表すことができます。次に、$\boldsymbol{x}$ の分散共分散行列を $V$、$\boldsymbol{u}$ による射影後の分散を $V_\boldsymbol{u}$ とすれば、
\begin{eqnarray}
V_\boldsymbol{u} &=&
\displaystyle\frac{1}{m}\sum^m_{k=1}\bigl(\boldsymbol{u}^Tx_k - \boldsymbol{u}^T\overline{x}\bigr)^2\\
&=&
\displaystyle\frac{1}{m}\sum^m_{k=1}\Biggl(u_1\biggl(x^{(1)}_k-\overline{x^{(1)}}\biggr)+\cdots+u_m\biggl(x^{(n)}_k-\overline{x^{(n)}}\biggr)\Biggr)^2\\
&=&
\displaystyle\frac{1}{m}\sum^m_{k=1}\Biggl(u_1^2\biggl(x^{(1)}_k-\overline{x^{(1)}}\biggr)^2+
\cdots+
u_m^2\biggl(x^{(n)}_k-\overline{x^{(n)}}\biggr)^2+
2u_1u_2\biggl(x^{(1)}_k-\overline{x^{(1)}}\biggr)\biggl(x^{(2)}_k-\overline{x^{(2)}}\biggr)+
\cdots+
2u_1u_m\biggl(x^{(1)}_k-\overline{x^{(1)}}\biggr)\biggl(x^{(n)}_k-\overline{x^{(n)}}\biggr)+
2u_2u_3\biggl(x^{(2)}_k-\overline{x^{(2)}}\biggr)\biggl(x^{(3)}_k-\overline{x^{(3)}}\biggr)+
\cdots+
2u_2u_m\biggl(x^{(2)}_k-\overline{x^{(2)}}\biggr)\biggl(x^{(n)}_k-\overline{x^{(n)}}\biggr)+
\cdots+
2u_{m-1}u_m\biggl(x^{(n-1)}_k-\overline{x^{n-1}}\biggr)\biggl(x^{(n)}_m-\overline{x^{(n)}}\biggr)\Biggr)\\
&=&
u_1^2\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(1)}_1-\overline{x^{(1)}}\biggr)^2+
\cdots+
\biggl(x^{(1)}_m-\overline{x^{(1)}}\biggr)^2\Biggr)\right\}+
\cdots+
u_m^2\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(n)}_1-\overline{x^{(n)}}\biggr)^2+
\cdots+
\biggl(x^{(n)}_m-\overline{x^{(n)}}\biggr)^2\Biggr)\right\}+
2u_1u_2\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(1)}_1-\overline{x^{(1)}}\biggr)\biggl(x^{(2)}_1-\overline{x^{(2)}}\biggr)+
\cdots+
\biggl(x^{(1)}_m-\overline{x^{(1)}}\biggr)\biggl(x^{(2)}_m-\overline{x^{(2)}}\biggr)\Biggr)\right\}+
\cdots+
2u_1u_m\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(1)}_1-\overline{x^{(1)}}\biggr)\biggl(x^{(n)}_1-\overline{x^{(2)}}\biggr)+
\cdots+
\biggl(x^{(1)}_m-\overline{x^{(1)}}\biggr)\biggl(x^{(n)}_m-\overline{x^{(n)}}\biggr)\Biggr)\right\}+
2u_2u_3\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(2)}_1-\overline{x^{(2)}}\biggr)\biggl(x^{(3)}_1-\overline{x^{(3)}}\biggr)+
\cdots+
\biggl(x^{(2)}_m-\overline{x^{(2)}}\biggr)\biggl(x^{(3)}_m-\overline{x^{(3)}}\biggr)\Biggr)\right\}+
\cdots+
2u_2u_m\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(2)}_1-\overline{x^{(2)}}\biggr)\biggl(x^{(n)}_1-\overline{x^{(n)}}\biggr)+
\cdots+
\biggl(x^{(2)}_m-\overline{x^{(2)}}\biggr)\biggl(x^{(n)}_m-\overline{x^{(n)}}\biggr)\Biggr)\right\}+
\cdots+
2u_{m-1}u_m\left\{\displaystyle\frac{1}{m}\Biggl(\biggl(x^{(n-1)}_1-\overline{x^{(n-1)}}\biggr)\biggl(x^{(n)}_1-\overline{x^{(n)}}\biggr)+
\cdots+
\biggl(x^{(n-1)}_m-\overline{x^{(n^1)}}\biggr)\biggl(x^{(n)}_m-\overline{x^{(n)}}\biggr)\Biggr)\right\}\\
&=&
u_1^2V^{(1,1)} + \cdots + u_m^2V^{(n,n)} +
2u_1u_2V^{(1,2)} + \cdots + 2u_1u_mV^{(1,n)} +
2u_2u_3V^{(2,3)} + \cdots + 2u_2u_mV^{(2,n)} +
\cdots + 2u_{m-1}u_mV^{(n-1,n)}
\end{eqnarray}
と表せます。ここで、
Cov(x^i, x^j) = Cov(x^j, x^i)
であるので、
2u_iu_jV^{(i,\hspace{2pt}j)} = u_iu_jV^{(i,\hspace{2pt}j)} + u_ju_iV^{(j,\hspace{2pt}i)}
のように分解できます。よって、
\begin{eqnarray}
V_\boldsymbol{u} &=& \boldsymbol{u}^T
\begin{pmatrix}
V^{(1,1)} & V^{(1,2)} & \cdots & V^{(1,n)}\\
V^{(2,1)} & V^{(2,2)} & \cdots & V^{(2,n)}\\
\vdots & \vdots & \ddots & \vdots\\
V^{(n,1)} & V^{(n,2)} & \cdots & V^{(n,n)}\\
\end{pmatrix}
\boldsymbol{u}\\
&=&
\boldsymbol{u}^TV\boldsymbol{u}
\end{eqnarray}
と表すことができます。
この式から明らかなように、$\boldsymbol{u}$ をスカラー倍すれば、際限なく $V_\boldsymbol{u}$ を大きくすることができます。例えば、
\boldsymbol{u}' = 2\boldsymbol{u}
として、$\boldsymbol{u}'$ による射影後の分散を $V_\boldsymbol{u}'$ とすれば、
V_\boldsymbol{u}' = 4V_\boldsymbol{u}
となり、分散が4倍になります。これでは、議論ができなくなってしまうので、
\vert\vert\boldsymbol{u}\vert\vert = 1\\
\leftrightarrows \boldsymbol{u}^T\boldsymbol{u} = 1
という条件を設定します。
ここまでをまとめると、
1次元への次元削減とは、$\boldsymbol{u}^T\boldsymbol{u} = 1$ の下で、$\boldsymbol{u}^TV\boldsymbol{u}$ を最大にする $\boldsymbol{u}$ を見つけ、その $\boldsymbol{u}$ で射影することであると言えます。
これは条件付きの最大化問題です。もうお気付きですね？ラグランジュの未定乗数法を使って最適解を求めることができます。
f(u) = \boldsymbol{u}^TV\boldsymbol{u},
\\g(u) = 1-\boldsymbol{u}^T\boldsymbol{u}
とおき、$\boldsymbol{u}$ とある定数 $\lambda$ についての関数 $\mathcal{L}(\boldsymbol{u}, \lambda)$ を以下のように定めます。
\mathcal{L}(\boldsymbol{u},\lambda) = f(\boldsymbol{u})+\lambda g(\boldsymbol{u})
ラグランジュの未定乗数法より、
\begin{cases}
\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\mathcal{L}(\boldsymbol{u}, \lambda) =0\hspace{1pt}\\
\displaystyle\frac{\partial}{\partial \lambda}\mathcal{L}(\boldsymbol{u}, \lambda) =0
\end{cases}
を解けば良いと分かります。
\begin{cases}
\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\mathcal{L}(\boldsymbol{u}, \lambda) =0\hspace{1pt}\\
\displaystyle\frac{\partial}{\partial \lambda}\mathcal{L}(\boldsymbol{u}, \lambda) =0
\end{cases}\\
\leftrightarrows
\begin{cases}
\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\boldsymbol{u}^TV\boldsymbol{u}+
\lambda\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\bigl(1-\boldsymbol{u}^T\boldsymbol{u}\bigr)
= 0\hspace{1pt}\\
1-\boldsymbol{u}^T\boldsymbol{u}=0
\end{cases}
スカラーをベクトルで微分すれば良いので、事前準備で証明した性質と定義から、
\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\boldsymbol{u}^TV\boldsymbol{u}+
\lambda\displaystyle\frac{\partial}{\partial \boldsymbol{u}}\bigl(1-\boldsymbol{u}^T\boldsymbol{u}\bigr)
= 0\\
\leftrightarrows
\bigl(V+V^T\bigr)\boldsymbol{u}+\lambda\cdot(-2\boldsymbol{u})=0
となります。さらに、$V$ は対称行列なので、
V = V^T
が成り立ちます。よって、
2V\boldsymbol{u}-2\lambda\boldsymbol{u}=0
となります。移項して2で割れば、
V\boldsymbol{u}=\lambda\boldsymbol{u}
という形になります。したがって、解くべき方程式は次のようになります。
\begin{cases}
V\boldsymbol{u}=\lambda\boldsymbol{u}\\
\boldsymbol{u}^T\boldsymbol{u}=1
\end{cases}
さて、この連立方程式の1つ目の式
V\boldsymbol{u}=\lambda\boldsymbol{u}
に見覚えがあるのではないでしょうか？固有値と固有ベクトルの定義式と同じになっていますよね？
つまり、求めたい $\boldsymbol{u}$ とは、$\boldsymbol{x}$ の分散共分散行列 $V$ の固有ベクトルであり、$\lambda$ は $V$ の固有値である、ということが分かります。
さらに、この式の両辺に左から $\boldsymbol{u}^T$ をかけると、解くべき連立方程式は、
\begin{cases}
\boldsymbol{u}^TV\boldsymbol{u} = \boldsymbol{u}^T\lambda\boldsymbol{u}\\
\boldsymbol{u}^T\boldsymbol{u}=1
\end{cases}
となるので、$\lambda$ について解けば、
\boldsymbol{u}^TV\boldsymbol{u} =\lambda
と表せます。つまり、$\lambda$ は $V$ の固有値であるのと同時に、射影後の分散をも表している、と分かります。
今回この射影後の分散を最大化したかったので、最も大きい $\lambda$ が射影後の分散を最大にすると分かります。すなわち、射影後の分散は、最大固有値のとき最大になり、その固有値に対応する固有ベクトルが今回求めたかった $\boldsymbol{u}$ になります。また、この $\boldsymbol{u}$ 、すなわち、最大固有値に対応する固有ベクトルを第1主成分といいます。
###n 次元から d 次元へ次元削減
1次元への次元削減の方法は分かったので、$d(\leq n)$ 次元への次元削減の方法を見てみましょう。
主成分分析の考え方は、分散が最大になるようなベクトルを見つけることでした。もし2次元へ次元削減するなら、射影後の分散を最も大きくするベクトルとその次に大きくするベクトルの2つを見つけてくれば良い訳です。3次元なら3個、4次元なら4個、$d$ 次元なら $d$ 個見つけてくれば良い訳です。
ここで、固有値と固有ベクトルの性質1から、$V$ は $n$ 個の固有値を持つので、それらの固有値を大きい順に、$\lambda_1,\lambda_2,\cdots,\lambda_n$ とします。式で表せば、以下のようになります。
\lambda_1>\lambda_2>\cdots>\lambda_n
また、それぞれの固有値に対応する固有ベクトルをそれぞれ $\boldsymbol{u}_1,\boldsymbol{u}_2,\cdots,\boldsymbol{u}_n$ とします。
改めて、次の式を見てみましょう。
\boldsymbol{u}^TV\boldsymbol{u} =\lambda
この式から明らかに、射影後の分散を2番目に大きくするのは $\lambda_2$ になります。$d$ 番目に大きくするのは $\lambda_d$ になります。すなわち、$d$ 次元へ圧縮するのであれば、固有値を大きいものから $d$ 個探し、それぞれに対応する固有ベクトルで射影させれば良い訳です。この2番目に大きい固有値に対応する固有ベクトル $\boldsymbol{u}_2$ を第2主成分、3番目に大きい固有値に対応する固有ベクトル $\boldsymbol{u}_3$ を第3主成分、$d$ 番目に大きい固有値に対応する固有ベクトル $\boldsymbol{u}_d$ を第 $d$ 主成分といいます。
ちなみに、固有値と固有ベクトルの性質2から、異なる固有値に対応する固有ベクトル同士は直交するので、$\boldsymbol{u}_1,\boldsymbol{u}_2,\cdots,\boldsymbol{u}_n$ はいずれも直交します。3次元空間なら、$x$ 軸、$y$ 軸、$z$ 軸のようなイメージです。（厳密には、$\boldsymbol{u}_k$ は大きさが1なので、軸ではないですが。）
これが主成分分析の理論です。お疲れ様でした。
#理論を理解する中で湧いた疑問
最後に、生じた疑問について共有しようかと思います。
##固有方程式が重解を持つ場合不都合はないのか
一般的に、$n$ 次正方行列の固有方程式が重解を持った場合、固有ベクトルが $n$ 個見つからない場合があります。固有ベクトルが少ないと、思っている次元数までデータの次元を落とせない訳です。具体的に言うと、10次元のデータを3次元へ落としたいが、固有ベクトルが2個しか見つからず、2次元までしか次元削減ができないという可能性がある訳です。
これについては問題ありませんでした。主成分分析で扱う $n$ 次正方行列は分散共分散行列なので、$n$ 次実対称行列です。実対称行列に関しては重解を持っていても必ず $n$ 個の固有ベクトルが見つかるので、固有ベクトルが足らないということはありません。
しかし、そうなった場合また1つ問題が生じます。事前準備の固有値と固有ベクトルの性質2で紹介したように、対称行列の異なる固有値に対応する固有ベクトル同士は直交しますが、同じ固有値に対応する異なる固有ベクトル同士は一般に直交しません。こうなるとあまり良くないことが起きます。今回紹介しなかったのですが、大きさ1の固有ベクトルは正規直交系になるので、直交しない固有ベクトルがあると、正規直交基底を成し得ない訳です。すると、データがしっかり射影されない場合が発生してしまうのです。
では、重解を持つようなデータでは主成分分析を行えないのでしょうか？
私は、どのようなデータでも主成分分析を行って大丈夫だと考えます。なぜなら、そもそも実際のデータが全く同じ固有値を持つことは考えづらいからです。例えば、自分と全く同じ身長の人は存在するでしょうか。「いや、普通にいますよ！」と思われるかもしれませんが、それは少数第1位や少数第2位で四捨五入しているからです。少数第100位まで全く同じ人はそうはいないでしょう。（そこまで精密に測定できる機器があるかは分かりませんが。）つまり、実際に観測したデータが重解を持つことは、ほぼほぼあり得ないので、考慮する必要はないと考えられます。
##射影前より射影後の分散が大きくなることはあるのか
主成分分析のアイデアはできるだけデータの情報を保持するために分散を最大化するというものでした。では、射影した結果、分散が射影前よりも大きくなってしまっても情報を保持できていると言えるのでしょうか。そもそも射影前のデータは分散共分散行列で、射影後の分散はスカラーな分散なので、比較できないのではないかという気もします。。。しかし、もし比較する方法があり、射影によって分散が射影前より大きくなる可能性があるのであれば、
\left|\boldsymbol{u}^TV\boldsymbol{u}-V\right|
が小さくなるような $\boldsymbol{u}$ を選択した方が情報を保持できていそうに思えます。
この疑問に関しては答えが出ていません。。しかし、この主成分分析が実際に利用されていることを考えれば、きっとこういった心配は不要なのでしょう。
#おわりに
かなり長くなってしまいましたが、最初から読んで頂くだけで主成分分析の理論が理解できる内容になっていると思います。読んでくださった方の理解の一助になれば嬉しいです。
#参考文献
[1] 後藤正幸・小林学 「入門 パターン認識と機械学習」
[2] 金谷健一 「これなら分かる最適化数学-基礎原理から計算手法まで-」