PRML 読書会 #15 「12.2 確率的主成分解析」補足資料

PRML 機械学習 読書会 PCA
「パターン認識と機械学習」(PRML)読書会 #15 の 12.2.1「最尤法による主成分分析(確率的主成分解析)」と 12.2.2 「EMアルゴリズムによる主成分分析」の必要最小限の補足資料です。資料本体は超手抜き仕様の予定。

確率的主成分解析(PPCA)の嬉しいところ

  • 多変量ガウス分布の自由パラメータの数を制限(節約)できる
  • EMアルゴリズムを適用できる
    • D が大きいと計算量で有利になる場合も
    • 欠損値を扱える
  • ベイズで ARD(関連度自動決定)。主成分部分空間の次元を適切に決めてくれる
  • 尤度関数で他モデルと比較できる
  • 分類問題にも
  • 生成モデルにも
  • ★「因子分析と密接な関係を」なんて うかつに言うと宗教論争に巻き込まれる？

いいことだらけ！

PPCA モデル

p(\bf{z})=\mathcal{N}(\bf{z}|\bf{0},\bf{I})
p(\bf{x}|\bf{z})=\mathcal{N}(\bf{x}|\bf{W}\bf{z}+\bf{\mu},\sigma^2\bf{I})
where N:データ点の数、D:データ空間の次元、M:主部分空間の次元

  • 周辺化

p(\bf{x})=\int p(\bf{x}|\bf{z})p(\bf{z})d\bf{z}=\mathcal{N}(\bf{x}|\bf{\mu},\bf{C})
where \bf{C}=\bf{W}\bf{W}^T+\sigma^2\bf{I}

最尤法で解く

  • 平均 : \bar{\bf{x}}=\frac{1}{N}\sum_n \bf{x}_n
  • サンプル共分散 : \bf{S}=\frac{1}{N}\sum_{n=1}^N(\bf{x}_n-\bar{\bf{x}})(\bf{x}_n-\bar{\bf{x}})^T
  • 固有方程式: \bf{S}\bf{u}_i=\lambda_i\bf{u}_i
    • 固有値 λ_i(i=1,...,D) は大きい順に並べられているとする
  • σ^2の最尤値: \sigma_{\rm{ML}}=\frac{1}{D-M}\sum_{i=M+1}^D \lambda_i
  • Wの最尤値: \bf{W}_{\rm{ML}}=\bf{U}_M(\bf{L}_M-\sigma^2\bf{I})^{1/2}\bf{R}
    • L_M: Sの固有値 λ_i のうち、M個を成分に持つ対角行列
    • U_M: Sの固有ベクトル u_i のうち、対応するM個を列ベクトルに持つ D×M 行列、
    • R: 任意の直交行列(回転・反転の自由度)

★実装したよ

Tipping and Bishop(1999b). "Probabilistic principal component analysis"

  • W_ML の導出
  • 全ての停留点が W_ML の式で表せる
  • 尤度関数の最大値は、固有値の上位 M 個に属する固有ベクトルを選ぶとき得られる

概略

  • 対数尤度 L(W,σ)=(PRML 12.44 式) とおき、これを最大化するWとσを求める
  • \frac{\partial L}{\partial \bf{W}}=N(\bf{C}^{-1}\bf{S}\bf{C}^{-1}\bf{W}-\bf{C}^{-1}\bf{W})=0 より \bf{S}\bf{C}^{-1}\bf{W}=\bf{W}
    • この計算は結構がんばる(行列偏微分が得意じゃあなかったら)
  • このとき (1) W=0, (2) C=S, (3) C≠S の３つの可能性がある
  • (3) のとき、W の特異値分解 \bf{W}=\bf{U}\bf{L}\bf{V}^T を使って、 \bf{C}=\bf{W}\bf{W}^T+\sigma^2\bf{I}=\bf{U}(\bf{L}^2+\sigma^2\bf{I})\bf{U}^T を得て、\bf{S}\bf{C}^{-1}\bf{W}=\bf{W} に代入すると W_ML の式を得る
    • V がそのまま R になる
    • ★論文では U を D×M 行列、L を M×M 行列となるように SVD しているが、U を D×D行列、L を D×M行列という(多分こっちが一般的な)SVDをしたほうが計算しやすいと思われる
  • (1)(2) は (3) で導出された式の特殊系として表現できる
  • σ_ML の式の導出は容易。
  • L(W_ML,σ_ML)の最大化は \ln\left(\frac{1}{D-M'}\sum_{i=M'+1}^D \lambda_i\right)-\frac{1}{D-M'}\sum_{i=M'+1}^D \ln\lambda_i の最小化と等価であることが言え、固有値を大きい方から取ればＯＫ。

EM アルゴリズムおさらい

  • X:観測変数、Z:隠れ変数、θ:パラメータ
  • ln p(X) を求めるのは難しいが、ln p(X,Z)は容易
  • E-step: p(Z|X, θ^old) を求める
    • M-step の更新式に現れる十分統計量の計算に落とし込む
  • M-step: 対数尤度の期待値 Q(θ,θ^old)=Σ p(Z|X, θ^old)ln p(X,Z) を最大化するθを求める
    • パラメータごとの更新式を導出する

PPCA with EMA

E-step:

  • \mathbb{E}[\bf{z}_n] = \bf{M}^{-1}\bf{W}^T(\bf{x}_n-\bar{\bf{x}})
  • \mathbb{E}[\bf{z}_n\bf{z}_n^T] = \sigma^2\bf{M}^{-1}+\mathbb{E}[\bf{z}_n]\mathbb{E}[\bf{z}_n]^T

M-step:

  • \bf{W}_{\rm{new}}=\left[\sum_{n=1}^N (\bf{x}_n-\bar{\bf{x}})\mathbb{E}[\bf{z}_n]^T\right]\left[ \sum_{n=1}^N\mathbb{E}[\bf{z}_n\bf{z}_n^T] \right]^{-1}
  • \sigma_{\rm{new}}^2 = \frac{1}{ND}\sum_{n=1}^N \left\{|\bf{x}_n-\bar{\bf{x}}|^2-2\mathbb{E}[\bf{z}_n]^T\bf{W}_{\rm{new}}^T(\bf{x}_n-\bar{\bf{x}})+\rm{Tr}(\mathbb{E}[\bf{z}_n\bf{z}_n^T]\bf{W}_{\rm{new}}^T\bf{W}_{\rm{new}})\right\}

★実装したよ

欠損データ

  • ランダム欠損を扱うことが出来る
    • ★「ランダム欠損」ってあまり発生しなさそう。
  • 欠損値を隠れ変数と見なして EM アルゴリズムを適用する
    • x_ni (n 番目のデータ点の i 番目のフィールド) が欠損している場合、この x_ni を Z と一緒に潜在変数とする
  • \mathbb{E}_{\bf{Z},x_{ni}}[\ln p(\bf{X},\bf{Z}|\bf{\mu},\bf{W},\sigma^2)] を展開して、M-step に必要な十分統計量を見いだす
    • \mathbb{E}[x_{ni}],\; \mathbb{E}[x_{ni}^2],\; \mathbb{E}[x_{ni}\bf{z}_n] を E-step の計算に追加すればいい(はず)
  • ★めんどくさくて実装してない……