この記事はこの記事は新歓ブログリレー2024、52日目の記事です。

本記事では1, 2年生に向けて機械学習に関連する線形代数のトピック，特に主成分分析に関連したトピックを紹介します。理論的なところから始め，最後にPythonで簡単な実装をしてみます。「線形代数学第二」(3Q or 4Qで履修) 程度の内容で1年生の線形代数のモチベーションになることを目標にして書きました。履修した2年生は難なく読めると思います。1年生は全部は読めないかもしれませんがぜひ眺めてみてください。線形代数学第一で僕を担当してくださった先生は以下で扱う線形代数のトピックを雑談で話されていて，その時に知っていればより面白かったのではとも思います。

主成分分析は流行りのいわゆる人工知能と直接的な関わりはありません。しかし応用として，例えば単語ベクトルの可視化やオートエンコーダーなど最近のいわゆる人工知能の技術が挙げられます。単語は高次元のベクトルとして扱うことができて，意味の近さがベクトルとしての近さに反映されます。これを可視化したのが下図です。（位置ベクトルとして点で表されています。）2次元まで削減しているにも関わらず，関係性を保ったままほとんど重なることなく図示できています。関係性の例として，右下の日本を表す点はオーストラリアやアメリカ、カナダといった文化的に近い国を表す点の近くにありますね。

(上図は[4]からの引用でt-SNEという手法を用いています。)
それでは実際に主成分分析の紹介に入っていきます。

射影行列

まず主成分分析の理論を説明するために便利な射影行列について解説します。

単位ベクトルの復習

単位ベクトルはとても便利な道具で，

• 内積を取ることで単位ベクトルの方向成分の符号付きの長さを求められる。
  • ベクトル(3, 5)と$x$軸の単位ベクトル(1, 0)とで内積をとると$x$軸方向成分である3が得られます。ただし負になることもあるので「符号付き」としました。
