目錄

• 引
• 主要內容
  • EM算法求解
• 附錄
  • 極大似然估計
• 代碼

Tipping M E, Bishop C M. Probabilistic Principal Component Analysis[J]. Journal of The Royal Statistical Society Series B-statistical Methodology, 1999, 61(3): 611-622.

引

PPCA 通過高斯過程給出了普通PCA一個概率解釋，這是很有意義的。論文還利用PPCA進行缺失數據等方面的處理（不過這方面主要歸功于高斯過程吧）。

\[ t = Wx + \mu + \epsilon \]
其中\(t \in \mathbb{R}^d\)為觀測變量，也就是樣本，而\(x \in \mathbb{R}^q\)是隱變量，\(W \in \mathbb{R}^{d \times q}\)將\(x,t\)二者聯系起來。另外，\(\epsilon\)是噪聲。

令\(S = \frac{1}{N} \sum \limits_{n=1}^N (t_n -\mu )(t_n - \mu)^T\)是樣本協方差矩陣，其中\(\mu\)是樣本均值。論文的主要工作，就是將\(S\)的列空間和\(W\)聯系起來。

主要內容

假設\(\epsilon \sim N(0, \sigma^2 I)\),\(\: x \sim N(0, I)\)，二者獨立。那么，容易知道\(t\)在\(x\)已知的情況下的條件概率為：
\[ t|x \sim N(Wx + \mu, \sigma^2I) \]
然后論文指出，通過其可求得\(t\)的邊際分布：
\[ t \sim N(\mu, C) \]
其中\(C = WW^T + \sigma^2 I\)。這個證明，在貝葉斯優化中有提過，不過我發現，因為\(t=Wx+\mu + \epsilon\)，是服從正態分布隨機變量的線性組合，所以\(t\)也服從正態分布，所以通過\(E(t)\)和\(E((t-E(t))(t-E(t))^T)\)也可以得到\(t\)的分布。

其似然函數\(L\)為：

將\(W，\sigma\)視為參數，我們可以得到其極大似然估計：

其中\(U_{q}\)的列是\(S\)的主特征向量，而\(\Lambda_q\)的對角線元素為特征向量對應的特征值\(\lambda_1, \ldots, \lambda_q\)（為所有特征值的前\(q\)個，否則\(W\)將成為鞍點），\(R \in \mathbb{R}^{q \times q}\)是一個旋轉矩陣。注意到，\(W_{ML}\)的列向量并不一定正交。

這部分的推導見附錄。

