Probabilistic PCA (PPCA) の実装

概要

• PRML12章の確率的主成分分析 (Probabilistic PCA; PPCA) について解説・実装しました。

確率的主成分分析 (PPCA)

確率的主成分分析（PPCA）は、主成分分析（PCA）をガウス分布に基づく生成モデルとして再定義したものです。

通常のPCAは、データの分散を最大化する軸を求める手法です。 これに対してPPCAは、観測されたデータ が低次元の潜在変数 の線形変換とガウスノイズ によって生成されると仮定し、データ から生成モデル（確率分布）のパラメータを推定します。

生成モデルとして定式化することで、尤度に基づくモデルの比較や欠損値の統計的な補完、さらに複数のPPCAを組み合わせた混合モデルへの拡張など、通常のPCAでは難しい拡張が可能になるという利点があります。

生成モデルの定式化

全体の流れを先に話すと、データ と潜在変数 とパラメータ集合 を元に、 と が陽に書き下せる状況を仮定して、データ分布 を求めます。 この後に、潜在変数 を消去して、パラメータ集合で条件付けられた確率分布 を求めます。 ここまでできたら、後は最尤推定法なり、EMアルゴリズムでパラメータ集合 を推定します。

それでは定式化を始めます。 観測された 次元データ が、 次元の潜在変数 から生成されると仮定します。 ただし、 はシンプルな正規分布に従う確率変数です。

具体的には、 が以下の線形関係で記述されるとします。

ここで、重み行列 は の行列、 は 次元ベクトル、 は平均 分散 の 次元正規分布に従うとします。 パラメータ集合は とします。 式を見て気づくとは思いますが、 は PCA における主成分ベクトルを並べた行列と直交行列（回転）の積のイメージを持つと良いです（後で を最尤推定し、主成分ベクトルとの関係について言及します）。

式の分布から を一点サンプリングして固定したときの条件付き確率 は以下のようになります。

次に、 式から を積分消去（周辺化）しましょう。

次に、 個のデータを並べた行列 に対する同時確率分布 の対数尤度関数 を求めます。

解析解による最尤推定

を について偏微分して が最大となる解 を求めます。 ここでの式展開のほとんどは The Matrix Cookbook¹ を参考にしました。

1) に関する偏微分

とおくと、

なので、

を満たす を求めれば良いので、 から

2) に関する偏微分

先程求めた を 式に代入して、

が得られます。ただし、 であり、式変形に と を用いました。

について の偏微分は

途中の式変形で と (参考文献 [2] の(124) 式) を用いました。 行列の微分で書き直すと、

と同様に を満たす を求めます。

i)

は自明な解の一つで、このときの対数尤度は

ここで、

• 固有値分解
  • は の固有値 を対角成分に並べた対角行列
  • は の対角成分 に対応する固有ベクトル を横に並べた行列

を用いました。

ii)

のとき 式は満たされます。 このとき

の特異値分解を と書くことにします。ここで はユニタリ行列, で は対角行列, はユニタリ行列です。この特異値分解から

が得られます。２つの結果を比較すると、

となります。

iii) かつ

特異値分解した を 式に代入すると、

これは の固有値・固有ベクトル の式そのものなので となり が得られます。 PRML の式に合わせるように と の次元を変えると最終的に

が得られます。ここで は の固有値を 個対角成分に並べた行列で、 は の固有ベクトルを 個横に並べた行列です。

3) に関する偏微分

と を に代入すると、

今までと同様に を満たす を求めると、

その他

PRML に書いてあることを雑にまとめました。

• 全体の計算量は
• 潜在変数空間は回転不変性を持つ
• のとき、 の各列 (つまり共分散行列 の固有ベクトル ) は でスケール(数値シミュレーションの重み行列の比較で確認)
  • から明らか

EM アルゴリズムによる最尤推定

ここでは反復法で最尤推定解を求める方法を説明します。 確率分布 を導入して を書き直すと、

と書けます。EMアルゴリズムは後述するEステップとMステップを交互に反復することで 式を最大化します。 の初期値は とします。

E-step

ステップ目において としたとき、 式の KL 項は0になるので

この E ステップで必要となるのは と の同時分布 であり、これは

この結果を 式に代入すると、

が得られます。

よって、

PRMLでは で としています（ の最尤解は解析的にすぐ分かるので妥当だと思います）。

M-step

E-step では とすることで 式の KL 項を 0 にした。

M-step では 式の を について最大化する。この最大化の過程で KL 項は 0 ではなくなるが、 より対数尤度は大きくなるので問題はない。 式を と について偏微分すると、

が得られます。

数値シミュレーション

実装したコードは PPCA.ipynb にある。 データは Digits データセット² を使った。

このデータセットは、データ数 、特徴量（次元数） です。
主成分の数 は に設定して PPCA の最尤推定を行いました。

また、 の特異値分解（SVD）を考えると、 を と分解できますが、この （ 次元のユニタリ行列）は任意に選ぶことができます。 今回は式をシンプルにするため、 を単位行列 とみなしました。

PCA と PPCA での潜在表現の比較

データを第一・第二主成分ベクトルに射影した結果を下図に示す。 横軸は 、縦軸は で、各色はデータ のラベル に対応している。

• (a) : の平均を使って を に変換して最初の2成分をプロット
• (b) : を PCA の主成分ベクトル を並べた行列 と思って各主成分ベクトルに射影し、最初の2成分をプロット
• (c) : PRML通りに自分で実装した PCA を使って第二主成分までプロット
• (d) : scipy.decomposition.PCA を使って第二主成分までプロット

(a) と (b) はほとんど同じように見えますが、値のレンジが大きく異なります。これは式を比較すると分かるように から に変換するときに (a) は標準化されていて (b) はされていないからです。

(c) と (d) を比較すると、原点対称関係にあることが分かります。これは scipy.decomposition.PCA の処理³の中で、左特異ベクトルを並べた と右特異ベクトルを並べた に対し絶対値が最大の列が常に正になるような符号変換を行うためだと考えられます。

PRMLでも説明されているように、PPCAは 潜在空間の回転不変性があるため特に問題にはなりません。

PCA と PPCA での重み行列の比較

PCA の主成分ベクトルを並べた行列 と PPCA の重み が対応するらしいので、下図の (a) と (b) に各行列をプロットしました。(c) は縦ベクトルのノルムをプロットしたグラフです。

で計算したため、重み行列の各列ベクトルが でスケールされていることが確認できました。

EM アルゴリズムの解と解析解の比較

条件下でEM アルゴリズムで得られる解と解析解について対数尤度の比較を行った。 下図は対数尤度の推移を表したグラフで、横軸が EMアルゴリズムの 1 step (E-step と M-step を 1 回実行) で、縦軸が対数尤度です。

EMアルゴリズムによる対数尤度は解析解に収束することが確認できました。 ただし、 の場合は収束までにかなり時間がかかりそうです。

Footnotes

1. K. B. Petersen and M. S. Pedersen, The Matrix Cookbook, (2012) ↩

2. Pen-Based Recognition of Handwritten Digits Data Set ↩

3. scikit-learn/sklearn/decomposition/_pca.py ↩