ホーム » 「PCA」タグがついた投稿

タグアーカイブ: PCA

PythonでBayesian PCAの実装に失敗した話

こんにちは。

本日はいつものように「実装しましたぁ」ではなく、「実装に失敗しましたぁ」という話になります。 こちらが論文でPRML(Pattern Recognition and Machine Learning)の作者のものです。非常に難しい。本当に難しかったです。ですがなんとかそのアルゴリズムがわかってコードも失敗しながらですが通ったので失敗談共有という形で記事にします。この記事を見て実装に成功した方がいればぜひコメントお待ちしております。

PCA

PCAを知らない方のために説明しますと Principal Component Analysis(主成分分析) と言いまして次元圧縮方法の一つです。恐らく最もポピュラー。過去にその実装の記事をあげましたのでリンクを貼っておきます。 固有ベクトルへの射影、以上です。このPCAの発展版として今回、BPCAをやります。

BPCA

みんな大好き「隠れ変数」のようなものを導入します。何かと言いますと
上の図においてt_nが観測データ、x_nが隠れ変数です。どうやって次元を削減しているかと言いますとProbablistic PCAという理論に基づいており、t_n, x_nの次元をそれぞれn>m>0として観測値t_nは「実はもっと小さい次元から生成されているはずであり、x_nの線形変換+ノイズという過程を経てt_nとして現れている」というものです。そこで論文ではガウス分布からサンプルを生成してヒントン図を用いて検証を行っています。

ヒントン図

アレイ(行列)の成分の大小を色とブロックサイズでの可視化を試みた図です。

実装編 Python

論文をベースに求める結果としては
となって欲しかったのですが綺麗に3列残りませんでした。。どうしてでしょう。。。コメントお待ちしております。
僕と違う実装ですが成功例がQiitaにあったのでそちら含め参考リンクをおいて置きます。でわ

References

ねぇPython、PCAって何?(実装編+)

こんにちは。では前回の内容を踏まえて空間を変えないPCAを行いましょう。今回使うのは

顔のデータ

各次元4096で400サンプル。これを空間を保ち64次元に落とす。つまり4096次元の空間で64本の固有ベクトルがはる空間へ射影する。

これでPCAは一件落着。

 

ニューラルネットワークに深入りしていく予定。

でわ。

ねぇPython、PCAって何?(理論編+2)

こんにちは。前回、射影行列を学びました。n次元ベクトルaを考えよう。行列Vの列ベクトルが作る空間の中でもっともaに近いものはVcという形で表された。ただし、cは縦ベクトルとすると

(1)   \begin{align*} Vc =  V(V^TV)^{-1}V^Ta \end{align*}

となりこれがもっともaに近いヤツで

    \[V(V^TV)^{-1}V^T \]

Vの列ベクトルがはる空間への射影行列

だった。今回はそのVはどんなやつ?か考える。つまりどんなベクトルが作る空間へ射影すればいいんだ?ということ

データ\in \mathbb{R}^m(m次元)

平面を貼るベクトルをU=[v_1, \cdots, v_n]

ただしv_i \in \mathbb{R}^m , i=1,\cdots,n

とする。目標はUを見つけ出すこと。計算の前にいくつか確認しておこう。

  1. 実対称行列の固有ベクトルは直行
  2. 実対称行列は対角化可能
  3. グラムシュミットで正規直交基底が作れる
  4. tr(AB) = tr(BA)
  5. x^T A x = tr(x^T A x)

まず、3番より射影行列P_wは次のようになる

    \[P_w = U(U^TU)^{-1}U^T = UU^T\]

ではどうやってUを決定するかだが、前回同様、分散を最大化させる作戦でいこう。射影後の分散が最大化されるようにUを決定したいので

(2)   \begin{align*} \sum_{i=1}^{N} \|P_w x^{(i)} \|^2 &=& \sum_{i=1}^{N} \|UU^Tx^{(i)} \|^2     \\ &=& \sum_{i=1}^{N} (UU^Tx^{(i)})^T(UU^Tx^{(i)})     \\ &=& \sum_{i=1}^{N} x^{{(i)}^T}UU^TUU^Tx^{(i)}     \\ &=& \sum_{i=1}^{N} x^{{(i)}^T}UU^Tx^{(i)}         \\ &=& \sum_{i=1}^{N} tr(x^{{(i)}^T}UU^Tx^{(i)})       \\ &=& \sum_{i=1}^{N} tr(U^Tx^{(i)}x^{{(i)}^T}U)     \\ &=& tr(U^TXX^TU)      \\ &=& tr(U^T P^TDPU)      \\ &=& tr((PU)^T D (PU)) \end{align*}

