主成分分析(Principle Component Analysis, PCA)は、機械学習の実務等でよく使われる手法です。
主成分分析をするタイミングは非常に多くありますが、基本的に主成分分析の目的は、高次元データを低次元に変換することが一番大きな目的になります。いわゆる、次元削減というやつですね。
次元削減をすることで、扱いたいデータの大まかな傾向を人間が可視化し確認することができます。例えば、100次元以上あるようなデータを扱うことになった際、これらの次元を1つ1つ確認することは大変ですが、この多数の次元を、データが持つ関係性や傾向をできる限り保持したまま、2次元や3次元に変換することで、データの傾向を掴むことができるようになります。
また機械学習で実際にモデルを組む際には、次元の呪いという問題があります。
次元の呪いは、機械学習モデルを学習させる際に、データの次元数が多くなればなるほど、モデルを学習するためのデータセット量が大量数必要になってしまう現象です。
この次元の呪いを回避するために、前処理として主成分分析を行い、データに含まれる基本的な特徴を保持したまま次元削減を行うことで、少ないデータセットで精度の高い機械学習モデルを構築したりします。
例えば、PCAによってデータの次元削減を行なった上で、回帰モデルやk-meansといったクラスタリング手法を用いてモデルを構築するようなことは業務でも多くあります。他にも、パターン認識、画像識別、時系列分析などの論文に登場するような手法となっています。
今回は、実際のデータ解析や機械学習モデルの構築の際に使われることが多い主成分分析(PCA)について、そのメリットや数式的な背景が理解できるよう、実際に具体例を交えながら解説をしていきます。
- 主成分分析の数学的な背景について解説
- Irisデータセットを用いて、実際に主成分分析を行う
主成分分析の基本的な考え方
主成分分析の概要
まず、最初に主成分分析の手法の概要について抑えましょう。
主成分分析について下記のポイントを押さえておくことが重要です。
- 主成分分析は多次元データを低次元に変換する手法
- 低次元に変換する際に、データが持っている情報量をできる限り損なわないように変換する
要約すると、今持っているデータが多次元（高次元）を、できる限り情報を損なわないように低次元にする手法なんだな、というイメージを最後まで持っておいてください。
ちなみに、主成分分析で低次元にする場合は色々なユースケースがありますが、高次元のデータをまず目視で確認したいような場合には、主成分分析を用いてまず2~3次元に落とすようなことはよくやります。いわゆるvisualizationというものです。
主成分軸の決め方
主成分分析の目的は、高次元データを低次元に変換する手法であると書きましたが、今持っているデータが$d$次元である時、主成分分析を用いて $q (q < d)$に落とす際には、分散が最大化するような軸をとることで、次元数を効率的に減らすことができる手法になっています。
この主成分分析を利用して、新たに取り直した軸のことを主成分軸といいます。
分散を最大化するような軸をとると書きましたが、分散を最大化することで、次元削減をする際に、できる限り失う情報(残差)が少なくなるからです。この辺りは、この記事後半の数式でも出てきますが、まずはどういうことなのか、図を用いて見ていきましょう。
主成分軸の取り方を図で見てみる
こちらの図は、横軸に体重、縦軸に身長をプロットした散布図になっています。データの次元は、体重と身長で2次元データとなっています。あるクラスの学生の体重と身長のデータだとします。
この時、上の散布図に対して、2次元の主成分分析を行うことで、新たな主成分軸としてオレンジ色の軸とピンク色の軸を取り直すことができます。
この様に軸を取り直すと、オレンジ色の軸は「体の大きさ」の様な尺度になり、ピンク色の軸は「肥満度」の様な尺度になります。このように、主成分分析によって、軸の取り方を変えることができます。このオレンジ色の軸が第一主成分に相当し、ピンク色の軸が第二主成分に相当する主成分軸となります。今回は体重と身長の二次元データで主成分分析を行い、第二主成分まで得ることができました。
当然、2次元データを2次元データに変更したので、この例では次元削減はできていません。ただ、主成分分析がやっている軸の取り方のイメージはできたのではないのでしょうか。
当然、次元削減をしたければ、第二主成分であるピンク色の軸を削除すれば、2次元データを一次元に変換することができます。これで、身長と体重の2次元データを1次元のオレンジ色の軸「体の大きさ」に相当する様なデータに変換することができました。なお、ここで注意して欲しいのは、データ列は二次元データから1次元データに変換できていますが、このオレンジ色のベクトル（主成分ベクトル）は、前のデータ列と同じ次元を持っていることは意識しておきましょう（次元削減で削減しているのは、あくまでデータの次元数であり、主成分軸の次元数ではありません）
さて、ここまでの話では主成分軸が取れた前提ですが、この主成分軸は、このデータを主成分軸に下ろしたとき（正射影した時）の値の分散が最大化される様な軸を取得します。
主成分分析の目的は一番最初に記載しましたが、一番の目的は次元削減をすることです。ただ、持っているデータそのものの情報はできるだけ失いたくないため、分散が最大化するような軸の取り方をします。
- d次元データをr次元データに変換する際の、r次元ベクトル
- r次元ベクトル空間での、変換後のデータの値
この2つを得ることができれば、無事主成分分析ができたことになります。
ここまでで、主成分分析の概念的なものを説明してきました。以降ではここまでの考え方を少し数学的に扱ってみましょう。
主成分分析の導出や証明に必要な数学的な知識
数学的な議論をする前に、数学的に必要な前提知識を整理しておきます。ここから先の数学的な話を理解するには、まず下記の項目を理解できていないと確実に振り切られるでしょう。
ここでの項目で不十分な人がいる場合は、まず数学的な勉強を少しした後で、ここに戻ってくるのが良いでしょう。もちろん、業務で主成分分析を扱うだけであれば、数学的な議論はある程度目を通すだけでよいとは思いますが、ちゃんと理解したい人は自分で手を動かして、コードを書いてみて理解するのが一番良いかと思います。
行列・ベクトルの演算に習熟している人はここをスキップしてください。
微分積分・解析学
制約条件つきの最大化問題を扱う際には、下記のラグランジュの不等方程式を用いる。
線形代数
行列の 微分や積分の演算においては、下記の法則が成り立つ。
\frac{\partial}{\partial \bm{x}} \bm{x^T} \bm{x} = 2 \bm{x}主成分分析の数学的な議論
ここからは、主成分分析の数学的な議論を始めていきましょう。
ここで、何の目的もなく、数学的な話を始めると頭が混乱してしまうので、ここでやろうとしている話を整理します。
現状の整理とやろうとしていること
今、扱いたいデータセット $\bm{\mathcal{D}}$ があるとします。データセット $ \bm{\mathcal{D}} $は、$\bm{d}$次元を持つ、$N$個のデータから構成されるデータセットである。
今このデータを解析しようにも、$\bm{d}$次元は十分に大きいため、$\bm{r}$次元にまで下げたいと考え、主成分分析を行うことにした。ただし、$\bm{r} \leq \bm{d}$) 主成分分析は、観測データ$\mathcal{D}$が存在する、$\bm{d}$次元空間において、分散を最大化する様な、軸を$\bm{r}$つ選び、その軸からなる$\bm{r}$次元空間に、観測データを写像する手法である。
上記の記述がこれからやろうとしていることの、日本語の説明である。
以降は数学的な議論に入っていく。
数学的な議論
まず、扱いたい$\bm{n}$個のデータセット$ \boldsymbol{\mathcal{D}} $が存在している。ここに含まれるデータは、$\bm{d}$個の観点から計測されたd次元データであり、この1つ1つを$\bm{x_i}$ で表すと、このデータセット$ \bm{\mathcal{D}} $は、
\bm{\mathcal{D}} = \{ \bm{x_{1}}, \bm{x_{2}} \dots, \bm{x_{i}} \}
と表すことができる。よりわかりやすい表現として、$ \bm{\mathcal{D}} $をベクトル化した変数 $\bm{X} $を導入する。
\bm{X} = (\bm{x_1}, \bm{x_2}, \dots, \bm{x_n})^\top
を考える。また、今回主成分分析で取得したい第一主成分ベクトル $\bm{u}$を次のように定義する。
\begin{equation}
\bm{u} = \begin{pmatrix}
u_1 \\
u_2 \\
\vdots \\
u_n \\
\end{pmatrix}
\end{equation}今回は、簡単のため一旦、d = 2 の二次元データの例で考えると、$\bm{x}$ の $\bm{u}$方向のベクトルの正射影成分 y は、この様になる。
y_i = u_1x_{i1} + u_2x_{i2} = \bm{u^{T}}\bm{x}
と書くことができる。ここで何度も書くが、主成分分析においては、射影成分 $\bm{y} = {y_1, y_2, \dots , y_i}$ の分散を最大化するようなベクトル、つまり、$\bm{u} = (u_1, u_2)^T $を決めることである。
ここで、$\bm{u}$ 方向の射影成分である、$\bm{y}$の分散を考えてみる。（分散という統計量はちゃんとご存知であろうか？）ここでこの射影成分のデータ $\bm{y}$の分散を $ S_{y}$とおくと、
S_y^2 = \frac{1}{n} \sum_{i} ( y_i - \bar{y})^2 である。まず、$\bar{y}$を求めると、
\begin{split}
\bar{y} & = \frac{1}{n} \sum_{i}^n y_i = \frac{1}{n} \sum_{i}^n \bm{u}^T・\bm{x_i} \\
& = \frac{1}{n} \sum_{i}^n (u_1 x_{i1} + u_2 x_{i2}) \\
& = u_1 \frac{1}{n} \sum_{i}^n x_{i1} + u_2 \frac{1}{n} \sum_{i}^n x_{i1} \\
& =u_1 \bar{x_{1}} + u_2 \bar{x_{2}}
\end{split}
となる。よって、求めたい分散 $S_y^2$は、
\begin{split}
S_y^2 & = \frac{1}{n} \sum_{i} ( y_i - \bar{y})^2 \\
& = \frac{1}{n} \sum_{i}^n ( u_1・x_{i1} + u_2・x_{i2} - u_1 \bar{x_{1}} + u_2 \bar{x_{2}} )^2 \\
& = \frac{1}{n} \sum_{i}^n (u_1 (x_{i1} - \bar{x_{1}}) + u_2 (x_{i2} - \bar{x_{2}}) )^2 \\
&= u_1^2 \frac{1}{n} \sum_{i}^n (x_{i1} - \bar{x_{1}})^2 + 2 u_1 u_2 \frac{1}{n} \sum_{i}^n (x_{i1} - \bar{x_{1}}) (x_{i2} - \bar{x_{2}}) ) + u_2^2 \frac{1}{n} \sum_{i}^n (x_{i2} - \bar{x_{2}})^2 \\
& = u_1^2 s_{11} + 2u_1 u_2 s_{12} + u_2^2 s_{22} \\
& = \bm{u}^T S \bm{u}
\end{split}
と表現される。
最後の式変形が少し難しい。逆に式展開して、
\bm{u}^T S \bm{u} = u_1^2 s_{11}+ 2u_1 u_2 s_{12} + u_2 ^2 s_{22}が成り立つことを確認してもらうのが良いと思う。
ここで、sやらSを勝手に導入したが、Sは観測によって得られた2次元データ $\bm{X} $ の分散・共分散行列である。S, sの定義を下記に記載する。
\begin{equation}
S =
\begin{pmatrix}
s_{11} & s_{12} \\
s_{21} & x_{22} \\
\end{pmatrix}
\end{equation}\begin{equation}
s_{jk} = \frac{1}{n} \sum_i^n (x_{ij} - \bar{x_j}) (x_{ik} - \bar{x_k})
\end{equation}と表現される。分散・共分散行列に慣れてないと、上式についてこれないかもしれないので、一度式展開して上の式変形が成り立つことを確認してほしい。
随分と長くなったが、ここで、(8)式より、この$\bm{u}$ 方向の射影成分 $y$の分散 を最大化するには、(8)式の帰結である、$ \bm{u}^T S \bm{u} $ を最大化するような、$\bm{u}$ を見つければ良い。ここで、観測データ$\bm{X}$の分散共分散行列 S は、観測データから導出できる統計量であり、既知であることに注意です。
ここで、$ \bm{u}^T S \bm{u} $を最大化するには、できるだけ大きい $\bm{u}$ を選べば良いという話になるが、今回$\bm{u}$の大きさには興味がなく、その方向に興味があるので、$||u|| = 1 $という制約条件を与える。
すると、上式の最大化問題は、下記の様に定式化される。
\begin{equation}
\text{maximize} \quad \bm{u}^T S \bm{u} \quad \text{subject to} \quad \bm{u}^T\bm{u} = 1
\end{equation}これは、一種の制約がついた最大化問題になり、この様な最大値をとる様な$\bm{u}$ は、ラグランジュの不等方程式を解くことで得ることができます。
ここで、ラグランジュの未定乗数を$\lambda$とすると、下記の様なラグランジュ関数を解くことで、この様な$\bm{u}$ を見つけることができる。
(ここにラグランジュ関数を書く)
ここから主成分分析を多次元に拡張する
今回、ある個体に関する $p$ 次元の変数を持つ観測データ $x_i$ において、
$ \bm{x_i} = ( x_{i1}, x_{i2} \dots, x_{ip})^\top $ があったとする。
このような変数を有するデータが$n$個観測された際には、このデータ群$X$は次のようなデータ行列で表現することができる。
\bm{X} = (\bm{x_1}, \bm{x_2}, \dots, \bm{x_n})^\top
=
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1d} \\
x_{21} & x_{22} & \cdots & x_{2d} \\
\vdots & \vdots & \vdots & \ddots \\
a_{n1} & x_{n2} & \cdots & x_{nd} \\
\end{pmatrix}
ここで、このXデータベクトルは、計画行列（design matrix)と呼ばれており、一般化線形モデルなどで登場する表現形式なので、この場で慣れておくとよと思います。
ここで、n個の観測データにおける平均ベクトル
\begin{equation}
\bar{\bm{x}} = (\bar{x_1}, \bar{x_2}, \dots, \bar{x_n})^\top
\end{equation}
において、すべてのデータから平均値を引いた、下記のようなデータ列$\bar{X}$を考える。
このような$n$個の$p$次元データにおいて、標本分散・共分散行列を考えます。
\begin{split}
\bar{\bm{X}} = (\bm{x_1}-\bar{\bm{x}}, \dots, \bm{x_n}-\bar{\bm{x}})^\top
=
\begin{pmatrix}
(\bm{x_1}-\bar{\bm{x}})^\top \\
(\bm{x_2}-\bar{\bm{x}})^\top \\
\vdots \\
(\bm{x_n}-\bar{\bm{x}})^\top \\
\end{pmatrix} \\
=
\begin{pmatrix}
x_{11} & x_{12} & \cdots & x_{1d} \\
x_{21} & x_{22} & \cdots & x_{2d} \\
\vdots & \vdots & \vdots & \ddots \\
a_{n1} & x_{n2} & \cdots & x_{nn} \\
\end{pmatrix}
\end{split}とすると、この平均化行列を用いることで、分散・共分散行列が下記のように書くことができる。
\begin{equation}
\sum = Var({\bar{X}}) = S = \frac{1}{n}\bar{X}^{\mathrm{T}}\bar{X}
\end{equation}さらにここで、$d$次元空間における単位ベクトル$\bm{a}=(a_1, a_2, \dots, a_d)^\top $
とすると、式（1.3）で示した平均データベクトル$\bar{X}$の$\bm{a}$方向の射影成分$\bm{s}$は、
s_i = (s_{1}, .. ,s_{d})^{\mathrm{T}} = \bar{Xa}
Var({{s}}) = \frac{1}{N}{s}^Ts = \frac{1}{N}(\bar{Xa})^{\mathrm{T}}(\bar{Xa}) = \frac{1}{N}{a}^T\bar{X}^{\mathrm{T}}\bar{Xa} = a^T \sum aここで、この分散が最大となるような単位ベクトル$\bm{a}$は、
\bm{a} \bm{a}^\top = 1
という制約があることを利用して、ラグランジェの未定乗数法より、
\begin{equation}
E(a) = a^TVar({{\bar{X}}})a - \lambda(a^T a - 1)
\end{equation}を最大にする$\bm{a}$を見つける最大化問題に帰着する。ここで、$\lambda$ はラグランジェ未定定数である。ここでaで微分することで、
主成分分析の出力のイメージを掴む
身長と体重のデータセットを利用しましょう。今回は、統計解析ソフトRなどに含まれるDavisのデータセットを利用しました。Davisのデータセットを散布図にプロットするとこのようになっています。（明らかに外れ値となっている異常データは外しています）
実際に、主成分分析を2回行なった結果がこちらです。
右側がある意味、「体の大きさ」、縦軸が「肥満度」を表すような尺度に変換することができた。この際の、基底$a$ ベクトル $ \bm{a} $ は、
\begin{equation}
\bm {a} =
\begin{pmatrix}
0.82863725 & 0.55978594 \\
-0.55978594 & 0.82863725 \\
\end{pmatrix}
\end{equation}しかし、これだと単なる回転に過ぎない。2次元のデータを、2次元の主成分分析よると回転することの証明を行おう。
元のデータ群を下記のように考える。
ちなみに二次元ベクトルにおける回転行列は下記のように定義される。
\begin{equation}
\begin{pmatrix}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta \\
\end{pmatrix}
\end{equation}となる。そのため、この回転行列によって、(x, y )がどこに移動するかというと、
\begin{equation}
\begin{pmatrix}
\cos \theta & - \sin \theta \\
\sin \theta & \cos \theta \\
\end{pmatrix}
\begin{pmatrix}
x \\
y \\
\end{pmatrix}
=
\begin{pmatrix}
x \cos \theta - y \sin \theta \\
x \sin \theta + y \cos \theta \\
\end{pmatrix}
\end{equation}で表せることになる。
さらに、この式と、(1)から、$ \theta = – 34.040995934522 $ となる。
つまり、大まかではあるが、2次元のデータを2次元のデータに主成分分析することは、回転行列によって軸を回転させることに他ならないことが理解できただろう。
Irisデータセットで主成分分析を行う
ここまでで、主成分分析の理論的な説明を長々とやってきました。
次は演習として、Irisデータセットを用いて、主成分分析を実際にやってみましょう。Irisデータセットは3種類のアヤメの花のデータセットで、花弁(petal)、がく(sepal)の長さ(length)・幅(width)の、系4種類のデータ系列からなるデータセットです。
三種類のアヤメは、setosa、versicolor、virginicaで、このような違うがある花です。
Irisデータセットの詳しい解説はこちらをご覧ください。
冒頭で、主成分分析の大きな目的の1つに次元削減と可視化の話をしました。
今回扱うIrisデータセットは4次元のデータなので、これらを主成分分析を用いて2次元に変換したのちに可視化をしてみます。
Irisデータセットの準備
先ほど提示したリンクにIrisデータセットの準備方法は記載しているので、こちらでは簡略化して説明します。
Irisデータセットを用意する方法は多数ありますが、今回はscilkit-learnを用いてデータセットを準備します。
まず必要なライブラリをインポートしたのちに、Irisデータセットを読み込みます。
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from sklearn.datasets import load_iris
from sklearn.decomposition import PCA
iris = load_iris()
df = pd.DataFrame(iris.data, columns=iris.feature_names)
df['label'] = iris.targetsklearn.datasetsモジュールに含まれるload_iris関数を用いることで、手早くIrisデータセットを準備することができます。これらのデータを一旦確認してみます。
len(df)
# => 150
df.head()
データ数が150個あり、それぞれ4系列と正解ラベル(label)から成るデータであることがわかりました。
scikit-learnでIrisデータセットの主成分分析
ここまででデータの準備が整ったので、実際に主成分分析をやっていきましょう。今回は4次元データを2次元に次元削減した上で、可視化してみます。
今回はPythonの科学計算ライブラリであるscikit-learnを用います。
scikit-learnのsklearn.decompositionモジュールに含まれるPCAクラスを利用することで、簡単に主成分分析を行うことができます。
PCAクラスのインスタンスを作成する際のイニシャライザのn_componentsの値に、主成分軸の次元数(今回の場合2)を指定し、できたpcaオブジェクトのfit_transformメソッドにデータを渡してあげることで実行することができます。
pca = PCA(n_components=2)
data = df[["sepal length (cm)", "sepal width (cm)", "petal length (cm)", "petal width (cm)"]].values
trans_pca = pca.fit_transform(data)結果を可視化するのには、matplotlibを利用します。
fig, ax = plt.subplots()
# colormapの設定
cmap = plt.get_cmap("tab10")
color = [cmap(label) for label in df.label]
ax.set_xlabel("PCA axis 1")
ax.set_ylabel("PCA axis 1")
ax.scatter(trans_pca[:, 0], trans_pca[:, 1], color=color, s=10, marker="o")結果はこのようになりました。
4次元のデータを2次元に落とし込めていることがわかります。さらに、3種類のアヤメ間の持つ特徴がうまく集約できていることがわかります。
実際に機械学習を用いる場合では正解データ（上図でいう所のプロット点1つ1つの色)に相当する部分が分からないため、この主成分分析を実行した結果に対し、k-meansやGMM(混合ガウスモデル)などのクラスタリングを行うことで、データの分類を行ったりします。
この辺りは本記事の範疇外となりますが、別の記事で解説しているのでぜひそちらをご覧ください。
主成分分析の課題
主成分分析は、入力空間上のデータが線形構造を有しいている場合は有効である場合がおおいが、データが非線形な構造を背後に有している場合は主成分を適切に見つけるのが難しいとされています。
このような場合にはカーネル関数を用いたカーネル主成分分析(Kernel PCA)を用いることがあります。カーネル関数を用いた主成分分析についてはまた別の記事で解説します。
おわりに
主成分分析を実際にやってみてお分かりかと思いますが、主成分分析によって高次元のデータを2次元に落とし込むことができ、人間によるデータの可視化を簡単にすることができることから、実際のデータ解析や機械学習の実務の現場で非常によく使われる手法です。
よく使われるがゆえに、主成分分析に関係するPythonやRなどのライブラリの機能が充実しており、今回やったようにコードにして数行で実行することができますが、主成分分析の数学的な背景や導出 に登場する数学の知識や行列の知識は非常に重要なので、ぜひ身につけてみてください。
参考文献
- Jolliffe I.T 1986, Principal Component Analysis. New York: Springer-Verlag