主成分分析を実装してみよう

こんにちは、小澤です。

当エントリは「Machine Learning Advent Calendar 2017」の3日目のエントリです。

今回は、次元削減の手法である主成分分析(以下PCA)を実装してみます。

次元削減とは

次元削減は、データの次元すなわち列数を減らすための手法です。

何のためにそんなことをするのかというと、世の中には「次元の呪い」という不吉なものがあります。 データ"サイエンス"というと科学的なことをしている雰囲気があるのに、"呪い"というオカルトな単語が出てくるのはなかなか面白いですね。

次元の呪いとは、高次元空間においてはデータ点のほとんどが外周付近に分布するような状態になってしまう現象のことです。 よくわからないですね。 よくわからないので、以下が参考になります(?)

• 次元の呪い、あるいは「サクサクメロンパン問題」

データ分析では、テーブル形式のデータの各列がそれぞれ1つの次元として扱われるので、特徴の数を増やしていくと数万次元空間のデータを扱うような状態になったりもします。 何でも入るドラえもんのポケットですら4次元です。 数万次元なんて世界、どうなっているのか私にもよくわかりません。

ですが、高次元空間におけるメロンパンはほとんどが皮とのことです。 これはちょっと嬉しいですね。

この他にも、マルチコ(多重共線性)という扱うデータの複数の列間で強い相関があると困るような問題があったり、 そもそも分析対象として、あっても意味のない列なんかは計算時間を増やすだけですので、省いておきたいですね。

こういったときに役立つのが次元削減と呼ばれる手法になります。

PCAとは

PCA概要

さて、ではPCAとはどういうものなのでしょうか？ 例えば以下のようなデータがあったとします。 明らかに2つの軸に相関がありますね。

これの軸を次の図のように赤ものに変換します。

これは、どのような変換を行っているかというと、1つ目の軸がもっとも分散が大きくなるように、2つ目の列がその次に分散が大きくなるように... としていきます。 n次元のデータに対して、軸を変えた新たなn次元のデータが出来上がります。

このままでは、次元が削減されてません。 ここがPCAの特徴なのですが、このようにして変換した軸は1つの軸で元のデータよりも多くの情報を表しています。 そのため、上記の図における第1主成分の軸のみで元のデータの何%を表現できるかといった寄与率という考え方ができます。

例えば、元データが10次元だったとして、

• 第1主成分の寄与率が80%
• 第1主成分, 第2主成分までで90%
• 第3主成分までで...

のようなことができます。 元が10次元のデータなので、第10主成分まで使うと100%になりますが、主成分分析がうまくいきやすいデータだと、より少ない次元数で多くを表現できることが分かります。

PCAの軸変換

PCAでは、n次元のデータに以下のような変換を行います

[latex] z_1 = a_{11} x_1 + a_{12} x_2 + ... + a_{1n} x_n \\ z_2 = a_{21} x_1 + a_{22} x_2 + ... + a_{2n} x_n \\ ... ... z_n = a_{n1} x_1 + a_{n2} x_2 + ... + a_{nn} x_n [/latex]

実装してみる

さて、実装してみましょう。

利用するデータ

今回は、線形回帰でも利用したBoston house-pricesを利用します。

実装

まずは、データをロードして共分散行列を求めます。 変数xとyの共分散は以下のような式で求められます。

[latex] s_{xy} = \frac{1}{N}\sum{N}{i} (x_i - \mu_x) (y_i - \mu_y) [/latex]

すべての列同士の共分散を求めたものが共分散行列となります。

boston = load_boston().data

# 各行の平均を求める
boston_mu = [sum(boston[:, row]) / len(boston[:, row]) for row in range(boston.shape[1])]

# 共分散行列
## 初期化
cov = np.zeros([boston.shape[1], boston.shape[1]])

## 共分散を求める
for i in range(boston.shape[1]):
    for j in range(boston.shape[1]):
        cov[i, j] = sum((boston[:, i] - boston_mu[i]) * (boston[:, j] - boston_mu[j])) / boston.shape[0]

これで共分散行列ができました。

では、次にこれを固有値分解してみましょう。 固有値分解は行列Aに対して、

では、これの固有値・固有ベクトルを求めてみましょう(※ 2017/12/04 コメントでご指摘いただいて修正)

[latex] \lambda A = \lambda {\bf x} [/latex]

となるようなλとベクトルxを求めるものになります。 λ, xはn　x　nの正方行列であればそれぞれn個出てきます。

では、実装...えーっ、とこのアドベントカレンダーはNumpyだけ使っていい縛りですが、Numpyでは関数が実装されてるので使っちゃいましょう(え？ずるい？)

# 固有値・固有ベクトルを求める
l, v = np.linalg.eig(cov)

# 固有値をソートした時のインデックスリスト
sl = np.argsort(l)
# インデックスリストの順で固有ベクトルをソート
sv = v[sl]

最後に変換をかけて、第2主成分までをプロットしてみます。

pca = np.dot(boston, sv)
plt.scatter(pca[:, -1], pca[:, -2])

おわりに

今回は、主成分分析の実装をしてみました。

え？固有値・固有ベクトルを求める部分を実装してないだろ！って？ まぁ、numpyに入ってますし...そのうち気が向いたらブログも書きますね。たぶん...

明日は、yoshimによる『TF-IDF理論編』の予定です。 お楽しみに！