ここで、対称行列XX^TP^TDPと変形しました。これを行列の対角化といいます。ここで、Pは固有ベクトルを基底とした空間から標準基底の空間への変換行列なんですが、重要なことは基底(軸)の取り方によって座標が変わるということでしたね。これを踏まえて

    \[tr((PU)^T D (PU))\]

を最大化させるUの形を考えよう!例えば

    \[  tr \left(  \begin{pmatrix} a  & b & c  \\ d  & e & f  \\  \end{pmatrix}   \begin{pmatrix} 12 & 0 & 0 \\ 0 & 7 & 0 \\ 0 & 0 & 1 \\   \end{pmatrix} \begin{pmatrix} a & d  \\ b & e  \\ c & f    \end{pmatrix}  \right)\]

において[a,b,c], [d,e,f]を正規直交とするときどんなものがトレースを最大化させるか?

e_1, e_2

しかないのは明らか。よってUの列ベクトルには固有ベクトルが並べばいいことがわかります。実際、例えばU=p_1が最大固有値 \lambda_1に対応する固有ベクトルとする時Pp_1e_1となる(Pが座標変換するから)

(3)   \begin{align*} tr((PU)^T D (PU)) = tr((Pp_1)^T D (Pp_1)) = tr(e_1^T D e_1) = \lambda_1 \end{align*}

となるからです。従って射影行列P_w=UU^TUは列ベクトルとして固有ベクトルを並べればいいことがわかりました。

次回はこれを実装で確かめましょう。

 

でわ

 

ねぇPython、PCAって何?(理論編+1)

こんにちは。理論編で主成分分析についてはみなさん理解できたと思います。が。本日は少し応用です。

とりあえず復習から始めましょう。

PCAは固有ベクトルへの射影

そうです。たとえば400次元から2次元へのPCAは固有値を大きい方から2つ選び、対応する固有ベクトルに射影です。つまり

軸を新しく定義している

つまり

空間が変わっている

 

まさにこの図、3次元のデータ(左)にPCAをかけて2次元データ(右)にしている。軸が異なるどころかその数も変わっていることより空間が変わったことは明らか。

しかし、今回は

次元を保ったまま空間を押しつぶす

をします。

重要なのは次の考え方のみ!

適当に作った例だと

    \[ [1,2,3] \rightarrow [4,-2]  \]

が前回のPCAに対し

    \[ [1,2,3] \rightarrow [4,-2, 0]\in H\]

