# Probabilistic PCA

> Organization contact [name= [ierosodin]([e-mail])]
###### tags: `deep learning` `學習筆記`

==[Back to Catalog](https://hackmd.io/@ierosodin/Deep_Learning)==

* PPCA 優點：
* 可以推導出EM算法進行迭代求解。
* 在主成份M定得小的情況下非常高效。
* 把PCA概率模型化結合EM算法可以處理數據缺失的情況。
* 多個PPCA模型可以混合，並用EM算法進行訓練。
* 子空間的維度可以自動選擇。建立了似然函數，便於模型比較。
* PPCA建立瞭如何由隱變量產生觀測變量的過程，可以進行採樣產生樣本數據。
* PPCA can be used to model class-conditional densities and hence be applied to classification problems.
* PPCA represents a constrained form of the Gaussian distribution in which the number of free parameters can be restricted while still allowing the model to capture thee dominant correlations in a data set.
* 假設
* $Z$：隱變量，代表低維子空間上的點。
* $P(z) = N(z|0,I)$
* 隱變量的分佈為高斯，協方差矩陣為單位陣。
* 假設觀測數據的產生過程為：$x = Wz + \mu + \epsilon$
* 其中$W$ 是$M$ 維到$D$ 維的投影矩陣； $\mu$ 使得觀測數據可以為非零均值； $\epsilon$ 是一個isotropic 的高斯噪聲$N(0, {\sigma} ^{2}I)$。
* 所以 $p(x|z)=N(x|Wz+\mu, {\sigma}^{2}I)$
* $z$ 和 $\epsilon$ 相互獨立。
* 現在已知 $p(x|z)$ 和 $p(z)$，可以求出 $p(x)$。
* $p(x) = \int p(x|z)(p(z)dz$
* $p(x)$ 仍然是高斯，$p(x) = N(x; \mu, C)$
* 其中
* $C = WW^T + \sigma I$
* 因為
* $E[x] = E[Wz + \mu + \epsilon] = \mu$
* $\begin{split}cov[x] &= E[(Wz + \epsilon)(Wz + \epsilon)^T] \\
&= E[Wzz^TW^T] + E[\epsilon \epsilon^T] = WW^T + \sigma^2I\end{split}$
* 已知 $p(z)$ 和 $p(x|z)$，根據 Baysian 公式可以算出 $p(z|x)$。
* $p(x, z) = p(z)p(x|z)$, so that $p(z|x) \propto p(z)p(x|z)$
$\begin{split}p(z|x) \sim p(z)p(x|z) &= cexp\{-\frac{1}{2}z^Tz - \frac{1}{2\sigma^2}{||Wz + \mu - x||}^2\} \\
&\sim exp{\{\sigma^{-2}{||Wz + \mu - x||}^2 + z^Tz\}} \\
&= exp\{\sigma^{-2}\lgroup (Wz + \mu)^T(Wz + \mu - 2(Wz + \mu)^Tx + x^Tx \rgroup + z^Tz\} \\
& \begin{split} =exp\{&\sigma^{-2}\lgroup z^TW^TWz + 2z^TW^T\mu + \mu^T\mu - 2z^TW^Tx - 2\mu^Tx + x^Tx \rgroup + \\
&z^Tz\} \end{split} \\
&= exp\{z^T(\sigma^{-2}W^TW + I)z -2z^T\lgroup \sigma^{-2}W^T(x - \mu) \rgroup + const\}
\end{split}$
* We could find that $p(z|x)$ is still in gaussian form.
* Compare exponatial term with gaussian quadratic form:
$(z - \bar \mu)^T\Sigma^{-1}(z - \bar \mu) = z^T\Sigma^{-1}z - 2z^T\Sigma^{-1}\bar \mu + \bar \mu^T \bar \mu$
$\therefore \Sigma^{-1} = \sigma^{-2}W^TW + I = \sigma^{-2}(W^TW + \sigma^{2}I)$
Let $M = W^TW + \sigma^{2}I \Rightarrow \Sigma = \sigma^{2}M^{-1}$
$\bar \mu = \Sigma\lgroup \sigma^{-2}W^T(x - \mu) \rgroup = M^{-1}W^T(x - \mu)$
$\Rightarrow p(z|x) = N(z; M^{-1}W^T(x - \mu), \sigma^{2}M^{-1})$, where $M = W^TW + \sigma^{2}I$
* Maximum likelihood for PCA
* $x$ 的 log likelihood function 可以寫作
* 
* 對 $\mu$ 求導為零，容易得 $\mu= \bar{x}$ ，帶入得
* 
* 對 $W$ 和 ${\sigma}^{2}$ 求導就不那麼容易了（推導見最後）。 Tipping&Bishop（1999b）給出了一個閉式解。
* 
* 其中 ${U}_{M}$ 是 $D*M$ 的矩陣，其列向量為協方差矩陣 $S$ 的的任意 $M$ 個特徵向量。 ${L}_{M}$ 是對角陣，元素值為對應的特徵值。 $R$ 是酉矩陣，理解為對投影矩陣進行旋轉。
* 對 $W_{ML}$ 的解釋：$R$ 是個酉矩陣，所以 $RR’=I$。所以不管 $R$ 怎麼變，$WW'$ 總不變，也就是說 $x$ 的分佈不變，也就是說不管隱變量對應的子空間怎麼旋轉，只要經過 $x = Wz + \mu + \epsilon$ 那樣的線性變換，得到的觀測變量 $x$ 的分佈都不變。在 $R=I$ 的情況下，$W$ 的列向量可以看成是 ${U}_{M}$ 的列向量被 $({\lambda}_{i}-{\sigma}^{2})$ 放縮的形式，考慮到原數據 $x$ 的產生式 $x = Wz + \mu + \epsilon$，這就比較好理解了，${\lambda}_{i}$ 是原數據在 ${u}_{i}$ 方向上的方差，而 ${\sigma}^{2}$ 是 $x = Wz + \mu + \epsilon$ 中加法模型的誤差的方差，這樣由隱空間產生觀測變量空間的產生過程就與實際觀測到的數據 ${{x}_{n}}$ 就對應上了。
* $\sigma_{ML}$ 的解釋也比較簡單，誤差的方差是被丟掉的那些坐標 $({u}_{M+1} … {u}_{D})$ 方向上的方差的平均。
* 實際操作時，有了觀測數據，也就是數據的協方差矩陣 $S$，如何估計出最可能的產生模型 $x = Wz + \mu + \epsilon$？首先用 $W_{ML}$ 算 ${\sigma}^{2}$ ，再使 $R=I$ 算 $W$。也可以用梯度法來對似然函數進行優化，但優化的結果對應的 $R$ 就不一定是單位矩陣了。
* 由 $p(z|x) = N(z; M^{-1}W^T(x - \mu), \sigma^{2}M^{-1})$ 看出，
* 
* 當 ${\sigma}^{2}->0$ 時，$z$ 的後驗均值：
* 
* 這就與 conventional PCA 聯繫起來了。
* EM algorithm for PCA
* 既然每一個觀測變量 $x$ 對應一個因變量 $z$，自然而然想到 EM 算法來算參數的極大似然估計。要估計的參數包括 ${\mu，W和{\sigma}^{2}}$，這裡的EM算法並不是估計所有的參數，而是只估計 $W$ 和 ${\sigma}^{2}$ ，而 $\mu$ 的求解在上一節給出。因為 $W$ 和 ${\sigma}^{2}$ 的閉式解涉及到對S進行特徵值分解，用 EM 算法來迭代求解有優勢，而 $\mu$ 僅僅是觀測數據的均值，沒什麼計算量。
* complete-data 的 log likelihood
* 
* 下面計算 complete-data 的 log likelihood 關於 $z$ given $x$ 的條件期望。並把 $\mu$ 的估計值帶入，得
* 
* 由可以得到 $p(z|x) = N(z; M^{-1}W^T(x - \mu), \sigma^{2}M^{-1})$
* 
* 均為z的後驗期望，參數均為上一步的。接下來 M 步就是最大化 log likelihood w.r.t. $W$ 和 ${\sigma}^{2}$ ，得到
* 
* 如果 ${\sigma}^{2} = 0$ 即沒有噪聲時，我們發現就變成了 conventional PCA 只不過求解方式不再是特徵值分解而是迭代。 E 步只需要算 $z$ 關於 $z$ given $x$ 的條件期望，M 步只需要更新 $W$。