• $k\mathbf{e}$で$\mathbf{e}$方向，符号付き長さ$k$のベクトルを作れる。
  これらを用いると，正射影ベクトルがすぐ導けます。$\mathbf{x}$の$\mathbf{e}$方向の正射影ベクトル${\mathbf{x}}_{//}$は

$${\mathbf{x}}_{//}=(\mathbf{x},\mathbf{e})\mathbf{e}$$

と書けます。ただし，$(\cdot ,\cdot )$は内積とします。

射影する関数は線形性を持ちそうです。よって行列として射影する関数を表せそうです。
線形性を持ちそうと言うのは，図形的に，ベクトルを定数倍してから射影しても射影してから定数倍しても等しそうで，和の射影と射影の和も等しそうということです。関数が線形性を持つとき行列で表すことができます。
実際，

$$(\mathbf{x},\mathbf{e})\mathbf{e}=\mathbf{e}(\mathbf{e},\mathbf{x})=\mathbf{e}{\mathbf{e}}^{\top }\mathbf{x}(\text{ここで}\top \text{は転置記号で}{A}^{\top }\text{で}A\text{の転置を表す})$$

となるので$P=\mathbf{e}{\mathbf{e}}^T$と置くと，${\mathbf{x}}_{//}=P\mathbf{x}$と書けます。 この$P$は直交射影行列と呼ばれます。ここでは直交でない射影行列は扱わないので単に射影行列と書くことにします。ところで$\mathbf{e}$方向に，射影したものをもう一度射影しても変わりません。また式の形からして対称行列です。すなわち次が成り立ちます。

🐾 直交射影行列の性質

$$P^2=P$$

$$P^{\top }=P$$

補足：こちらの2つを満たす行列を直交射影行列と定義する場合が多いです。本記事では直感的なわかりやすさのために逆側から紹介しています。

問題
${\mathbf{e}}_1=(1,0,0{)}^{\top },{\mathbf{e}}_2=(0,1,0{)}^{\top },{\mathbf{e}}_3=(0,0,1{)}^{\top }$とする。

1. ${\mathbf{e}}_1$の張る直線に射影する射影行列を求めてください。
   (直線は向きだけではなく通る点も決めないと定まらないですが，今回は射影を考えるだけなので向きさえ決まれば問題ないです。)

2. ${\mathbf{e}}_1,{\mathbf{e}}_2$の張る平面に射影する射影行列を求めてください。
   ヒント: ベクトル $\mathbf{x}$の平面への正射影はどのように書けるでしょうか。

3. ${\mathbf{e}}_1,{\mathbf{e}}_2,{\mathbf{e}}_3$は${\mathbb{R}}^3$の正規直交基底です。${\mathbf{e}}_1,{\mathbf{e}}_2$の張る部分空間$\mathcal{U}$の直交補空間${\mathcal{U}}^{\perp }$への射影行列$P_{{\mathcal{U}}^{\perp }}$を求めてください。

主成分分析

主成分分析は機械学習の手法の一つで，簡単に言うと情報をなるべく減らさずにデータの次元を下げる手法です。

データのパラメータが大きいときに，それを可視化したり扱ったりするのは大変になりがちです。まず，少なくとも2~3次元まで落とさないと可視化は難しそうです。しかも，高次元での振る舞いは奇妙で1,000,000次元の単位超立方体で無作為に2つの点を選ぶとその平均距離は約408.25(3次元は約0.66)。 データの点同士がより離れているので低次元のときと比べて予想の信頼性が落ちてしまいます。
実際にどのように次元を落とすかですが，下の画像のデータは直線を当てはめて射影することで1次元に落とせそうです。

このように直線や平面の当てはめて射影して次元を下げるのが主成分分析です。

どのように当てはめるか　分散は情報の大きさ

さてどのような指標を使って最もフィットする直線を決めれば良いでしょうか？結論としては主成分分析では射影行列を使って，ある直線に射影し，その射影と元のベクトルとの差の長さ($\Vert \mathbf{x}-P_{\mathcal{U}}\mathbf{x}\Vert =\Vert P_{{\mathcal{U}}^{\perp }}\mathbf{x}\Vert$)の2乗をデータ全体で和をとったもの$E$を最小化するような直線を求めます。

それではベクトル$\mathbf{x}$を直線にあてはめるとき，なぜこのような指標で考えるのでしょうか？
それは分散が最も大きく維持される当てはめを見つけることができるからです。分散というのは情報の大きさと見ることができます。例えば，全人類というデータがあったときに染色体の本数という一次元の指標に当てはめると，確かに次元は下がりましたが，ほとんどの人類が46本のところに集中して(つまり分散が小さい)まるっきり情報がありません。ほとんどの点が46付近に集中するため，射影した後のデータを見てグループ分けなんかもできないです。また，最頻値や平均値であろう46という数字がありさえすればよく，あえて分布として考える必要もありません。一方，身長や体重，年齢などそれなりにばらつき，分散が大きいので分布が意味を持ちます。それではどのようにあてはめる直線を求めるのか見ていきましょう。

ここで，補足として，当てはめる直線や平面と述べてきましたが，正確にはアフィン空間の当てはめです。アフィン空間とは部分空間を並行移動したものです。例えば，原点を通らない直線も並行移動して原点を通るようにすると部分空間になるのでアフィン空間です。
アフィン空間を求めるとはつまり，通る1点と基底を求めるということで，直線の当てはめであれば、通る1点と方向ベクトルです。
まず，求めるアフィン空間$\mathcal{A}$はデータの重心を通ることを示します。すなわち($\mathcal{A}$の基底$\mathbf{e}$を固定したとき)$\mathbf{c}={\mathbf{x}}_g$で$E$は極小値を取ることを示します。
$\mathbf{e}$で張られる空間を$\mathcal{U}$とします。
射影行列の性質を使って

$$\begin{array}{rl}{\nabla }_{\mathbf{c}}E& ={\nabla }_{\mathbf{c}}\left(\sum _n^N\Vert P_{{\mathcal{U}}^{\perp }}({\mathbf{x}}_n-\mathbf{c}){\Vert }^2\right)\\ & ={\nabla }_{\mathbf{c}}\left(\sum _n^N({\mathbf{x}}_n-\mathbf{c}{)}^{\top }{P}_{{\mathcal{U}}^{\perp }}({\mathbf{x}}_n-\mathbf{c})\right)\\ & =-2P_{{\mathcal{U}}^{\perp }}\sum _n^N({\mathbf{x}}_n-\mathbf{c})\end{array}$$

${\nabla }_{\mathbf{c}}E=0$とおくと，$P_{{\mathcal{U}}^{\perp }}{\sum }_n({\mathbf{x}}_n-\mathbf{c})=0$
ところで一般にベクトル空間$\mathcal{U}$とその直交補空間${\mathcal{U}}^{\perp }$に対して

$$\ker P_{\mathcal{U}}=\mathrm{im}{P}_{{\mathcal{U}}^{\perp }}$$

であることが示せます。(やってみよう)
これらより，${\nabla }_{\mathbf{c}}E=0$を満たす${\sum }_n({\mathbf{x}}_n-\mathbf{c})$の値は$\mathcal{U}$全体で，$\mathcal{U}$のどの点をとっても同じですが($\mathcal{A},\mathcal{U}$の基底が等しいため，例えばこの例では二つの平行な直線であるため)

例えば$0$を取ると$\mathbf{c}=\frac{1}{N}\sum _n^N{\mathbf{x}}_n={\mathbf{x}}_g$よって$\mathcal{A}$は重心を通ることが示されました。
どの点を取っても同じということについて：上の例だと$k\mathbf{e}\in \mathcal{U}$という点を任意に取ると，$\mathbf{c}=\frac{1}{N}\sum _n^N{\mathbf{x}}_n-\frac{k}{N}\mathbf{e}$

このときアフィン空間の定義から$\mathbf{c}+\frac{k}{N}\mathbf{e}\in \mathcal{A}$すなわち${\mathbf{x}}_g\in \mathcal{A}$となって結局重心は求めるアフィン空間に属するということです。


さて今度は基底$\mathbf{e}$を求めます。
まず，${\mathbf{x}}_i$に${\mathbf{x}}_i-{\mathbf{x}}_g$を代入して中心化しておきます。これで当てはめるアフィン空間は原点を通ります。

三平方の定理より

$$\Vert \mathbf{x}{\Vert }^2=\Vert P_{{\mathcal{U}}^{\perp }}\mathbf{x}{\Vert }^2+\Vert P_{\mathcal{U}}\mathbf{x}{\Vert }^2$$

三平方の定理という図形的考察に頼らなくとも前述の射影行列の性質を使うことでも求められます。(やってみよう)

一旦，分散のことは置いておいて，単に当てはめた直線と点の距離の2乗和の総和を小さくすることを考えてみます。これは化学実験や物理学実験で使う最小二乗法に似ていますが，最小二乗法は残差(予測値とデータの値の差)の2乗であって，当てはめた直線と点の距離の2乗ではないです。

$$\begin{array}{rl}\sum _n\Vert P_{{\mathcal{U}}^{\perp }}{\mathbf{x}}_n{\Vert }^2& =\sum _n\Vert {\mathbf{x}}_n{\Vert }^2-\sum _n\Vert P_{\mathcal{U}}{\mathbf{x}}_n{\Vert }^2\\ & =\sum _n\Vert {\mathbf{x}}_n{\Vert }^2-\sum _n\Vert (\mathbf{e},{\mathbf{x}}_n){\Vert }^2\\ & =\sum _n\Vert {\mathbf{x}}_n{\Vert }^2-\sum _n{\mathbf{e}}^{\top }{\mathbf{x}}_n{\mathbf{x}}_n^{\top }\mathbf{e}\\ & =\sum _n\Vert {\mathbf{x}}_n{\Vert }^2-\sum _n{\mathbf{e}}^{\top }{\mathbf{x}}_n{\mathbf{x}}_n^{\top }\mathbf{e}\\ & =\sum _n\Vert {\mathbf{x}}_n{\Vert }^2-{\mathbf{e}}^{\top }X{X}^{\top }\mathbf{e}\end{array}$$

${\sum }_n\Vert {\mathbf{x}}_n{\Vert }^2$は定数であるので，${\sum }_n\Vert P_{{\mathcal{U}}^{\perp }}{\mathbf{x}}_n{\Vert }^2$の最小化問題は${\mathbf{e}}^{\top }X{X}^{\top }\mathbf{e}$の最大化問題に帰着されました。これについて次々節で詳しく見ていきます。

統計的な意味

${\mathbf{x}}_i$は${\mathbf{x}}_i-{\mathbf{x}}_g$だったことを思い出すと，最大化した${\mathbf{e}}^{\top }X{X}^{\top }\mathbf{e}={\sum }_n\Vert P_{\mathcal{U}}{\mathbf{x}}_n{\Vert }^2$は射影したデータの分散のデータの個数倍と見れますね。分散の最大化をしていた訳です。

固有値と2次形式の最大最小

さて主成分分析の理論は結局，2次形式の最大値問題に帰着されました。
${\max }_{\mathbf{e}:\Vert \mathbf{e}\Vert =1}(\mathbf{e},X{X}^{\top }\mathbf{e})$
求める$\mathbf{e}$は固有値を使うことで簡単に求めることができます。なお，この$\mathbf{e}$方向の符号付きの長さを主成分，特に1番目の主成分なので第1主成分と呼び，$\mathbf{e}$に直交する平面でもう一度直線の当てはめを行って得られた主成分は第2主成分と呼ばれます。これは二つの直交する直線を当てはめているのでつまり平面の当てはめです。
高次元データでは第3主成分, 第4, 第5, ...と必要な分だけ増やしていきます。

定理 2次形式の最大最小

対称行列$A\in {\mathbb{R}}^{n\times n}$
の固有値を${\lambda }_i$ (だだし${\lambda }_1\le {\lambda }_2\le \cdots {\lambda }_n$とする)
単位ベクトルに対する2次形式$(\mathbf{v},A\mathbf{v})$の最大値は${\lambda }_n$, 最小値は${\lambda }_1$である。

(証明)

対称行列の固有ベクトルは正規直交系をなす様に取ることができます。(対称行列について，異なる固有値に対応する固有ベクトルたちは互いに直交しているので正規化することで正規直交系とすることができる)
この正規直交系をなすベクトルを${\mathbf{u}}_1,{\mathbf{u}}_2,\cdots ,{\mathbf{u}}_n$ただし固有値が${\lambda }_i$と${\mathbf{u}}_i$とが対応するとします。
(対称行列$A$は固有ベクトルをただ$n$個だけ持つ)
また${\mathbf{u}}_i$たちは$n$本の一次独立なベクトルなので${\mathbb{R}}^n$の正規直交基底です。
よって，任意に与えられた単位ベクトル$\mathbf{v}$を${\mathbf{u}}_i$たちの線形結合で表すことができて，

$$\mathbf{v}=\sum _i{c}_i{\mathbf{u}}_i$$

これを代入して

$$\begin{array}{rl}(\mathbf{v},A\mathbf{v})& =\left(\sum _i{c}_i{\mathbf{u}}_i,\sum _i{c}_i{\lambda }_i{\mathbf{u}}_i\right)\\ & =\sum _i{\lambda }_i{c}_i^2\le \sum _i{\lambda }_n{c}_i^2\le {\lambda }_n\end{array}$$

最後は$\mathbf{v}$は単位ベクトルだから${\sum }_i{c}_i^2=1$であることを用いました。
等号成立は$\mathbf{v}={\mathbf{u}}_n$のとき。
よって2次形式の最大値は${\lambda }_n$．同様にして最小値は${\lambda }_1$．

特異値分解

主成分分析では$X{X}^T$の固有ベクトルを求めましたが，特異値分解という手法を用いると少ない計算量で$X{X}^T$の固有ベクトルを求めることができます。

一般の行列$A\ne O$について，

$$\begin{array}{r}\{\begin{array}{rl}& A\mathbf{v}=\sigma \mathbf{u}\\ & A^{\top }\mathbf{u}=\sigma \mathbf{v}\end{array}\\ \sigma >0,\mathbf{u}\ne 0,\mathbf{v}\ne 0\end{array}$$

なる組$(\sigma ,\mathbf{u},\mathbf{v})$は$\mathrm{rank}A$個だけ存在する。
$\mathbf{u},\mathbf{v}$をそれぞれ左特異ベクトル, 右特異ベクトルと呼び，左特異ベクトルは$A{A}^T$の固有ベクトルになっています。また，右特異ベクトルは$A^{\top }A$の固有ベクトルで，$\sigma$は固有値の平方根になっています。（自分で確かめてみよう）
特異値分解をする方が対角化をするよりも一般に高速なので実際には特異値分解で$X{X}^T$の固有ベクトルを求めます。

この特異値，特異ベクトルを用いると行列$A$を次のように分解することができます。

${\sigma }_i,{\mathbf{u}}_i,{\mathbf{v}}_i$を同じ順番に，

$$U=[{\mathbf{u}}_1{\mathbf{u}}_2\cdots {\mathbf{u}}_r]$$

$$V^{\top }=[{\mathbf{v}}_1^{\top }{\mathbf{v}}_2^{\top }\cdots {\mathbf{v}}_r^{\top }]$$

$$\Sigma =\mathrm{diag}({\sigma }_1,{\sigma }_2,\cdots ,{\sigma }_r)$$

$$(\text{ここで}r=\mathrm{rank}A)$$

と並べて

$$A=U\Sigma V^{\top }\text{(特異値分解)}$$

と表せることを指すことが多い。($U,\Sigma ,V$の定義にはいくつか流儀があるようですが，次の節のコラム考慮するなら上の定義のが良いと思う。)
次のようにも表せます。計算する際にはこちらのが便利です。

$$A=\sum _{i=1}^r{\sigma }_i\mathbf{u}{\mathbf{v}}^{\top }$$

またこの記事では${\sigma }_1\ge {\sigma }_2\ge \cdots \ge {\sigma }_r$と並び替えてあるとします。

補足: $r=\mathrm{rank}A$について

まず$\ker (A)=\ker (A^{\top }A)$を示す。
任意に$\mathbf{x}\in \ker (A)$が与えられたとする。$A\mathbf{x}=0$から$A^{\top }A\mathbf{x}=0$となって$\mathbf{x}\in \ker (A^{\top }A)$　よって $\ker (A)\subseteq \ker (A^{\top }A)$

任意に$\mathbf{x}\in \ker (A^{\top }A)$が与えられたとすると，

$$(\mathbf{x},A^{\top }A\mathbf{x})=(A\mathbf{x},A\mathbf{x})=\Vert A\mathbf{x}{\Vert }^2$$

ここで$A^{\top }A\mathbf{x}=0$より$\Vert A\mathbf{x}{\Vert }^2=0$すなわち $A\mathbf{x}=0$
これより，$\mathbf{x}\in \ker (A)$ よって $\ker (A^{\top }A)\subseteq \ker (A)$

ここでは数ベクトルを扱っているので$\dim (\mathrm{im}(A))=\mathrm{rank}A$

次元定理より

$$n=\mathrm{rank}A+\dim (\ker A)またn=\mathrm{rank}(A^{\top }A)+\dim (\ker (A^{\top }A))$$

よって$\mathrm{rank}A=\mathrm{rank}(A^{\top }A)$

ところで右特異ベクトルは$A^{\top }A$の固有ベクトルでした。($A^{\top }A$の固有ベクトルであることが必要。)次に$A^{\top }A$の非負固有値の個数は$\mathrm{rank}(A^{\top }A)$に等しいことを示します。
$A^{\top }A$は半正定値行列であることが示せます(略)
$A^{\top }A$は対称行列なので直交行列で対角化可能で，つまり正規直交系をなす${\mathbf{v}}_1,{\mathbf{v}}_2,\cdots ,{\mathbf{v}}_n$
を列ベクトルに持つ$V$で対角化できて，（ただし，固有値が降順になるように\bm{v}_iたちをとります。）
${\mathbb{R}}^n\ni \mathbf{c}=c_1{\mathbf{v}}_1+\cdots +c_n{\mathbf{v}}_n$に対して

$$\begin{array}{rl}{A}^{\top }A\mathbf{c}& =({\lambda }_1{\mathbf{v}}_1{\mathbf{v}}_1^{\top }+\cdots +{\lambda }_n{\mathbf{v}}_n{\mathbf{v}}_n^{\top })\mathbf{c}\\ & ={\lambda }_1{c}_1{\mathbf{v}}_1+\cdots +{\lambda }_n{c}_n{\mathbf{v}}_n\\ & ={\lambda }_1{c}_1{\mathbf{v}}_1+\cdots +{\lambda }_n{c}_r{\mathbf{v}}_r({\lambda }_i=0(i=r+1,\dots ,n))\end{array}$$

と表せるので$0$でない固有値の数と$\mathrm{rank}(A^{\top }A)$がちょうど等しいことがわかります。対応する固有値が$0$でない固有ベクトルが$A$の右特異ベクトルです。
半正定値行列だったので$0$でない固有値は全て正です。よって${\sigma }_i:=\sqrt{{\lambda }_i}$とすれば良いですね。
最後に，右特異ベクトルと特異値から左特異ベクトルが得られます:

$${\mathbf{u}}_i:=\frac{A{\mathbf{v}}_i}{{\sigma }_i}$$

と定義すると特異値・特異ベクトルの定義を満たしていますね。

コラム　ムーア・ペンローズ型の一般逆行列

逆行列は少なくとも正方行列でないと定義されませんでした。一般の$A\in {\mathbb{R}}^{m\times n}$行列に逆行列は定義されないのでしょうか？
ところで射影行列は射影した空間のベクトルに対しては単位行列$E$のように働きます。

$$\begin{array}{rl}{P}_{\mathcal{U}}{\mathbf{x}}_{//}& ={\mathbf{x}}_{//}\\ E{\mathbf{x}}_{//}& ={\mathbf{x}}_{//}\end{array}$$

正方行列$A\in {\mathbb{R}}^{n\times n}$の逆行列とは

$$AB=E_n$$

なる$B$のことでした。

これを拡張し，$A\in {\mathbb{R}}^{m\times n}$に対し，
$A$の列空間($A$を列ベクトルに分解した時，その列ベクトルの線型結合で張られる部分空間です)への射影行列を
$P_{\mathcal{U}}$とした時

$$AB=P_{\mathcal{U}}$$

なる$B$をムーア・ペンローズ型の一般逆行列と呼びます。
ところで$P_{\mathcal{U}}$は実は左特異ベクトルの張る空間への射影行列と一致します。
簡単に説明すると，特異ベクトルは対称行列の固有ベクトルでもあったので正規直交系に選べます。
${\mathbb{R}}^n\ni \mathbf{c}=c_1{\mathbf{v}}_1+\cdots +c_n{\mathbf{v}}_r$に対して

$$A\mathbf{c}=U\Sigma V^{\top }\mathbf{c}=\sum _{i=1}^r{\sigma }_i{c}_i{\mathbf{u}}_i$$

であり，$\sigma \ne 0,c_i$は任意であるからです。
さて一般逆行列$A^{-}$は次のよう定義されます。

$$A^{-}:=B=V{\Sigma }^{-1}{U}^{\top }\text{ここで}{\Sigma }^{-1}=\mathrm{diag}({\sigma }_1^{-1},{\sigma }_2^{-1},\cdots ,{\sigma }_r^{-1})$$

$U$は直交行列であることに注意すると，実際に，

$$A{A}^{-}=U{U}^{\top }={\mathbf{u}}_1{\mathbf{u}}_1^{\top }+{\mathbf{u}}_2{\mathbf{u}}_2^{\top }+\cdots +{\mathbf{u}}_r{\mathbf{u}}_r^{\top }=P_{\mathcal{U}}$$

途中で

$$U{U}^{\top }=\left[\begin{array}{cccc}{\mathbf{u}}_1& {\mathbf{u}}_2& \cdots & {\mathbf{u}}_n\end{array}\right]\left[\begin{array}{c}{\mathbf{u}}_1^{\top }\\ {\mathbf{u}}_2^{\top }\\ ⋮\\ {\mathbf{u}}_n^{\top }\end{array}\right]=\sum _i{\mathbf{u}}_i{\mathbf{u}}_i^{\top }$$

を使いました。${\mathbf{u}}_i$は列ベクトルなので自明ではありませんが成分計算で簡単に示せます。

この一般逆行列を用いて最小二乗法の解を直接書き下すことができます。また，最小二乗解を図形的に解釈すると射影を用いて直感的に捉えることが出来ます。
詳しくは参考文献[1]や関連するウェブサイトを探してみてください。

主成分分析って実際何をすればいいの

これまで理論よりに解説してきました。主成分分析をする手順を実際にまとめてみます。

1. 縦ベクトルのデータ${\mathbf{x}}_1,{\mathbf{x}}_2,\cdots ,{\mathbf{x}}_N$を並べて行列$X=[{\mathbf{x}}_1{\mathbf{x}}_2\cdots {\mathbf{x}}_N]$を作ります。

2. $X$を特異値分解し，左特異ベクトルを求めます。

3. 第k主成分まで欲しい時は 対応する特異値が大きい順に左特異ベクトルを$k$個とり並べて行列$A^{\top }=[{\mathbf{u}}_1{\mathbf{u}}_2\cdots {\mathbf{u}}_k]$を作ります。

4. ベクトル${\mathbf{x}}_{\mathrm{proj},k}=A{\mathbf{x}}_k$の$i$番目の成分が第$i$主成分になっています。まとめてデータ全体を射影するときは$X_{\mathrm{proj}}=AX$とすれば良い。

実践編

最後に簡単にPythonで主成分分析をしてみましょう。
まずは2次元のデータで行います。

$c_1$と書いてあるのが第1主成分の軸の方向ベクトルです。（$c_2$は第2主成分の軸の方向ベクトル。データを$c_2$方向に射影するよりも$c_1$方向に射影する方が分散が保たれていますね。）

かなり人工的な設定ですが，標準正規分布に従う乱数を各要素にもつ$2\times 200$行列$X$を作り，最初は大体円形に分布しているのを$x$軸方向に引き伸ばし回転させることで斜め楕円型の分布を得ます。上の理論パートの説明同様に$X$の列ベクトルがデータの点です。
次に3次元のデータで主成分分析してみます。3次元から2次元に削減します。

3本のベクトルが主成分の軸です。$c_i$が第$i$主成分の軸になっています。
第2主成分まで抽出する(灰色の平面に射影する)と次図のように情報を損なわずに次元が減らせていそうです。

実際に各主成分の分布が次です。第1主成分から第3主成分に行くに従って分散が小さくなっています。

Pythonのコードは

import numpy as np
import matplotlib.pyplot as plt

angle = np.pi / 5
stretch = 5
m = 200

np.random.seed(3)
X = np.random.randn(2, m) / 10
X = np.array([[stretch, 0],[0, 1]]) @ X # Xを伸縮
X = np.array([[np.cos(angle), -np.sin(angle)], [np.sin(angle), np.cos(angle)]]) @ X # Xを回転
X_g = np.mean(X, axis=1) # 重心
# センタリング。　ブロードキャストの都合で転置同士で計算してさらに転置を取っているが
# 気持ちは X_centered = X - X_g 
X_centered = (X.T - X_g.T).T 

U, s, Vt = np.linalg.svd(X_centered) # 特異値分解
u1, u2 = -U[:, 0], U[:, 1]# 直線の基底なのでどちらを向いていても良い

plt.plot(X_centered[0, :], X_centered[1, :], "bo", alpha=0.7)
plt.plot(0, 0, "o", color="red")
plt.arrow(0, 0, u1[0], u1[1], head_width=0.1, linewidth=3, length_includes_head=True, head_length=0.1, fc='k', ec='k')
plt.arrow(0, 0, u2[0], u2[1], head_width=0.1, linewidth=3, length_includes_head=True, head_length=0.1, fc='k', ec='k')

plt.text(u1[0] + 0.1, u1[1] - 0.05, r"$c1$", fontsize=22)
plt.text(u2[0] + 0.1, u2[1], r"$c_2$", fontsize=22)
plt.xlabel("$x_1$", fontsize=18)
plt.ylabel("$x_2$", fontsize=18, rotation=0)
plt.grid(True)
plt.show()

3次元版

import numpy as np
import matplotlib.pyplot as plt

angle1 = np.pi / 6
angle2 = -np.pi / 12
stretch1 = 25
stretch2 = 10
m = 200

np.random.seed(3)
X = np.random.randn(3, m) / 15
X = np.array([[stretch1, 0, 0],
              [0, stretch2, 0],
              [0, 0, 1]]) @ X  # 伸縮
X = np.array([[1, 0, 0],
              [0, np.cos(angle1), -np.sin(angle1)],
              [0, np.sin(angle1), np.sin(angle1)]]) @ X  # x軸周りに回転

X = np.array([[np.cos(angle2), 0, np.sin(angle2)],
              [0, 1, 0],
              [-np.sin(angle2), 0, np.cos(angle2)]]) @ X  # y軸周りに回転
              
X_g = np.mean(X, axis=1)  # 重心

# センタリング。　ブロードキャストの都合で転置同士で計算してさらに転置を取っているが
# 気持ちは X_centered = X - X_g
X_centered = (X.T - X_g.T).T

U, s, Vt = np.linalg.svd(X_centered)

# 左特異ベクトルを抽出
u1 = -U[:, 0]  # 逆を向いていても良い。図の見やすさのために逆にする。
u2 = U[:, 1]
u3 = U[:, 2]
normal_vector = np.cross(u1, u2) # 平面の点を得るために使う

xmin,xmax=-4.0, 4.0
ymin,ymax=xmin,xmax
zmin,zmax=xmin,xmax

x = np.linspace(xmin - 0.5, xmax + 0.5, 1000)
y = np.linspace(ymin + 1.0, ymax - 1.0, 1000)
xx, yy = np.meshgrid(x, y)
zz = - 1 / normal_vector[2] * (normal_vector[0] * xx + normal_vector[1] * yy)

fig = plt.figure()
ax = fig.add_subplot(projection='3d', aspect='equal')
ax.set_xlim3d(xmin, xmax)
ax.set_ylim3d(ymin, ymax)
ax.set_zlim3d(zmin, zmax)
ax.plot_surface(xx, yy, zz, color="gray", alpha=0.4)

ax.scatter(X[0, :], X[1, :], X[2, :], color="blue", alpha = 0.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
# 主成分の軸の描画
ax.quiver(*[0, 0, 0], *u1, color="black", length=3, linewidth=2.5)
ax.quiver(*[0, 0, 0], *u2, color="black", length=3, linewidth=2.5)
ax.quiver(*[0, 0, 0], *u3, color="black", length=3, linewidth=2.5)
plt.savefig("01.png", format="png", dpi=800)

# 第1, 2主成分をプロット
A = U.T
X_proj = A @ X_centered
plt.figure()
plt.plot(X_proj[0, :], X_proj[1, :], "bo")
plt.savefig("02.png", format="png", dpi=800)

# 各主成分の分布
xmin2, xmax2 = -4, 4
fig = plt.figure(figsize=(5, 6))
ax1 = fig.add_subplot(3, 1, 1)
plt.xlim(xmin2, xmax2)
ax1.plot([-4, 4], [0, 0], "k-", linewidth=1)
ax1.plot(X_proj[0, :], np.zeros(m), "bo", alpha=0.3)

ax2 = fig.add_subplot(3, 1, 2)
plt.xlim(xmin2, xmax2)
ax2.plot([-4, 4], [0, 0], "k-", linewidth=1)
ax2.plot(X_proj[1, :], np.zeros(m), "bo", alpha=0.3)

ax3 = fig.add_subplot(3, 1, 3)
plt.xlim(xmin2, xmax2)
ax3.plot([-4, 4], [0, 0], "k-", linewidth=1)
ax3.plot(X_proj[2, :], np.zeros(m), "bo", alpha=0.3)
plt.savefig("03.png", format="png", dpi=800)

補足:
本記事では[1]に合わせてデータ点を列ベクトルとして横に並べて行列$X$を作りましたが、[3]をはじめとして転置したものを$X$とすることも多いです。また後者の方がPythonで書く際にブロードキャストが使える範囲が広く便利です。

参考文献

線形代数は以下2冊の本，特に前者の本を参考にしました。どちらもかなり面白い本です。

[1]『線形代数セミナー』https://www.kyoritsu-pub.co.jp/book/b10003952.html

[2]『これならわかる深層学習入門』https://www.kspub.co.jp/book/detail/1538283.html

実装例はこの本を参考にしました。

[3]『scikit-learn，Keras，TensorFlowによる実践機械学習 第2版』https://www.oreilly.co.jp/books/9784873119281/


t-SNEの例はこちらからの引用です。
[4] https://qiita.com/segavvy/items/fe530927df30732e2a46

次は統計について書こうかな