Probabilistic Principal Component Analysis

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 +ϵ$$
其中$t\in {\mathbb{R}}^d$为观测变量，也就是样本，而$x\in {\mathbb{R}}^q$是隐变量，$W\in {\mathbb{R}}^{d\times q}$将$x,t$二者联系起来。另外，$ϵ$是噪声。

令$S=\frac{1}{N}\sum _{n=1}^N(t_n-\mu )(t_n-\mu {)}^T$是样本协方差矩阵，其中$\mu$是样本均值。论文的主要工作，就是将$S$的列空间和$W$联系起来。

主要内容

假设$ϵ\sim N(0,{\sigma }^{2}I)$,$x\sim N(0,I)$，二者独立。那么，容易知道$t$在$x$已知的情况下的条件概率为：
$$t|x\sim N(Wx+\mu ,{\sigma }^{2}I)$$
然后论文指出，通过其可求得$t$的边际分布：
$$t\sim N(\mu ,C)$$
其中$C=W{W}^T+{\sigma }^{2}I$。这个证明，在贝叶斯优化中有提过，不过我发现，因为$t=Wx+\mu +ϵ$，是服从正态分布随机变量的线性组合，所以$t$也服从正态分布，所以通过$E(t)$和$E((t-E(t))(t-E(t){)}^T)$也可以得到$t$的分布。

其似然函数$L$为：

将$W，\sigma$视为参数，我们可以得到其极大似然估计：