同樣的，我們可以推導出，\(x\)在\(t\)已知的情況下的條件分布：
\[ x|t \sim N(M^{-1}W^T(t-\mu),\sigma^2 M^{-1} \]
其中\(M = W^TW+\sigma^2I\)
這個推導需要利用貝葉斯公式：
\[ p(x|t) = \frac{p(t|x)p(t)}{p(t)} \]

為什么要提及這個東西，因為可以引出一個很有趣的性質，注意到\(x|t\)的均值為：
\[ M^{-1}W^T(t-u) \]
令\(W = W_{ML}\),且假設\(\sigma^2 \rightarrow 0\)，那么均值就成為：
\[ (W_{ML}^TW_{ML})^{-1}W_{ML}^T(t-u) \]
實際上就是\((t-u)\)在主成分載荷向量上的正交投影，當然這里不要計較\(W_{ML}^TW_{ML}\)是否可逆。這就又將PPCA與普通的PCA聯系在了一起。

EM算法求解

論文給出了\(W\)的顯式解（雖然有點地方不是很明白），也給出了如何利用EM算法來求解。

構造似然估計：

對\(x_n\)求條件期望（條件概率為\(p(x_n|t_n,W,\sigma^2)\)）：

\(M\)步是對上述\(W,\sigma\)求極大值，注意\(<\cdot>\)里面的\(M, \sigma\)是已知的（實際上，用\(M', \sigma'\)來表述更加合適）：

有更加簡練的表述形式：

符號雖然多，但是推導并不麻煩，自己推導的時候并沒有花多大工夫。

附錄

極大似然估計

已知對數似然函數為：

先考察對\(W\)的微分：
\[ \begin{array}{ll} \mathrmze8trgl8bvbqL = -\frac{N}{2}\{\frac{\mathrmze8trgl8bvbq|C|}{|C|} + \mathrm{dtr}(C^{-1}S)\} \end{array} \]
\[ \begin{array}{ll} \frac{\mathrmze8trgl8bvbq|C|}{|C|} &= \mathrm{tr}(C^{-1}\mathrmze8trgl8bvbqC) \\ &= \mathrm{tr}[C^{-1}(\mathrmze8trgl8bvbqWW^T+W\mathrmze8trgl8bvbqW^T)] \\ &= 2\mathrm{tr}[W^TC^{-1}\mathrmze8trgl8bvbqW] \\ \end{array} \]
\[ \begin{array}{ll} \mathrm{dtr}(C^{-1}S) &= \mathrm{tr}(\mathrmze8trgl8bvbqC^{-1}S) \\ &= -\mathrm{tr}(C^{-1}[\mathrmze8trgl8bvbqC]C^{-1}S) \\ &= -\mathrm{tr}(C^{-1}SC^{-1}\mathrmze8trgl8bvbqC) \\ &= -2\mathrm{tr}(W^TC^{-1}SC^{-1}\mathrmze8trgl8bvbqW) \\ \end{array} \]
所以，要想取得極值，需滿足：
\[ C^{-1}W = C^{-1}SC^{-1}W \\ \Rightarrow \quad SC^{-1}W=W \]

論文說這個方程有三種解：

\(W=0\)，0解，此時對數似然函數取得最小值（雖然我沒看出來）。

\(C=S\):
\[ WW^T + \sigma^2 I = S \\ \Rightarrow WW^T = S-\sigma^2 \]
其解為：
\[ W = U_{S} (\Lambda_S-\sigma^2I)^{1/2}R \]
其中\(S = U_S \Lambda U_S^T\)。

第三種，也是最有趣的解，\(SC^{-1}W=W\)但是\(W \ne 0, C \ne S\)。假設\(W=ULV^T\),其中\(U \in \mathbb{R}^{d \times q}\), \(L \in \mathbb{R}^{q \times q}\)為對角矩陣，\(V \in \mathbb{R}^{q \times q}\)。通過一系列的變換(我沒有把握能完成這部分證明，感覺好像是對的)，可以得到：
\[ SUL = U(\sigma^2L+L^2)L \]
于是\(Su_j = (\sigma^2I + l_j^2)u_j\)，其中\(u_j\)為\(U\)的第j列，\(l_j\)為\(L\)的第j個對角線元素。因此，\(u_j\)就成了\(S\)的對應特征值\(\lambda_j = \sigma^2 + l_j^2\)的特征向量（注意到這個特征值是必須大于等于\(\sigma^2\)）。于是，有：
\[ W = U_q(K_q-\sigma^2I)^{1/2}R \]
其中：
\[ k_j = \left \{ \begin{array}{ll} \lambda_j & 對應特征值u_j \\ \sigma^2 \end{array} \right . \]
實際上就是\(k_j=\lambda_j\)

注意，上面的分析只能說明其為駐定解，事實上\(U_q\)只說明了其為\(S\)的特征向量，而沒有限定是哪些特征向量。

將解\(W = U_q(K_q-\sigma^2I)^{1/2}R\)代入對數似然函數可得（\(C = WW^T+\sigma^2 I\)）：

其中\(q'\)是非零\(l_1,\ldots,l_{q'}\)的個數。
上面的是蠻好證明的，注意\(\{\cdot\}\)中第2項和第4項和為\(\ln |C|\)，第3，5項構成\(\mathrm{tr}(C^{-1}S)\)。
對\(\sigma\)求極值，可得：

且是極大值，因為顯然\(\sigma \rightarrow 0\)會導致\(L \rightarrow - \infty\)。代入原式可得：

最大化上式等價于最小化下式：

注意到，上式只與被舍棄的\(l_j=0\)的\(\lambda_j\)有關，又\(\lambda_i \ge \sigma^2, i=1,\ldots, q\)，再結合(18)，可以知道最小的特征值一定是被舍棄的。但是論文說，應當是最小的\(d-q'\)個特征值作為被舍棄的（因為這些特征值必須在一塊？）。

仔細想來，似然函數可以寫成：
\[ L = -\frac{N}{2} \{d \ln (2\pi) + \sum \limits_{j=1}^q \ln (\lambda_j) + \frac{1}{\sigma^2}\sum \limits_{j=q+1}^d \lambda_j + (d-q)\ln (\sigma^2) + q\} \]

好吧，還是不知道該如何證明。

代碼

""" 瞎寫的，測試結果很詭異啊 """ import numpy as npclass PPCA:def __init__(self, data, q):self.__data = np.array(data, dtype=float)self.__n, self.__p = data.shapeself.__mean = np.mean(self.data, 0)self.q = qassert q < self.__p, "Invalid q"@propertydef data(self):return self.__data@propertydef n(self):return self.__n@propertydef p(self):return self.__pdef explicit(self):data = self.data - self.__meanS = data.T @ data / self.nvalue, vector = np.linalg.eig(S)U = vector[:, :self.q]sigma = np.mean(value[self.q:])newvalue = value[:self.q] - sigmareturn U * newvaluedef EM(self):data = self.data - self.__meanS = data.T @ data / self.nW_old = np.random.randn(self.p, self.q)sigma = np.random.randn()count = 0while True:count += 1M = W_old.T @ W_old + sigmaM_inv = np.linalg.inv(M)W_new = S @ W_old @ np.linalg.inv(sigma + M_inv @ W_old.T @ S @ W_old)sigma_new = np.trace(S - S @ W_old @ M_inv @ W_new.T) / self.pif np.sum(np.abs(W_new - W_old)) / np.sum(np.abs(W_old)) < 1e-13 and \np.abs(sigma_new - sigma) < 1e-13:return W_newelse:W_old = W_newsigma = sigma_new

轉載于:https://www.cnblogs.com/MTandHJ/p/10795355.html

總結

以上是生活随笔為你收集整理的Probabilistic Principal Component Analysis的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。