のように次元を保ったまま平面上の点に置き換える。上の図だとc_1,c_2 \in \mathbb{R}に対し

    \[ a' =  c_1v_1 + c_2 v_2  = c_1 \begin{pmatrix} v_{11} \\ v_{12} \\ v_{13}  \end{pmatrix} + c_2 \begin{pmatrix} v_{21} \\ v_{22} \\ v_{23}  \end{pmatrix} \]

とできる。ただし、a,a',v_1,v_2\in \mathbb{R}^3である。

 

まとめると、

あるベクトルたちの線型結合で全部書き直そう

ということ。上では3次元ベクトル2本(v_1,v_2)用いて平面を構成した。aを平面H上の点a'としようということ、これを全データ点で行う。計算してみよう。V = [v_1,v_2]とするとVcv_1,v_2が貼る平面状の点を表すことができるので

    \[\|a - Vc\|^2  \]

を最小化するようなcを見つければいい。このときVc=a'となる。ここで

    \[\frac{\partial}{\partial x} x^T A x = (A+A^T)x\]

    \[ \frac{\partial}{\partial x} a^Tx = a^T, \frac{\partial}{\partial x} x^Ta = a\]

これらの公式を用いて

(1)   \begin{align*} \frac{\partial}{\partial c} \|a - Vc\|^2 &= \frac{\partial}{\partial c} (a - Vc)^T(a - Vc) \\ &= \frac{\partial}{\partial c} \left(c^TV^TVc - c^TV^Ta - a^TVc +a^Ta  \right) \\ &= \frac{\partial}{\partial c} \left(c^TV^TVc - 2a^TVc +a^Ta  \right) \end{align*}

微分を0とする。

(2)   \begin{align*} \frac{\partial}{\partial c} \|a - Vc\|^2 &= \frac{\partial}{\partial c} \left(c^TV^TVc - 2a^TVc +a^Ta  \right) \\ &= (V^TV + (V^TV)^T)c - 2V^Ta = 2V^TVc - 2V^Ta = 0 \end{align*}

従って逆行列をかけることで

(3)   \begin{align*} c =  (V^TV)^{-1}V^Ta \end{align*}

つまり!

(4)   \begin{align*} Vc =  V(V^TV)^{-1}V^Ta \end{align*}

a'である。

    \[V(V^TV)^{-1}V^T \]

Vの列ベクトル空間がはる空間への射影行列

という。

では数学の時間です。

データ\in \mathbb{R}^m(m次元)

平面を貼るベクトルをV=[v_1, \cdots, v_n]

ただしv_i \in \mathbb{R}^m , i=1,\cdots,n

とする。目標はVを見つけ出すこと。見つかれば上で導出した射影行列を全データにかければ良い。なんだけど、ちょっと長くなるから次の記事で。

 

でわ。

 

READMORE

  1. http://python.astrotech.io/machine-learning/principal-component-analysis.html#id1
  2. http://www.asahi-net.or.jp/~jb2y-bk/NaturalSci/math/daenmen.htm

ねぇPython、PCAって何?(実装編)

こんにちは。英語が非常に難しい」そんな日々が続きます。

主成分分析とは?

共分散行列の固有ベクトルへの射影による次元削減でした。固有ベクトルと固有値をパパッとほしいですね。

これ便利です。Numpyにその計算機あるんです。これを使って実装へレッツゴー

PCAってすごい簡単ですよね。ただ次元とデータが多くなると共分散行列と固有ベクトルの計算がかなり大変そうですが、、、

では簡単にまとめ

  • 共分散行列が大事
  • 各軸は固有ベクトル

最後に

軸30個 の世界から 軸30個の世界を軸2個の世界と見た 世界へ移動

はどうなる?

これは空間を維持した押しつぶしによる射影である。

と意味わからないことが書いていますがこれについてはまた今度。簡単な例だと、りんごを押しつぶして平たくすると(平面化)これは3次元空間を保ったまま2次元に落ちたと考えられるということです。

 

でわ。

ねぇPython、PCAって何?(理論編)

こんにちは。いまだにたまに暑いと感じるkzです。今回はなにをやりましょうかねえ。

このデータは見ての通り2次元(平面)データです。これを1次元(直線)に置き換えたいです。

各点の配置は?

置き換える時に重要なのは距離感です。離れていた点通しが置き換えた先でご近所さんになると元のデータをうまく表現できてるとは言えませんよね。それが次の図です。

離れていたはずがご近所さんになってしまっています。一方で次の図はどうでしょう、

さっきよりはうまく元のデータを表せていますよね。ではこれを先ほどの青いデータで考えると

この赤矢印上に先ほどと同様にして点を置き換えれば\mathbb{R}^2 \Rightarrow \mathbb{R}^1つまり2次元のデータを1次元で説明できることになります。

では本日の本題に入りましょう。

上の赤矢印(軸)

どうやっていい軸を選びましょうか。じーっと見てみましょう。広がっている方向に矢印が伸びていますね。つまりいい軸とはデータを矢印だけで表現したものと捉えることができそうです。では具体的にいい軸(ベクトル)をどう求めるかを考えましょう。

とりえあず求めるベクトルをvとします。今長さは別に興味ないので| v |=1とします。

正射影ベクトルを考えます。下の緑(a_1)はaのbへの射影です。

  1. https://www.geogebra.org/m/mkV7F8Jf

上の図はaをデータ点bを求めたいベクトルとした時にa_2の最小化もしくわa_1の最大化がしたいことになります。この事実は次のgifより明らかです。

なので今回はデータをx_i、ベクトルをbとするとその正射影ベクトルは

    \[ \frac{x_i \cdot v}{|v|}  v\]

です。これを全データで最大化したいので次のようになります。\frac{1}{N-1}は都合上のものです!

    \[\frac{1}{N-1} \sum_{i=1}^{N}  (x_i^T v)^2 = \frac{1}{N-1}  \sum_{i=1}^{N}  (x_i^T v)^T(x_i^T v)\]

    \[=   v^T \left(\frac{1}{N-1}\sum_{i=1}^{N} ( x_i x_i^T)\right) v\]

となります。\frac{1}{N-1} \sum_{i=1}^{N} ( x_i x_i^T)これって何でしょう?一応ですがx_i,vは共に縦ベクトルです。

そうです共分散行列です。これを\Sigmaと書くことにします。すると

    \[v^T \Sigma v  \]

が今回の目的関数です。今まで通りならこれを最大化するだけだったんですが今回は制約があります。

    \[| v | = 1  \]

です。この状況を条件付き極値問題と言います。これを解くには技が必要です。

    \[\texttt{Maximize} ~~~ f(x,y)  \]

    \[\texttt{Subject to} ~~~ g(x,y)=0  \]

を考える、f,g\mathbb{C}^1級(微分可能かつ導関数が連続)とする。点Pが極値をとるならば

    \[\mathcal{L} (x,y, \lambda) = f(x,y) - \lambda g(x,y)  \]

    \[\nabla f (P) - \lambda \nabla g(P) = 0\]

を満たす。あくまで候補点です。要点だけいうと勾配が平行になる点です。

証明や詳しい解説はここにあります。

  1. https://en.wikipedia.org/wiki/Lagrange_multiplier
  2. http://www.yunabe.jp/docs/lagrange_multiplier.html
  3. https://ocw.mit.edu/courses/mathematics/18-02sc-multivariable-calculus-fall-2010/2.-partial-derivatives/part-c-lagrange-multipliers-and-constrained-differentials/session-40-proof-of-lagrange-multipliers/MIT18_02SC_notes_22.pdf

これをラグランジュの未定乗数法と言います。これを使って計算しましょう。

    \[  \mathcal{L} (v, \lambda)  =  v^T\Sigma v - \lambda (| v|)\]

    \[\frac{\partial \mathcal{L}}{\partial v}= (\Sigma+\Sigma^T)v - 2 \lambda v = 0  \]

    \[ \Sigma v = \lambda v\]

おっと、これわ。。。。

固有ベクトルと固有値そのもの

要点をいうとx\neq 0のベクトルが行列Aの固有値であり\lambdaがその固有ベクトルとは

    \[Ax = \lambda x\]

となることをいい、行列Aによる作用が伸縮になるベクトルを固有ベクトルといいその伸縮率を固有値といいます。

もう少し別の言い方をすると。行列をかけた時に向きが変わらないベクトルとその伸縮率です。

詳しくは下のリンクを見てください。

  1. https://en.wikipedia.org/wiki/Eigenvalues_and_eigenvectors
  2. https://qiita.com/kenmatsu4/items/2a8573e3c878fc2da306

したがって今回私たちが求めていたベクトル、いい軸は

データの共分散行列の固有ベクトル

とわかりました。つまりデータを固有ベクトルという軸へ射影することで得た新たなデータ点は

データの次元を落としつつ最大限に元のデータを表現している新しいデータ

ということになります。これがPCAのやってることです。

「PCAは射影である」

今回は2次元から1次元へのPCAでしたが一般にn次元へのPCAも同じです。固有ベクトルを軸としてとってきて射影により座標がデータのPCA後の座標です。2次元へのPCAなら2本の固有ベクトルにデータを射影して新たな座標を定義します。

最後に用語の説明をしておきます。

  • 第n主成分: 固有値が大きい方から数えてn番目のものに対応する固有ベクトルのことをいう。
  • 寄与率: (第n番目固有値/固有値の総和)
  • 累積寄与率: 寄与率の累積

寄与率?固有値?

はい、説明します。固有値は伸縮率といいました。つまり対応する固有ベクトル方向へのデータの散らばり具合です。よって固有値とはデータの広がり具合(分散)を説明します。なので「元データを説明している具合」を分散という観点から考えて寄与としています。

本日はここまで。でわ

参考;

  1. http://cs229.stanford.edu/notes/cs229-notes10.pdf

ねぇPython、前処理って何?

こんちには。良いバックパックがほしい、そんな日々が続きます。

僕は機械学習を勉強してきていろんなアルゴリズムを使ってきました。その都度いろいろなデータ触れてきました。

そもそもデータをそのままアルゴリズムに突っ込むのは正しいのか?

たとえばPCAの場合データの平均を0にしてからPCAにかける必要があります。

ところでデータは不完全な形でやってきます。今回は次のような車のデータを考えます。

みての通り「文字」「NaN」が混在しています。まずはこれをどうにかする必要があります。対処法は主に次の3パターンかと思います。

LabelEncorder

例えば[dog,cat,dog,mouse,cat]を[1,2,1,3,2]に変換する

OneHotEncorder

例えばdog cat monkey dogの4つのデータを[dog, cat, monkey]で表しそれぞれ

    \[[1,0,0], [0,1,0], [0,0,1],[1,0,0]\]

の4つにします

LabelBinarizer

OneHotEncorderとほぼ同じ。

 

さて、本題のデータの前処理なんですが、

  • 標準化(standardization)
  • 正規化(normalization)

の2種類がメインです。これらの手法をスケーリングといいます。重要なのは違いと必要性ですがその前に定義をみてみましょう。

標準化

    \[ x' = \frac{ x - \mu }{ \sigma }  \]

ただし\mu, \sigmaは平均、分散とする。

こうすると平均0、分散1の分布へと変換されます。

正規化

    \[ x' = \frac{ x - min }{ max - min }  \]

ただしmax,minはデータの最大値、最小値とする。

こうすると新たなデータは[0,1]へとスケーリングされます。標準化と違い範囲が固定されるわけではありません。

で、なんでスケーリングする必要があんの?

もちろんです。

スケーリングなんかせんでもいい

という場合ももちろんあります。例えば、生徒会長を決めるための投票データは「支持する」or「支持しない」のバイナリになります。この場合スケーリングをする必要がないのは明らかでしょう。他には男女別の五教科の点数のデータを分析する際、僕は前処理をする必要性を感じません。なぜならデータの範囲は[0,100]ですし、数値がめちゃめちゃ大きいわけでもないからです。

では、いつスケーリングやねん

これは正直難しい問題だと思います。なにを分析したいのかという状況・目的に大きく左右されると思います。ですが、冒頭でも言ったようにPCAなど特定のアルゴリズムは特定の前処理を必要とします

僕の経験上、正規化よりも標準化の方がポピュラーです。というか僕は正規化はしたことないです。理由としては標準化にはセンタリングが含まれるからです。また、そもそも世の中のほとんどデータは正規分布に従っているんです。IQのデータとか身長とか標高でさえ正規分布に従うようです。加えて数学の世界にはいい定理があります。「大数の法則」と「中心極限定理」です。

  1. 大数の法則:        データ超多かったら、その分布の真の平均と実測値の平均が一致!!
  2. 中心極限定理:     データ超多かったら、正規分布に従うと思ってOK!!

これらも標準化がポピュラーな理由の1つでしょう。(1の証明は簡単、マルコフとチェビシェフの不等式を使う)

では次にスケーリングのメリットについてです。1つは学習速度でしょう。次の図をみてください。

 

各軸でデータの範囲が大きく異なる場合、図左のように学習が遅くなることがあります。一方で右図はスラリと学習できています。これより勾配法を使うアルゴリズム

  • ロジスティック回帰/ソフトマックス回帰
  • ニューラルネットワーク/SVM

とかは恩恵を受けるでしょう。

他にはK-meansは一般にはユークリッド距離を用いてクラスタリングします。この際、各特徴を同じ程度重要視して距離を測りたい場合は標準化は効果を発揮するでしょう。(一方、マハラノビス距離は特徴量ごとに距離に重み付けする)

しかし、やはり簡単な問題ではないので何が目的か何をしたいのかを明確にして使い分ける必要がありそうですね。。ここでは標準化・正規化の2つでしたが無相関化(独立主成分分析)などもあります。

 

最後に簡単な実装を載せておきます。

 

READ MORE