【主成分分析（PCA）まとめ】分散最大化・確率的主成分分析・ベイズ主成分分析（MAP推定）まで

はじめに

PCAは色々と勉強になるので、今回ここでザッとまとめておこうと思います。実はこれまでにも何度も取り上げていて、

www.hellocybernetics.tech

www.hellocybernetics.tech

等などあるのですが、その場でその場で場当たり的に比較しているだけなので、今回しっかり色々な立場のPCAを記述しておこうと思います。

データの前提

まず観測データ $\{ {\bf x}_1, \cdots, {\bf x}_N \}$に関して、その平均を $${\bf \bar x} = \frac{1}{N}\sum_{n=1}^N {\bf x}_n $$ として計算しておきましょう。その後、各データ ${\bf x}_n$に関して下記のように変換を施して新しいデータ ${\bf x'}_n$を取り扱うことにします。これは単に平均が原点に来るように平行移動しただけにすぎません。 $$ {\bf x'}_n = {\bf x}_n - {\bf \bar x} $$ このような処理はいつでも簡単に実施できるので、以後、平均が原点に移動されているデータ${\bf x'}$ の方を普通に ${\bf x}$ と記述して扱っていきます。ここで今更ながら各文字の使われ方について、最低限かいておきます。

ボールド体の小文字 $\bf x$ はベクトルです。

ボールド体の大文字 $\bf A$ は行列です。

小文字 $a$ や 大文字 $N$ はスカラーを表します。

${\bf x}_i$ のように下付き文字 $i$ は何かしらのindexを表しています。何のindexかは文脈で。

分散最大化（KL展開）

取り出したい成分について

観測データに対して$\{ {\bf x}_1, \cdots, {\bf x}_N \}$とある単位ベクトル $\bf u$ を使って下記のように直行射影することを考えます。

$$ z_i = {\bf u^T} {\bf x}_i $$

データ点 ${\bf x}_i $ 1つ1つは単なるスカラー $z_i$ になり、 $\{ z_1, \cdots, z_N \}$というデータが得られます。この時 $ \bf u$ をデータ集合$\{z_n\}$ の分散が最大化されるように決定しましょう（もしも分散が小さくなったら、それぞれのデータの違いを数値から見いだせないことになる。例えば全部のデータが仮にほとんど $0$ に近い値となったらツマラナイだろう）。

問題の定式化

すると

$$ V_z = \frac{1}{N} \sum_{n=1}^N z_n^2 $$

を最大化したいという問題になります（元のデータ $\{{\bf x}_n \}$の平均は $0$ に移動されているので、線形変換である直行射影後のデータ$\{z_n\}$の平均も $0$ です。なので分散の計算が楽！）。ここで $\bf u^Tx$ という変換が入っていたことをちゃんと思い出せば、下記のように書き換えることができます。

$$ \begin{align} V_z & = \frac{1}{N} \sum_{n=1}^N ({\bf u^T} {\bf x}_n)^2 \\ & = {\bf u^T} \left( \frac{1}{N} \sum_{n=1}^N {\bf x}_n {\bf x^T}_n \right) {\bf u} \\ & = {\bf u^T} {\bf V_x} {\bf u} \end{align} $$

ここで $\bf V_x$ は $\{{\bf x}_i\}$ の分散共分散行列と呼ばれるものです。これはデータから計算できるものなので、最適化問題にとっては単なる与えられた定数であることに注意しましょう。

問題を解く

後はこの式を $\bf u$ に関して最大化するだけです。私達が求めようとしているのは 単位ベクトル $\bf u$ であったことを思い出し（もしも単位ベクトルに制限しなかったら、上記の式は $\bf u$ をやたらと大きなベクトルにするだけで際限なく増加してしまう）、ラグランジュ未定定数 $\lambda$ を導入して下記のように最適化の目的関数を作ります。

$$ L(\lambda, {\bf u}) = {\bf u^T} {\bf V_x} {\bf u} + \lambda (1 - \bf u^T u) $$

$\bf u$ に関して微分をして $0$ と置くことで

$$ \begin{align} &&{\bf V_x u} - \lambda {\bf u} &= 0 \\ \Leftrightarrow && {\bf V_x u} &= \lambda {\bf u} \end{align} $$

という固有値問題を獲得します。一方で $L(\lambda, {\bf u})$ の $\lambda$による微分を $0$ と置いたものは、元々仮定していた単位行列であれという制約 $\bf u^Tu= 1$ になっています。この制約を上記の固有値問題の式に適応してやるべく、$\bf u^T$ を左から掛けてやれば、

$$ {\bf u^T V_x u} = \lambda $$

という式になり、これは元々最大化したいと考えていた $V_z = \bf u^T V_x u$ に他ならないので、私達が求めている最適化問題は、固有値問題を解いたときの最大固有値が最大値となり、その最大値を満たす $\bf u$ は最大固有値に属する固有ベクトル（第一主成分という）であることが判明しました。また、分散共分散行列（対称行列）の固有ベクトルは互いに直交しているので、2番目に大きな固有値に属する固有ベクトル…、3番目に大きな固有値に属する固有ベクトル…と選んでいけば、複数個の成分を獲得することが出来ます。

元々の多次元信号

$\{ {\bf x}_1, \cdots, {\bf x}_N \}$ の インデックスを観測時間だとみなして5次元信号を見てみましょう。これは正弦波に基づくものと、単なるガウスノイズである信号がそれぞれ入っています。

次元削減後の信号

3次元信号に変換すると、正弦波に基づく信号は一箇所に集約されています。正弦波を重ねれば当然振幅は大きくなるので、この成分は分散が大きくなります。他の成分はどうも雑多な信号です。

寄与率

ここで主成分分析を使う主な応用事例として「次元削減」というものに着目します。 次元削減では、データの本質的な部分を損なわないようにデータを軽量化することが求められます。ひとまずPCAで言えば $\bf x \in \mathcal R^{D}$ であるならば、 $M (< D)$ 個の固有ベクトルを選出して新しい変数を作り出せば、ひとまず次元削減を達成することができるということになります。

ここで何個の固有ベクトルを取り出せば良いのだろうか？ということが疑問になります。 仮に固有値問題解いた結果、$D$個の固有値と固有ベクトルのセット $\{ (\lambda_1, {\bf u}_1), \cdots, (\lambda_D, {\bf u}_D) \}$ と手に入ったとしたら、これを固有値が大きい方から順に何個選べば良いのかという問題です。

PCAでは多くの場合「寄与率」という概念が用いられます。寄与率として固有値 の値 $\lambda_i $を用います。この固有値は、対応する固有ベクトルを用いて射影を行った場合の射影先での分散を表す値でした。分散が異様に小さい（極めつけは常に $0$ である）場合はほとんど意味のない成分であるという考えに基づけば、分散の小さな成分は削ぎ落としても構わないと考えられます。一方で分散が大きな成分は「そのデータを表現するための寄与が大きい」と考えられます。

そのような考えのもとで、$M (< D)$ 個の固有値と固有ベクトルのセットを取り出すことで、

$$ \frac {\sum_{i=1}^M \lambda_i}{\sum_{i=1}^D \lambda_i} $$

を計算したものが（取り出した部分空間の）寄与率と言われるものです。もしもある1つの主成分だけを見積もりたければ、分子に見積もりたい固有値だけを持ってこればいいです。

これは要するに取り出した射影先の分散の和が、全く次元を落とさなかった場合の分散の和に対してどの程度保持されているのかを表しているということになります。もしも $100$ 次元のデータから $3$ 次元のデータ（3個の固有ベクトル）を取り出すだけで、寄与率が $0.9$ などになったらどうでしょうか。ほとんどのデータの情報を残したまま、劇的にデータを軽量化出来ます。ただし残りの $0.1$ が意味のないデータであることなど誰も保証してくれないので、ここはちゃんと自分で考えなければなりません。

各成分が無相関になるような射影先を選ぶ

観測データ を${\bf x}_i \in \mathcal R^D $とします。このデータを並べた行列を ${\bf X} = ({\bf x}_1, \cdots, {\bf x}_N)^{\bf T}$ としましょう（すなわち $\bf X \in \mathcal R^{N\times D}$となっています）。観測データの分散共分散行列は

$$\bf V_x = X^T X \in \mathcal R^{D\times D}$$

で計算ができます。一方で、仮に $\bf U\in \mathcal R^{D \times D}$ という行列でデータを変換した場合には $\bf Z = XU$ という出たを並べた行列を考えることが出来ます。この新しいデータ行列に関して、分散共分散行列は

$$ \begin{align} \bf V_z &= \bf Z^T Z \\ \bf &=\bf U^TX^TXU \\ \bf & =\bf U^TV_xU \in \mathcal R^{D\times D}\\ \end{align} $$

と表すことが出来ます。仮に $Z$ の各成分が互いに無相関であるとするならば$V_z$は対角行列になっていなければいけません。すなわち、適当な対角行列 $\Lambda$ を使って

$$ \bf U^TV_xU = \Lambda $$

となるように $\bf U$ を選びましょうということです。これは「対角化」あるいは「固有値分解」と呼ばれる問題そのものです。それを下記に見ていきましょう。分散共分散行列は対称行列なので$\bf U$を直交行列として扱ってよく（対称行列は必ず直交行列で対角化可能）、

$$ \bf V_x = U\Lambda U^T $$

と固有値分解できます。ここで $\bf \Lambda =diag(\lambda_1,\cdots ,\lambda_D) $ で、$\lambda_i$は $\bf V_x $ の固有値です（縮退も含めれば次元の数 $ D$だけ得られます）。

5次元の観測データ

下記は適当にサンプリングしたデータです。対角に並んでいるのが分散の情報、非対角に並んでいるのが共分散の情報であり、共分散の散布図が斜めに伸びたりしている場合は相関があるという事になります。下記の散布図では相関がある成分、無い成分様々見られますね。

無相関化したデータ

無相関化したデータは共分散を表す散布図が円に近い状態になっています。

そうした後で、$\bf U$から適当な数 $M (<D)$ 個の縦ベクトルを取り出すことで無事次元削減も実施できます。このとき、どのように縦ベクトルを選出するかは、すでに述べた「寄与率」の考え方に従う場合が多いです。すなわち、この中から分散が大きくない（成分としてあまり情報を持たない）データを削ぎ落としてしまいます。

例えば、固有値をバープロットしてみて著しく小さな固有値に属する成分をなくしてしまえばいいという事になります。

実際の計算上では、わざわざ全部の成分を保持した変換を行ってから上記のような散布図を獲得するのではなく、予め寄与率のバープロットなどを見てから、いらない成分の固有ベクトルを省いた行列を準備して、必要な潜在変数の成分だけを計算するほうが当然、計算量的に有利です（GPUの時代にそんなことはどうでも良いかもしれないが…）。

確率的主成分分析（最尤推定）

定式化

確率的主成分分析では、潜在変数 $\bf z \in \mathcal R^M$ があり、こちらが線形変換 $\bf W$ を受けて $\bf x \in \mathcal R^D$ になって私達の目に入る（観測する）という問題設定になります。今まで主成分分析ではあくまで観測データ $\bf x$ が主体で、 $\bf z$ をうまく作りたいという設定だったのですが、これが逆転した発想になっていることに注意してください。これは生成モデルと呼ばれるものの基本的な考え方の1つで、「手に入ったデータ $\bf x$ が生まれてくる過程を確率モデルで記述する」という方針になります。

観測をするときには、$ \sigma^2$ の当方的分散をもつガウス分布としてバラツくと考えます。式としては下記のとおりです。

$$ \bf x \sim \mathcal N\bf (Wz, \sigma^2I) $$

潜在変数$\bf z$ はどこからやってくるのかというと、平均 $\bf 0$、分散 $\bf I$の当方的分散を持つガウス分布から生起していると考えます。

$$ \bf z \sim \mathcal N \bf (0, I) $$

解法

この確率変数が $\bf W$ によって変換されて $\bf x = Wz$ となるストーリーです。ここでパラメータ$\bf W$を最尤推定するのが確率的主成分分析になります。ただし、確率的主成分分析では、$\bf z$に関して周辺化を行った確率分布について最尤推定を実施します。$\bf z$ に関する周辺化というのは $ \bf z $ の発生するパターン全てを加味した場合に得られる $\bf x$ を考えるということです。よりフォーマルには、 $\bf p(z)$ に関する期待値を取ります。今、 $\bf x$ の発生過程そのものは下記のようなガウス分布でした。

$$ p({\bf x}\mid {\bf W,z}, \sigma) = \mathcal N\bf (Wz, \sigma^2I) $$

今回は $p({\bf z})=\mathcal N \bf (0, I)$に関する期待値を取るので

$$ p({\bf x} \mid {\bf W}, \sigma) = \int_{\bf z} \mathcal N {\bf (Wz, \sigma^2I)} \mathcal N (0, I) d{\bf z} $$

となります。こいつの計算結果は

$$ p({\bf x} \mid {\bf W}, \sigma) =\mathcal N {\bf (0, WW^T + \sigma^2I)} $$

となります。いきなり天下って申し訳ありませんが、ガウス分布とガウス分布の積は、ガウス分布になり、それぞれの平均と分散の計算も公式として調べれば出てくるので興味があれば検索してみてください（古い手法にガウス分布利用が目立つのは、ガウス分布が現実に現れやすい分布であることもさることながら、計算のしやすさも大きく寄与している）。この分布に関して最尤推定を行えば $\bf W$の解が得られるというわけです。具体的には $D = \{{\bf x}_1,\cdots, {\bf x}_N \}$が実際に手元に入る確率は、上記の $p({\bf x}_n \mid {\bf W}, \sigma)$ が $n= 1, \cdots, N$ に関して起こり続けた場合の確率であり、下記で表されます。

$$ \prod_{n=1}^{N} p({\bf x}_n \mid {\bf W}, \sigma) $$

これを尤度関数と呼び、$\bf W$ について最大化を行えばいいです。このままではかなり扱いにくいので、例のごとく対数尤度を最大化することに問題をすり替えます（対数関数は単調増加かつ、元の関数と極値が同じなので問題ありません）。ちなみに、これは解析的に求まります（ただし、$\bf W$の直交行列による変換の任意性を残すことになりますが）。

ベイズ主成分分析（MAP推定）

ベイズ主成分分析では最尤推定時にパラメータとして点推定した $\bf W$の方も何らかの確率分布位から発生していると考えます。

また、確率変数 $x$ が確率分布がガウス分布 $\mathcal N(\mu, \sigma^2)$から生起することを明記するために $\mathcal N(x\mid \mu, \sigma^2)$ などと表現します。より一般的には $p(x \mid \theta_1, \theta_2, \cdots)$ などと書かれたりします。

さて、ベイズ主成分分析の問題設定の一例は下記のようになります。

$$ \begin{align} \bf x & \sim \mathcal N\bf (x \mid Wz, \sigma^2I) \\ \bf z & \sim \mathcal N \bf (z \mid 0, I)\\ \bf W & \sim \sim \prod_{i=1}^M \left(\frac{\alpha_i}{2\pi} \right)^{D/2} \exp\left(-\frac{1}{2}\alpha_i {\bf w}_i^{\bf T} {\bf w}_i\right) \end{align} $$

ここで$\bf W$ に関しては、各縦ベクトル ${\bf w}_i$ 毎に異なる多次元正規分布を事前分布が当てられていることに注意してください。

補足：事前分布は任意性を持つ

事前分布はモデルの設計者が自ら考えれば良いものであり、任意性があります。例えば下記のように

$$ {\bf W} \sim \prod_i^D \prod_j^M \mathcal N ( w_{ij} \mid 0, \sigma^2) $$

という事前分布を与えても、ベイズ推論を実行することは可能です。ここで $w_{ij}$ とは $\bf W$ の $(i,j)$ 成分です。すなわち上の設定では、「$\bf W$ の各成分が独立に同じガウス分布から生起している」という事前分布を考えたことになります（こういう設定を主成分分析と呼ぶかは別として、ベイズ推論はデータの生成過程を確率モデルで記述し、事前分布を適当に当てるという手順を取ります）。

ベイズ学習の基本

さて、一旦ベイズ学習の基本的な話をします。確率的主成分分析の最尤推定では尤度関数を最大化することで$\bf W$をただ1つ求めていました（そしてそれは解析的に求められるのでした）。ベイズ主成分分析では $\bf W$ をただ1つに求めるのではなく、観測データ $D = \{{\bf x}_1, {\bf x}_2, \cdots, {\bf x}_N\}$ からみた場合には $\bf W$ はこれくらいだろうなということを表す事後分布 $p({\bf W} \mid D)$ を求めます。

これはベイズ主成分分析がというわけではなく、ベイズ学習全体に共通して言えることです。線形回帰でもロジスティック回帰でも、点推定していたパラメータ$\bf W$を確率変数だと考え、そして持ちうる観測データ $D $ から分布を推論します。この事後分布は条件付き分布と同時分布の関係から簡単に下記のように求まります。

$$ p({\bf W} \mid D) = \frac{p(D , {\bf W})}{p(D)} $$

右辺の分子は手元のデータとパラメータの同時分布であり、分母に関しては色々な呼び名がありますが、確率の体裁を整えるための正規化定数だと思っておいてください（汎化誤差の議論をしたい場合には、個々の部分は結構主役級の存在感を示しますが、学習を実行する段階では特に意識しなくても問題ないです）。結局のところベイズ推論を行う上で我々が「意図的に設計する」部分というのは右辺の分子である同時分布になります。「どのようなパラメータ（あるいは潜在変数）が手元の観測データと関係しあっているのだろうか」ということを考えて同時分布を条件付き分布の形で表現するのです。自分自身の一切の考えを導入しないのであれば、

$$ p(D , {\bf W}) = p(D \mid {\bf W})p({\bf W}) $$

という数学的に常に正しい式を入れることになります。

また、事後分布を求める場合において、 これをうまく求めるためには、尤度関数や事前分布を上手に設計する必要（共役事前分布を使うなど）がありましたが、現在はコンピュータの力で近似を実行できるようになってきました。うまく手計算がしやすいようにモデルを設計して行く必要は、今の環境ではそれほど重要ではなく、ちゃんとデータを表すモデルを設計することが大事になります。

またデータ $D$ から「何の」事後分布を求めたいのかは問題によることに注意してください。いつでも観測データを手に入れる手順を確率分布で記述し（生成モデルを記述し）、あとは同時分布と条件付き分布の関係や周辺化などのテクニックで欲しい事後分布を書き下すというのが一連の流れです。

MAP推定

$\bf W$ を確率変数だと考える。そしてその確率分布を真剣に考えてみましょう。と意気込んでみたは良いものの、これは結構難しい問題です。もしもデータが十分に多くあるならば「MAP推定」と呼ばれる点推定の手法で学習を終わらせることもしばしばあります。

そのような場合は

$$ p({\bf W} \mid D) = \frac{p(D , {\bf W})}{p(D)} $$

において右辺の分母は定数のようなものなので

$$ p({\bf W} \mid D) \propto p({\bf W}, D) = p(D \mid {\bf W})p({\bf W})$$

ということにだけ着目して、この事後分布を最大化する $\bf W$ を学習して終わるということを行います。 これは事後分布をちゃんと求めてから、事後確率が最大になるような$\bf W$を代表して予測モデルとして使うということと等価です。事後分布はデータが十分に多いときには、（単峰性なら）ある1つの値に収束していくので、別に最初から分布を求めることを諦めてMAP推定でも良いやと考えているわけです。

右辺は「尤度」と「事前分布」から構成されているので、対数を取って最大化することを考えると

$$ \log p(D \mid {\bf W}) + \log p({\bf W}) $$

という式を最大化することになり、最尤推定と比べて $\log(\bf W)$ が余分に増えています。この項は通常の機械学習手法て呼ばれるところの正則化項に相当します。例えば事前分布がガウス分布に設定されているのであれば、L2正則化と同じ機能を果たします。一方で無情報事前分布を用いる場合には最尤推定と同じ結果を示します。これは単なる最適化問題ですのでTensorFlowの自動微分機能を使えば簡単に実装が出来ます。例えばコードの一部は下記のようになるでしょう。

主成分分析のMAP推定コード（TensorFlow Probability）

インポートとデータの準備

とりあえずデータ数は200個で、次元は5、潜在変数の次元は3ということにしておきます。

import tensorflow as tf
import tensorflow_probability as tfp
tfe = tf.contrib.eager
tfd = tfp.distributions

def toy_model(N, data_dim, latent_dim, scale):
    z_rv = tfd.Normal(loc=tf.zeros([latent_dim]),
                      scale=tf.ones([latent_dim])).sample(200)
    W_rv = tfd.Normal(loc=tf.zeros([latent_dim, data_dim]),
                      scale=tf.ones([latent_dim, data_dim])).sample()
    x_rv = tfd.Normal(
        loc=tf.matmul(z_rv, W_rv),
        scale=scale
    ).sample()
    return x_rv

N = 200
data_dim = 5
latend_dim = 3


x_train = toy_model(200, 5, 3)

モデル

MAP推定は実質「単なる最適化問題」なので普通のTensorFlowの要領で記述できます。 今回は、パラメータ$\bf W$ と 潜在変数$\bf z$ が推定したい変数なのでこいつを tf.Variable で学習対象として初期化しておきます。その後は、確率変数$\bf x, z, W$ の同時分布を書き下します。$\propto p({\bf W, z, x}) = p({\bf x \mid W, z})p(\bf z)p(\bf W)$ としておきます。（最右辺は、$\bf z$ と $\bf W$ を個別に事前分布を置いたということです。これが、いやいや、両者は関係しあっている！と思うなら、両者をサンプリングするような適当な事前分布を選択しましょう）

z_var = tf.Variable(tf.random_normal([N, latent_dim]))
W_var = tf.Variable(tf.random_normal([latent_dim, data_dim]))

def joint_log_prob(x, z, W):
    z_rv = tfd.Normal(loc=tf.zeros([N, latent_dim]),
                      scale=tf.ones([N, latent_dim]))
    
    W_rv = tfd.Normal(loc=tf.zeros([latent_dim, data_dim]),
                      scale=tf.ones([latent_dim, data_dim]))
    x_rv = tfd.Normal(
        loc=tf.matmul(z, W),
        scale=0.5
    )
    return (
        tf.reduce_mean(tf.reduce_sum(x_rv.log_prob(x), axis=1)) +
        tf.reduce_mean(tf.reduce_sum(z_rv.log_prob(z), axis=1))+
        tf.reduce_sum(W_rv.log_prob(W))
    )

loss = lambda z, W: - joint_log_prob(x_train, z, W)

最後の loss = lambda z, W: - joint_log_prob(x_train, z, W)の部分は、同時分布の中の $\bf x$ に関しては既に手元の x_train で確定しているので先に代入しておきましょうということです。負を付けているのはTensorFlowが通常は最小化問題を扱うので、そのために符号を入れ替えているだけで本質的ではありません。あとは $\bf z, W$の実現値を代入し、loss を計算し、この lossが下がる方へ $\bf z, W$を更新することを繰り返せば、 同時確率が最大（事後確率が最大）になる $\bf W, z$が求まります。

最適化ループ

普通に学習するだけです。

optimizer = tf.train.AdamOptimizer()

grads_and_value = tfe.implicit_value_and_gradients(loss)

for i in range(1000):
    loss, grads = grads_and_value(z_var, W_var)
    optimizer.apply_gradients(grads)
    if i % 100 == 0:
        print("loss{}: {}".format(i, loss.numpy()))

学習されたモデルからのデータの生成

学習されたモデルの $\bf W$ を用いて、適当に正規分布から発生させた潜在変数を変換して、観測データを生成してみましょう。 まずは学習に用いたデータの可視化です。

import pandas as pd
import seaborn as sns

sns.pairplot(pd.DataFrame(x_train.numpy()))

続いてデータを生成して100サンプルだけプロットしてみます。

x_generated = tf.matmul(tf.random_normal([100, latent_dim]), W_var)
sns.pairplot(pd.DataFrame(x_generated.numpy()))

概ね同じようなデータが生成されましたね。

下記の項目については記事を分けて書きます。