其中$U_q$的列是$S$的主特征向量，而${\Lambda }_q$的对角线元素为特征向量对应的特征值${\lambda }_1,\dots ,{\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^{T}W+{\sigma }^{2}I$
这个推导需要利用贝叶斯公式：
$$p(x|t)=\frac{p(t|x)p(t)}{p(t)}$$

为什么要提及这个东西，因为可以引出一个很有趣的性质，注意到$x|t$的均值为：
$$M^{-1}{W}^T(t-u)$$
令$W=W_{ML}$,且假设${\sigma }^2\to 0$，那么均值就成为：
$$(W_{ML}^T{W}_{ML}{)}^{-1}{W}_{ML}^T(t-u)$$
实际上就是$(t-u)$在主成分载荷向量上的正交投影，当然这里不要计较$W_{ML}^T{W}_{ML}$是否可逆。这就又将PPCA与普通的PCA联系在了一起。

EM算法求解

论文给出了$W$的显式解（虽然有点地方不是很明白），也给出了如何利用EM算法来求解。

构造似然估计：

对$x_n$求条件期望（条件概率为$p(x_n|t_n,W,{\sigma }^2)$）：

$M$步是对上述$W,\sigma$求极大值，注意$<\cdot >$里面的$M,\sigma$是已知的（实际上，用$M^{\prime },{\sigma }^{\prime }$来表述更加合适）：

有更加简练的表述形式：

符号虽然多，但是推导并不麻烦，自己推导的时候并没有花多大工夫。

附录

极大似然估计

已知对数似然函数为：

先考察对$W$的微分：
d L = − N 2 { d ∣ C ∣ ∣ C ∣ + d t r ( C − 1 S ) } \begin{array}{ll} \mathrm{d}L = -\frac{N}{2}\{\frac{\mathrm{d}|C|}{|C|} + \mathrm{dtr}(C^{-1}S)\} \end{array} dL=−2N​{∣C∣d∣C∣​+dtr(C−1S)}​
$$\begin{array}{ll}\frac{\mathrm{d}|C|}{|C|}& =\mathrm{t}\mathrm{r}(C^{-1}\mathrm{d}C)\\ & =\mathrm{t}\mathrm{r}[C^{-1}(\mathrm{d}W{W}^T+W\mathrm{d}{W}^T)]\\ & =2\mathrm{t}\mathrm{r}[W^T{C}^{-1}\mathrm{d}W]\end{array}$$
$$\begin{array}{ll}\mathrm{d}\mathrm{t}\mathrm{r}(C^{-1}S)& =\mathrm{t}\mathrm{r}(\mathrm{d}{C}^{-1}S)\\ & =-\mathrm{t}\mathrm{r}(C^{-1}[\mathrm{d}C]C^{-1}S)\\ & =-\mathrm{t}\mathrm{r}(C^{-1}S{C}^{-1}\mathrm{d}C)\\ & =-2\mathrm{t}\mathrm{r}(W^T{C}^{-1}S{C}^{-1}\mathrm{d}W)\end{array}$$
所以，要想取得极值，需满足：
$$\begin{gathered} C^{-1}W=C^{-1}S{C}^{-1}W \\ \Rightarrow S{C}^{-1}W=W \end{gathered}$$

论文说这个方程有三种解：

1. $W=0$，0解，此时对数似然函数取得最小值（虽然我没看出来）。

2. $C=S$:
   $$\begin{gathered} W{W}^T+{\sigma }^{2}I=S \\ \Rightarrow W{W}^T=S-{\sigma }^2 \end{gathered}$$
   其解为：
   $$W=U_S({\Lambda }_S-{\sigma }^{2}I{)}^{1/2}R$$
   其中$S=U_S\Lambda U_S^T$。

3. 第三种，也是最有趣的解，$S{C}^{-1}W=W$但是$W̸=0,C̸=S$。假设$W=UL{V}^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 }^{2}L+L^2)L$$
   于是$S{u}_j=({\sigma }^{2}I+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 }^{2}I{)}^{1/2}R$$
   其中：
   $$k_j=\{\begin{array}{ll}{\lambda }_j& 对应特征值u_j\\ {\sigma }^2& \end{array}$$
   实际上就是$k_j={\lambda }_j$

注意，上面的分析只能说明其为驻定解，事实上$U_q$只说明了其为$S$的特征向量，而没有限定是哪些特征向量。

将解$W=U_q(K_q-{\sigma }^{2}I{)}^{1/2}R$代入对数似然函数可得（$C=W{W}^T+{\sigma }^{2}I$）：

其中$q^{\prime }$是非零$l_1,\dots ,l_{q^{\prime }}$的个数。
上面的是蛮好证明的，注意 { ⋅ } \{\cdot\} {⋅}中第2项和第4项和为$\ln |C|$，第3，5项构成$\mathrm{t}\mathrm{r}(C^{-1}S)$。
对$\sigma$求极值，可得：

且是极大值，因为显然$\sigma \to 0$会导致$L\to -\infty$。代入原式可得：

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

注意到，上式只与被舍弃的$l_j=0$的${\lambda }_j$有关，又${\lambda }_i\ge {\sigma }^2,i=1,\dots ,q$，再结合(18)，可以知道最小的特征值一定是被舍弃的。但是论文说，应当是最小的$d-q^{\prime }$个特征值作为被舍弃的（因为这些特征值必须在一块？）。

仔细想来，似然函数可以写成：
L = − N 2 { d ln ⁡ ( 2 π ) + ∑ j = 1 q ln ⁡ ( λ j ) + 1 σ 2 ∑ j = q + 1 d λ j + ( d − q ) ln ⁡ ( σ 2 ) + 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\} L=−2N​{dln(2π)+j=1∑q​ln(λj​)+σ21​j=q+1∑d​λj​+(d−q)ln(σ2)+q}

好吧，还是不知道该如何证明。

代码

"""
瞎写的，测试结果很诡异啊
"""
import numpy as np

class PPCA:

    def __init__(self, data, q):
        self.__data = np.array(data, dtype=float)
        self.__n, self.__p = data.shape
        self.__mean = np.mean(self.data, 0)
        self.q = q
        assert q < self.__p, "Invalid q"
    @property
    def data(self):
        return self.__data

    @property
    def n(self):
        return self.__n

    @property
    def p(self):
        return self.__p

    def explicit(self):
        data = self.data - self.__mean
        S = data.T @ data / self.n
        value, vector = np.linalg.eig(S)
        U = vector[:, :self.q]
        sigma = np.mean(value[self.q:])
        newvalue = value[:self.q] - sigma

        return U * newvalue

    def EM(self):
        data = self.data - self.__mean
        S = data.T @ data / self.n
        W_old = np.random.randn(self.p, self.q)
        sigma = np.random.randn()
        count = 0
        while True:
            count += 1
            M = W_old.T @ W_old + sigma
            M_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.p

            if 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_new
            else:
                W_old = W_new
                sigma = sigma_new