Table of contents
Open Table of contents
概要
- PRML12章の確率的主成分分析 (Probabilistic PCA; PPCA) について解説・実装しました。
確率的主成分分析 (PPCA)
確率的主成分分析（PPCA）は、主成分分析（PCA）をガウス分布に基づく生成モデルとして再定義したものです。
通常のPCAは、データの分散を最大化する軸を求める手法であり、決定論的な手法です。
これに対しPPCAは、観測されたデータ が低次元の潜在変数 の線形変換とガウスノイズ によって生成されると仮定し、データ から生成モデル（確率分布）のパラメータを推定します。
生成モデルとして定式化することで、尤度に基づくモデルの比較や欠損値の統計的な補完、さらに複数のPPCAを組み合わせた混合モデルへの拡張など、通常のPCAでは難しい拡張が可能になるという利点があります。
生成モデルの定式化
全体の流れを先に話すと、データ と潜在変数 とパラメータ集合 を元に、 と が陽に書き下せる状況を仮定します。
この後に、潜在変数 を積分消去して、パラメータ集合で条件付けられた確率分布 を求めます。
後は最尤推定法やEMアルゴリズムでパラメータ集合 を推定します。
それでは定式化を始めます。
観測された 次元データ が、 次元の潜在変数 から生成されると仮定します。
ただし、 は平均 、分散 （ 次元の単位行列）の正規分布に従う確率変数です。
具体的には、 が以下の線形関係で記述されるとします。
ここで、重み行列 は の行列、 は 次元ベクトル、 は平均 分散 の 次元正規分布に従うとします。
パラメータ集合は とします。
式を見て気づくとは思いますが、 は PCA における主成分ベクトルを並べた行列と直交行列（回転）の積のイメージを持つと良いです（後で を最尤推定し、主成分ベクトルとの関係について言及します）。
式の分布から を一点サンプリングして固定したときの条件付き確率 は以下のようになります。
次に、 式から を積分消去（周辺化）しましょう。
これで が求まりました。
次に、 個のデータを並べた行列 に対する同時確率分布 の対数尤度関数 を求めます。
解析解による最尤推定
を について偏微分して が最大となる解 を求めます。
ここでの式展開のほとんどは The Matrix Cookbook1 を参考にしました。
1) \boldsymbol{\mu} に関する偏微分
を解くと、
導出の詳細
とおくと、
なので、
かつ は正則なので、 から結果が得られます。
2) \mathbf{W} に関する偏微分
先程求めた を 式に代入して、
が得られます。ただし、 であり、式変形に と を用いました。
行列の微分を計算すると、
となります。
トレースの等式の証明
を示します。
なので成り立ちます。
に関する偏微分の成分計算
途中の式変形で と (参考文献 [2] の(124) 式) を用いました。
と同様に を満たす を求めます。
式を満たすには、 である必要があります。
この条件を満たす解を求めるために、以下の場合分けを行います：
- の場合（自明な解）
- 、すなわち の場合
- それ以外の場合（ かつ ）
2-1) \mathbf{W}=0
は自明な解の一つで、このときの対数尤度は
ここで、
- 固有値分解
- は の固有値 を対角成分に並べた対角行列
- は の対角成分 に対応する固有ベクトル を横に並べた行列
を用いました。
2-2) \mathbf{C}=\mathbf{S}
のとき 式は満たされます。
このとき
の特異値分解を と書くことにします。ここで はユニタリ行列, で は対角行列, はユニタリ行列です。この特異値分解から
が得られます。２つの結果を比較すると、
となります。
2-3) \mathbf{W}\neq0 かつ \mathbf{C}\neq\mathbf{S}
特異値分解した を 式に代入すると、
これは の固有値・固有ベクトル の式そのものなので となり が得られます。
PRML の式に合わせるように と の次元を変えると最終的に
が得られます。ここで は の固有値を 個対角成分に並べた行列で、 は の固有ベクトルを 個横に並べた行列です。
3) \sigma に関する偏微分
と を代入して を解くと、
PPCAの最尤推定まとめ
PRML に書いてあることを整理しました。
- 全体の計算量は
- 潜在変数空間は回転不変性を持つ
- のとき、 の各列 (つまり共分散行列 の固有ベクトル ) は でスケール(数値シミュレーションの重み行列の比較で確認)
EM アルゴリズムによる最尤推定
ここでは反復法で最尤推定解を求める方法を説明します。
確率分布 を導入して を書き直すと、
と書けます。EMアルゴリズムは後述するEステップとMステップを交互に反復することで 式を最大化します。
の初期値は とします。
式の導出
E-step
ステップ目において とすると、 式の KL 項は 0 になり、
計算する必要があるのは と で、これらは事後分布 から求められます。
E-step の詳細な式展開
同時分布を展開すると、
式に代入して期待値を計算すると、
途中で および を用いました。
よって、
PRMLでは で としています（ の最尤解は解析的にすぐ分かるので妥当だと思います）。
M-step
E-step では とすることで 式の KL 項を 0 にしました。
M-step では 式の を について最大化します。この最大化の過程で KL 項は 0 ではなくなりますが、 より対数尤度は大きくなるので問題はありません。 式を と について偏微分すると、
が得られます。
数値シミュレーション
実装したコードは PPCA.ipynb にあります。
データは Digits データセット2 を使いました。
このデータセットは、データ数 、特徴量（次元数） です。
主成分の数 は に設定して 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 の処理3の中で、左特異ベクトルを並べた と右特異ベクトルを並べた に対し絶対値が最大の列が常に正になるような符号変換を行うためだと考えられます。
PRMLでも説明されているように、PPCAは潜在空間の回転不変性があるため特に問題にはなりません。
PCA と PPCA での重み行列の比較
PCA の主成分ベクトルを並べた行列 と PPCA の重み が対応するので、下図の (a) と (b) に各行列をプロットしました。(c) は縦ベクトルのノルムをプロットしたグラフです。
で計算したため、重み行列の各列ベクトルが でスケールされていることが確認できました。
EM アルゴリズムの解と解析解の比較
の条件下で EM アルゴリズムで得られる解と解析解について対数尤度の比較を行いました。
下図は対数尤度の推移を表したグラフで、横軸が EMアルゴリズムの 1 step（E-step と M-step を 1 回実行）で、縦軸が対数尤度です。
EMアルゴリズムにより、 の場合は対数尤度が解析解に収束することが確認できました。
ただし、 の場合は局所解にはまったためか、解析解の対数尤度に収束しない結果となりました。