主成分分析PCA
Published:
主成分分析PCA相关内容，从最简单的PCA，再到核PCA、概率PCA等。
给定数据$X$，由于维数灾难（Curse of Dimensionality）的存在，对于高维数据的处理存在很多难题，此时通常需要对减低数据的维数，也即降维，而PCA就是经典的降维方式。
PCA
最简单的假设是数据服从正态分布，$X\sim \mathcal{N}(\mu,\Sigma)$
考虑对方差矩阵做特征值分解，并且取前$k$个特征值和其所对应的特征向量，可以看作忽略掉小方差所对应的变量，由于我们常常认为小方差对应的随机变量有更大的可能为噪声，因此此时可以取到去噪或者降维的作用。此时特征向量构成的矩阵$Q$，可以看作将数据$X$进行线性变换，变换到一组独立的随机变量，也即，$QX \sim \mathcal{N}(Q \mu, Q \Lambda Q^T)$
由于实际的方差矩阵$\Sigma$是未知的，可以考虑用样本方差矩阵$S = \frac{1}{n} X^T H X,H = I - \frac{1}{n}ee^n$ 代替，此时相当于对 $X^T HX$ 进行了特征值分解，也即对 $HX$ 进行了奇异值分解。
对于不服从正态分布的数据，PCA的合理性在于可以将此时的降维任务看作一个优化问题：最佳低秩逼近问题，也即寻找到一个低秩的矩阵逼近原本的高秩的数据矩阵，对于该部分的理论证明可以移步至 特征值不等式与最佳低秩逼近
也可以将PCA看作是一个自编码器（Auto-Encoder），寻找一个矩阵$Q$对数据$X$进行变换（编码，Encode），利用$Q^T$对变换后的数据进行恢复（解码，Decode），而最佳低秩逼近理论告诉我们如果使用F范数度量重构误差，PCA是使得重构误差最小化的解。
PCoA
主成分分析（PCA，Principal Component Analysis）基于数据协方差矩阵 $X^TX$ 对数据进行降维，而主坐标分析（PCoA，Principal Coordinate Analysis）基于数据距离矩阵 $XX^T$ 对数据进行降维，从奇异值分解的角度两者都是一样的，从下式可以看出，
\[\begin{align} X & = U \Sigma V^T \\ X^T X &= V \Sigma^2 V^T \\ X X^T &= U \Sigma^2 U^T \end{align}\]但在数据计算上，PCoA和PCA计算的矩阵大小不一样，对于样本数据为图像、基因这种高维数据（图像数据的特征维数为所有像素点的个数，基因的数据维数为所有基因的个数），计算距离矩阵 $XX^T$是更为方便的，且PCoA的方法也可以很好地推广到核方法。
Kernel PCA
核方法是一种很重要的技巧，核方法的洞见在于在低维空间线性不可分的数据，在高维空间内有更高的概率线性可分，因此很多问题在高维空间内更高处理。同时，一个低维到高维空间的映射也在模型中引入了非线性因素，使模型具有更强的表达能力。
Kernel Method
问题是低维空间到高维空间的映射函数$\phi$通常比较复杂，显示地将其表达非常困难，但核方法的关键是仅需要知道$\phi$的内积表示矩阵$K$，称为核（Kernel），通常就可以解决大部分的问题，而直接显示地表示$K$是容易得多的，形式化地，$K$有如下定义：
\[\begin{align} K(x,y) &= \langle \varphi(x), \varphi(y) \rangle ,\text{Kernel Function}\\ K_{ij} &= \langle \varphi(X_i), \varphi(X_j) \rangle ,\text{Kernel Matrix}\\ \end{align}\]$K$定义了一个关于内积的函数，称为核函数（Kernel Function），而对于样本$X$，可以计算得到核矩阵$K_{ij}$，表示样本$X_i$之间内积的结果。 对$K$有一定的限制，通常要求其为对称半正定矩阵，称为正定核，本质上在要求$K$是一个Hermit算子，根据泛函分析的理论，Hermit算子一定可以进行谱分解，其中$\phi_i(x)$为正交基：
\[K(x,y) = \sum_{i=1}^{\infty} \gamma_{i} \varphi_i(x) \varphi_i(y)\]该谱分解实际上说明了，存在一个映射$\varphi$，此时可以将核函数看作$\varphi$映射后空间上的内积，
\[K(x,y) = \langle\varphi(x), \varphi(y) \rangle\]在此基础上，可以定义函数的内积和范数，
\[\begin{align} \langle g,g' \rangle &= \langle \alpha_i \varphi_i(x), \alpha_i' \varphi_i(x) \rangle = \sum_{i=1}^{\infty} \frac{\alpha_i \alpha_i'}{\gamma_{i}} \\ \Vert g \Vert &= \langle \alpha_i \varphi_i(x) , \alpha_i \varphi_i(x) \rangle = \sum_{i=1}^{\infty} \frac{\alpha_i^2}{\gamma_{i}} \end{align}\]据此，可以得到$K$的如下性质：
性质1 $\langle K(x,y), g(x) \rangle = g(y)$
将核函数$K$和函数$g$都在正交基下展开，
\[\begin{align} \langle K(x,y), g(x) \rangle &= \langle \sum_{i=1}^{\infty} \gamma_{i} \varphi_i(x) \varphi_i(y), \sum_{i=1}^{\infty} \alpha_i \varphi_i(x) \rangle \\ &= \sum_{i=1}^{\infty} \alpha_i \phi_i(y) \\ &=g(y) \end{align}\]性质2 $\langle K(x’,x), K(x’,y) \rangle = K(x,y)$
证明只需在性质1的基础上，代入$g(x’) = K(x’,y)$ 即可。
根据该性质，可以基于样本数据$X_i$定义如下函数，
\[g(x) = \sum_{i=1}^N \alpha_i K( x, X_i )\]该函数在核岭回归（Kernel Ridge Regression），支持向量机（SVM，Support Vector Machine) 等都扮演着重要的角色。
性质3 $\Vert g \Vert = \sum_{ij} \alpha_i \alpha_j K(X_i,X_j) $
证明用到性质2及内积的线性性质，
\[\begin{align} \Vert g \Vert &= \langle \sum_{i}\alpha_i K(x,X_j) , \sum_{i}\alpha_i K(x,X_j) \rangle \\ & =\sum_{ij} \alpha_i \alpha_j \langle K(x,X_i) , K(x,X_j) \rangle \\ &= \sum_{ij} \alpha_i \alpha_j K(X_i,X_j) \end{align}\]给出了核函数的一些性质后，我们可以发现常见的核函数，如多项式核函数、高斯核函数等，都满足对称正定性质。而在现有的核函数的基础上，我们可以得到更多的核函数，需要用到矩阵的如下性质：
性质1 若$K_1,K_2$为正定阵，则$K_1+K_2$也为正定阵，也即加法运算保持正定性质，
\[\begin{align} x^T (K_1 +K_2) x & = x^T K_1 X + X^T K_2 x \ge 0 \end{align}\]性质2 若$K_1,K_2$为正定阵，则$K_1 \otimes K_2$也为正定阵 , 也即Kronecker积保持正定性质，证明用到了正定阵的Cholesky 分解，
\[\begin{align} x^T (K_1 \otimes K_2) x & =x^T (AA^T \otimes BB^T) x \\ &= x^T (A\otimes B )(A^T \otimes B^T) x \\ &= y^T y \ge0,\text{ Let } y = (A^T \otimes B^T) x \end{align}\]性质3 若$K_1,K_2$为正定阵，则$K_1 \odot K_2$也为正定阵 , 也即Hadamard积保持正定性质，
第一种证明方法是观察到Hadamard积为Kronecker积的主子阵，因此该结论是显然的， \(K_1 \odot K_2 = P (K_1 \otimes K_2) P^T , \exists P\) 第二种证明方法是同样进行Cholesky分解，
\[\begin{align} x^T (K_1 \odot K_2) x &= x^T(AA^T \odot K_2) x \\ &= \sum_{ij} (\sum_k A_{ik} A_{jk} ) K_{ij} x_i x_j \\ &= \sum_k \sum_{ij}A_{ik}x_i A_{jk}x_j K_{ij} \\ &= \sum_k y_k^T K y_k \ge 0,\text{ Let } y_k = A_{k}^Tx \end{align}\]Kernel PCA
有了上面核方法的基础，下面可以进入核PCA（Kerne PCA）的世界中。
核PCA也即在$\varphi(X)$定义的空间内做PCA，考虑用矩阵$H= I - \frac{1}{n}ee^T$做归一化， $X$映射后的矩阵为$F$, 定义核矩阵 $K = FF^T$
核PCA需要计算 $H F$的奇异值分解（SVD，Singular Value Decomposition），也即计算$F^T H F$的特征值分解，而类似于PCoA，上述问题又可以转化为 $ HFF^T H = H KH$ 的特征值分解，这样就可以绕开了$F$而计算得到降维后的矩阵。
Probabilistic PCA
概率PCA从因子模型的角度看待PCA中的降维问题，概率PCA假定存在标准高斯因变量$Z$，和线性变换$W$ ,使得
\[X = WZ + \mu + \epsilon, Z \sim \mathcal{N}(0,I), \epsilon \sim \mathcal{N}(0,\sigma^2)\]不失一般性，可以仅考虑 $\mu=0$ 的情况，对$\mu \ne 0$ 的情况对$X$进行中心化处理即可，也即我们考虑，
\[X = WZ +\epsilon ,\epsilon \sim \mathcal{N}(0,\sigma^2)\]与线性回归不同之处在于，此时$Z$是未被观测到的隐变量，隐变量的存在也使得该问题比线性回归难求解得多。
Prediction
$(X,Z,\epsilon)$ 满足联合正态分布，其均值显然为$0$，可以计算其协方差矩阵，
\[\begin{align} \begin{pmatrix} X \\ Z \\ \epsilon \end{pmatrix} \sim \mathcal{N} ( \begin{pmatrix} 0 \\ 0 \\ 0 \end{pmatrix} , \begin{pmatrix} WW^T + \sigma^2 I_p & W & \sigma^2 \\ W^T & I_q & O \\ \sigma^2 & O & \sigma^2 \end{pmatrix} ) \end{align}\]利用正态分布的条件分布可以得到，关于该部分内容，可以移步至 高斯无向图模型
\[Z \vert X \sim \mathcal{N} (W^T(WW^T+\sigma^2 I_p)^{-1} X,I - W(WW^T +\sigma^2 I_p)^{-1} W^T)\]上式可以利用Shermann–Morrison–Woodbury公式化简，首先该公式为，
\[(A+UV^T)^{-1} = A^{-1} - A^{-1}U (I + V^T A^{-1}U)^{-1}V^T A^{-1}\]代入可以得到，
\[Z \vert X \sim \mathcal{N}(M, \sigma^2D^{-1}), \text{ Let } D= W^T W +\sigma^2I_q ,M =W^T(WW^T+\sigma^2 I_p)^{-1} X\]此时得到的隐变量$Z$就是一个低维的向量，利用其可以对$X$进行降维，可以取$Z$的条件期望进行预测，
\[E[Z \vert X] = W^T(WW^T +\sigma^2I_p)^{-1} X = (W^T W+\sigma^2 I_q)^{-1}W^TX\]利用上式，可以对条件分布进行进一步化简，
\[Z \vert X \sim \mathcal{N}(D^{-1}W^T X, \sigma^2D^{-1}), \text{ Let } D= W^T W +\sigma^2 I_q\]概率PCA本质上是一种生成模型，给出了高维数据$X$和低维数据$Z$之间的联系，因此其不仅可以用来做降维，在该生成模型上也有很多其他的用途，例如可以推广到变分自编码器上面。
上述说明了概率PCA模型的含义，以及如何使用该模型进行降维，但该模型的核心是未知参数的估计。
该问题的求解分为两种方法，第一种是基于极大似然估计，另一种基于EM算法，下面分别介绍。
MLE
由于我们知道,
\[X \sim \mathcal(\mu, C), \text{ Let } C = WW^T +\sigma^2 I_p\]给定了$X$的前提下，可以对其进行极大似然估计，对其对数似然函数，
\[L(X \vert \sigma^2,\mu,W) = \sum_{i=1}^N[-\frac{1}{2} \log \det C - \frac{1}{2} (X_i -\mu)^T C^{-1}(X_i-\mu) ]\]求导寻找其极值点，可以得到$\hat \mu$的极大似然估计为其样本均值，
\[\hat \mu = \bar X = \sum_{i=1}^N X_i\]将其代入上式，并且令$S = \frac{1}{n}(X - \bar X)^T (X -\bar X)$为样本方差矩阵，可以得到，
\[\begin{align} L &= -\frac{n}{2}\log \det C- \frac{1}{2} tr (X-\bar X)^T C^{-1} (X - \bar X) \\ &= -\frac{n}{2}\log \det C- \frac{n}{2} tr (C^{-1} S) \end{align}\]求解上式需要使用到矩阵微分的方法，参见了 知乎’矩阵求导与矩阵微分
用到了几个常用的矩阵微分的结论，
\(\begin{align} d \log \det C &= tr(C^{-1} dC) \\ d tr(C^{-1}S) &= tr(dC^{-1}S) \\ d C^{-1} &= -C^{-1} dC C^{-1} \end{align}\) 并且根据$C$的定义，可以得到，
\[\begin{align} dC &= dW W^T +W dW^T \\ dC &= d\sigma^2 I \end{align}\]Estimate $W$
首先来估计$W$，计算稍显繁琐，稍安勿躁，因为其结果非常简洁直观
通过令$\frac{dL}{dW}=0 $ 计算$W$的极大似然估计，
\(\begin{align} \frac{1}{2}tr(C^{-1} dC) + \frac{1}{2}tr(-C^{-1} dC C^{-1} S) &= \frac{1}{2}tr(C^{-1} dC) + \frac{1}{2}tr(-C^{-1}SC^{-1} dC )\\ &= \frac{1}{2}tr(C^{-1}dW W^T +C^{-1}W dW^T ) + \frac{1}{2}tr(-C^{-1}SC^{-1} (dW W^T +W dW^T) ) \\ &= tr(C^{-1} WdW^T) - tr(C^{-1}S C^{-1}W dW^T)\\ &= tr([C^{-1} W-C^{-1}S C^{-1}W ]dW^T)\\ \frac{dL}{dW^T} &= C^{-1} W-C^{-1}S C^{-1}W =0 \\ W&=SC^{-1}W \end{align}\) 下面利用特征值分解对上述方程进行求解，对$W^TW$作谱分解，对应的特征向量记作$V$,也即 $W^T W = V \Lambda V^T$，
\[\begin{align} W &= SC^{-1}W \\ W &= S(WW^T +\sigma^2 I_p)^{-1} W \\ W &= SW(W^TW +\sigma^2 I_q)^{-1} \\ W &=SW V(\Lambda +\sigma^2 I_q)^{-1}V^T \\ WV &= SWV(\Lambda +\sigma^2 I_q)^{-1} \\ (\Lambda +\sigma^2 I_q) WV &= SWV \\ (\Lambda +\sigma^2 I_q) WV \Lambda^{-\frac{1}{2}} &= SWV \Lambda^{-\frac{1}{2}}, \text{With } (WV)^T(WV) = \Lambda \end{align}\]可以观察到上式等价于对$S$做了一个特征值分解，最后一步的作用在于归一化特征向量为正交阵，
\[\begin{align} \text{Let }S &= \Phi \Gamma\Phi^T ,S \Phi = \Phi \Gamma \\ \text{Then } \Gamma &=(\Lambda +\sigma^2 I_q) ,\Phi = WV \Lambda^{-\frac{1}{2}} \end{align}\]可以解得$W$的估计，
\[\begin{align} W &= \Phi \Lambda^{\frac{1}{2}} V^T = \Phi (\Gamma - \sigma^2 I_q)^{\frac{1}{2}} V^T \end{align}\]上式中$V$为未知量，但由于$W$也是未知参数，而$V$仅代表特征向量的方向，而该方向并不会影响模型的效果，因此$V$实际上可以随意指定，例如为了方便可以直接指定$V=I$,注意上上式中$\sigma^2$仍未计算出来，下面进行计算。
Estimate $\sigma^2$
通过令$\frac{dL}{d \sigma^2}=0$获得其估计，
\[\begin{align} tr(C^{-1} dC) + tr(-C^{-1} dC C^{-1} S) &= tr(C^{-1} dC) + tr(-C^{-1}SC^{-1} dC )\\ &= tr(C^{-1}d\sigma^2) -tr(-C^{-1}SC^{-1} d\sigma^2) \\ &= tr(C^{-1}-C^{-1}SC^{-1})d\sigma^2 \\ \frac{dL}{d \sigma^2} &= tr(C^{-1}-C^{-1}SC^{-1}) = 0 \end{align}\]利用之前的Shermann–Morrison–Woodbury公式和在Prediction一节中得到的结论，
\[\begin{align} C^{-1} &= \frac{1}{\sigma^2} (I - W(WW^T +\sigma^2 I_p)^{-1} W^T) \\ SC^{-1} &= \frac{1}{\sigma^2} (S - SW(WW^T +\sigma^2 I_p)^{-1} W^T) \\ &= \frac{1}{\sigma^2} (S - S(W^TW +\sigma^2 I_p)^{-1} WW^T) \\ &= \frac{1}{\sigma^2}(S - SC^{-1}WW^T) \\ &= \frac{1}{\sigma^2}(S - WW^T) ,\text{ With } W= SC^{-1}W \\ C^{-1} - C^{-1}SC^{-1} &= C^{-1} (I-SC^{-1}) \\ &= \frac{1}{\sigma^2}C^{-1} (\sigma^2I-S+WW^T) \\ &=\frac{1}{\sigma^2}C^{-1}(C-S) \\ tr(C^{-1} - C^{-1}SC^{-1}) &=\frac{1}{\sigma^2} (tr(I_p) - tr(C^{-1}S)) \\ &=\frac{1}{\sigma^2} (tr(I_p) - tr(S- WW^T)) = 0 \\ \sigma^2 p &=tr(S) - tr(WW^T) \\ &=tr(S) - tr(\Lambda) \\ &=tr(S) - tr(\Gamma_q-\sigma^2 I_q) \\ \sigma^2 (p-q) &=tr(S) - tr(\Gamma_q) = \sum_{i=q+1}^N \lambda_i \\ \sigma^2 &= \frac{\sum_{i=q+1}^N \lambda_i}{p-q} \end{align}\]经过一番推导得到最终的结论，上式的含义是明显地，也即$\sigma^2$的估计为$W$的最后$p-q$个特征值的平均值，这也与最佳低秩逼近中的结论具有内在的一致性，可以移步至 特征值不等式与最佳低秩逼近
综上，我们总结极大似然估计的结果，
\[\begin{align} W &= \Phi \Lambda^{\frac{1}{2}} V^T = \Phi (\Gamma - \sigma^2 I_q)^{\frac{1}{2}} V^T \\ \sigma^2 &= \frac{\sum_{i=q+1}^N \lambda_i}{p-q} \end{align}\]取$V=I$的时候，上述的含义在于首先用最小的几个特征值估计噪声，然后将奇异值矩阵减去估计得到的噪声项得到投影矩阵，与PCA有异曲同工之妙。
EM
极大似然估计的结果精巧，但计算过程稍显复杂，而如果使用EM算法，可以很好地简化计算。
对于EM算法不熟悉的读者，可以移步至 EM算法
与极大似然估计中计算$X$的似然函数不同，EM算法中需要计算完全样本$(X,Z)$的似然函数，用对数似然函数表示为，
\[\begin{align} \max L(X,Z) &= \max [\sum_{i=1}^N \log P(X \vert Z) + \log P(Z)] \\ &= \min \sum_{i=1}^N [\frac{p}{2} \log \sigma^2 + \frac{(X_i-\mu- WZ_i)^T (X_i- \mu -WZ_i)}{2\sigma^2}+ \frac{1}{2} Z_i^TZ_i] \\ &=\min \frac{np}{2} \log \sigma^2+ \frac{1}{2\sigma^2} \sum_{i=1}^N (X_i-\mu)^T(X_i-\mu) - \frac{1}{\sigma^2} \sum_{i=1}^N(X_i-\mu)^T WZ_i +\frac{1}{2\sigma^2} \sum_{i=1}^N Z_i^T W^T WZ_i +\frac{1}{2}Z_i^TZ_i \\ &=\min \frac{np}{2} \log \sigma^2+ \frac{n}{2\sigma^2} tr (S) - \frac{1}{\sigma^2} \sum_{i=1}^N(X_i-\mu)^T WZ_i +\frac{1}{2\sigma^2} \sum_{i=1}^N tr(W^T WZ_i Z_i^T) +\frac{1}{2} tr(Z_iZ_i^T) \\ \end{align}\]上式的相反数在$p(Z \vert X)$的条件期望下即为EM算法中的Q函数，首先计算该条件期望，根据条件分布，
\[Z \vert X \sim \mathcal{N}(D^{-1}W^T X, \sigma^2D^{-1}), \text{ Let } D= W^T W +\sigma^2 I_q\]可以计算得到$Z$的条件期望和条件方差，由此得到其在该条件下的一阶矩和二阶矩，表示为，
\[\begin{align} \langle Z\rangle &= E_{Z \vert X}[Z] = D^{-1}W^T X \\ \langle ZZ^T \rangle &= E_{Z \vert X}[ZZ^T] = \sigma^2D^{-1}+ D^{-1}W^T X X^T W D^{-1} \end{align}\]在E-Step中计算出上式中的结果后，嵌入M-Step的Q函数中，
\[\begin{align} Q &= \frac{np}{2} \log \sigma^2+ \frac{n}{2\sigma^2} tr (S) - \frac{1}{\sigma^2} \sum_{i=1}^N tr(\langle Z_i \rangle (X_i-\mu)^T W) +\frac{1}{2\sigma^2} \sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) +\frac{1}{2} tr \langle Z_iZ_i^T \rangle \end{align}\]为了最大化Q函数，对其求导，
\[\begin{align} \frac{dQ}{dW} &= -\frac{1}{\sigma^2} \sum_{i=1}^N \langle Z_i \rangle (X_i - \mu)^T + \frac{1}{ \sigma^2} \langle Z_i Z_i^T \rangle W^T = 0 \\ \frac{dQ}{d\sigma^2} &=\frac{np}{2 \sigma^2} - \frac{n}{2 \sigma^4} tr(S) + \frac{1}{\sigma^4} \sum_{i=1}^N tr(\langle Z_i \rangle (X_i-\mu)^T W) -\frac{1}{2\sigma^4} \sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) =0 \\ \end{align}\]据此可以得到$W$的解，
\[\begin{align} W \sum_{i=1}^N \langle Z_iZ_i^T \rangle &= \sum_{i=1}^N (X_i -\mu) \langle Z_i^T \rangle \\ W&= \sum_{i=1}^N (X_i -\mu) \langle Z_i^T \rangle \sum_{i=1}^N\langle Z_iZ_i^T \rangle^{-1} \end{align}\]再代入关于$\sigma^2$的式子中求解，
\[\begin{align} np \sigma^2 &= ntr(S) - 2\sum_{i=1}^N tr(\langle Z_i \rangle (X_i-\mu)^T W) +\sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) \\ np \sigma^2 &=ntr(S) - 2\sum_{i=1}^N tr(\langle Z_iZ_i^T \rangle W^T W) +\sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) \\ np \sigma^2 &=ntr(S) - \sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) \\ \sigma^2 &= \frac{1}{p}tr(S) -\frac{1}{np}\sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) \end{align}\]总结EM算法，EM算法迭代地更新参数，
在M-Step中，最大化Q函数,得到参数估计，
\(\begin{align} \hat W&= \sum_{i=1}^N (X_i -\mu) \langle Z_i^T \rangle \sum_{i=1}^N\langle Z_iZ_i^T \rangle^{-1} \\ \hat \sigma^2 &= \frac{1}{p}tr(S) - \frac{1}{np}\sum_{i=1}^N tr(W^T W \langle Z_i Z_i^T \rangle) \end{align}\) 而在E-Step中，求隐变量的条件期望，
\[\begin{align} \langle Z\rangle &= D^{-1}W^T (X -\mu) \\ \langle ZZ^T \rangle &= \sigma^2D^{-1}+ D^{-1}W^T (X-\mu) (X-\mu)^T W D^{-1} \\ \text{With } D&= (W^TW+\sigma^2I_q) \end{align}\]可以看到在E-Step中，并不需要显式地求解条件概率 $P(Z \vert X)$， 而只需要求解上述的一阶矩和二阶矩，本质原因是概率PCA中的分布为高斯分布，而高斯分布的充分统计量正好就是其一阶矩和二阶矩。更一般地，对于指数族分布EM算法中均只需要关于其充分统计量求期望即可，关于这部分的理论推导，感兴趣的读者可以移步至 指数族分布
而且可以证明，EM算法将收敛到极大化似然函数的点，证明只需要代入两个方法的表达式进行验证即可，但相关计算同样较为繁琐，感兴趣的读者可以参见 张志华‘机器学习导论 中该部分对应内容的讲解。