正如我们已经看到的那样，概率PCA模型可以根据连续潜在空间$z$上的积分或求和来表示， 其中对于每个数据点$x_n$，都存在一个对应的潜在变量$z_n$。于是，我们可以使用EM算法来找到模型参数。这看起来似乎相当没有意义，因为我们已经得到了最大似然参数值的一个精确的解析解。然而，在高维空间中，使用迭代的EM算法而不是直接计算样本协方差矩阵可能会有一些计算上的优势。这个EM的求解步骤也可以推广到因子分析模型中，那里不存在解析解。最后，它使得我们可以用一种有理有据的方式处理缺失的数据。

我们可以使用一般的EM框架来推导用于概率PCA的EM算法。因此，我们写出完整数据对数似然函数，然后关于使用旧的参数值计算的潜在变量的后验概率分布求期望。最大化完整数据 对数似然函数的期望就会产生新的参数值。因为我们假定数据点是独立的，因此完整数据对数似然函数的形式为

$$\ln p(X,Z|\mu,W,\sigma^2) = \sum\limits_{n=1}^N\{\ln p(x_n|z_n) + \ln p(z_n)\} \tag{12.52}$$

其中矩阵$Z$的第$n$行由$z_n$给出。我们已经知道$\mu$的精确的最大似然解是式（12.1）定义的样本均值$\bar{x}$。在这个阶段将$\mu$替换掉是比较方便的。分别使用式（12.31）和（12.32）给出的潜在概率分布和条件概率分布的表达式，然后关于潜在变量上的后验概率分布求期望，我们有

$$\begin{eqnarray} \mathbb{E}[\ln p(X,Z|\mu,W,\sigma^2)] &=& -\sum\limits_{n=1}^N\Bigg\{\frac{D}{2}\ln(2\pi\sigma^2)+\frac{1}{2}Tr(\mathbb{E}[z_nz_n^T]) \\ & & +\frac{1}{2\sigma^2}\Vert x_n - \mu \Vert^2 - \frac{1}{2\sigma^2}\mathbb{E}[z_n]^TW^T(x_n - \mu) \\ & & +\frac{1}{2\sigma^2}Tr(\mathbb{E}[z_nz_n^T]W^TW) + \frac{M}{2}\ln(2\pi)\Bigg\} \tag{12.53} \end{eqnarray}$$

注意，上式仅仅通过高斯分布的充分统计量对后验概率分布产生依赖。因此在E步骤中，我们使用旧的参数计算

$\begin{eqnarray} \mathbb{E}[z_n] &=& M^{-1}W^T(x_n - \bar{x}) \tag{12.54} \\ \mathbb{E}[z_nz_n^T] &=& \sigma^2M^{-1} + \mathbb{E}[z_n]\mathbb{E}[z_n]^T \tag{12.55} \end{eqnarray}$

这可以直接从后验概率分布（12.42）以及标准的结果$\mathbb{E}[z_nz_n^T] = cov[z_n] + \mathbb{E}[z_n]\mathbb{E}[z_n]^T$中得出。这里，$M$由式（12.41）定义。

在M步骤中，我们关于$W$和$\sigma^2$进行最大化，保持后验统计量固定。关于$\sigma^2$的最大化很容易。对于关于$W$的最大化，我们可以使用（C.24）。求得的M步骤方程为

$$\begin{eqnarray} W_{new} &=& \left[\sum\limits_{n=1}^N(x_n - \bar{x})\mathbb{E}[z_n]^T\right]\left[\sum\limits_{n=1}^N\mathbb{E}[z_nz_n^T]\right]^{-1} \tag{12.56} \\ \sigma_{new}^2 &=& \frac{1}{ND}\sum\limits_{n=1}^N\{\Vert x_n - \bar{x} \Vert^2 - 2\mathbb{E}[z_n]^TW_{new}^T(x_n - \bar{x}) \\ & & + Tr(\mathbb{E}[z_nz_n^T]W_{new}^TW_{new})\} \tag{12.57} \end{eqnarray}$$

概率PCA的EM算法的执行过程为：对参数进行初始化，然后交替地在E步骤中使用式（12.54）和（12.55）计算潜在空间的后验概率分布的充分统计量，以及在M步骤中使用式（12.56）和（12.57）来修正参数的值。

用于PCA的EM算法的一个好处是对于大规模应用的计算效率（Roweis，1998）。与传统的基于样本协方差矩阵的特征向量分解的PCA不同，EM算法时迭代的，因此似乎没有什么吸引力。 然而，在高维空间中，EM算法的每次迭代所需的计算量都要比传统的PCA小得多。为了说明这一点，我们注意到，对协方差矩阵的特征分解的计算复杂度为$O(D^3)$。通常我们只对前$M$个特征向量和它们的特征值感兴趣，这种情况下我们可以使用$O(MD^2)$的算法。然而，计算协方差矩阵本身需要$O(ND^2)$的计算量，其中$N$是数据点的数量。有一些能够避免直接计算协方差矩阵的算法，例如快照方法（snapshot method）（Sirovich， 1987）假设特征向量是数据向量的线性组合，但是这种算法的计算复杂度为$O(N^3)$，因此不适用于大规模数据。这里描述的EM算法也没有显式地建立协方差矩阵。相反，计算量最大的步骤是涉及到对数据集求和的操作，计算代价为$O(NDM)$。对于较大的$D,M \ll D$，这与$O(ND^2)$相比，计算量极大地降低，因此可以抵消EM算法的迭代本质。

注意，这个EM算法可以用一种在线的形式执行，其中每个$D$维数据点被读入、处理，然后在处理下一个数据点之前丢弃这个数据点。为了说明这一点，注意在E步骤中需要计算的量（一个$M$维向量和一个$M \times M$的矩阵）可以分别对每个数据点单独计算，在M步骤中，我们需要在数据点上累积求和，这个可以增量地完成。如果$N$和$D$都很大，那么这种方法会很有优势。

由于我们现在对PCA建立了一个完全的概率模型，因此我们可以通过对未观测变量进行积分或求和的方式，处理缺失的数据，假设数据的缺失是随机的。与之前一样，这些缺失值可以使用EM算法进行处理。我们在图12.11中给出了使用这种方法进行数据可视化的一个例子。

EM算法的另一个特征是，我们可以取极限$\sigma^2 \to 0$，对应于标准的PCA，仍然可以得到一 个合法的类似EM的算法（Roweis， 1998）。根据式（12.55），我们看到我们在E步骤中需要计算的唯一的量是$\mathbb{E}[z_n]$。此外，M步骤可以得到简化，因为$M = W^TW$。为了强调算法的简化，让我们将$\tilde{X}$定义为一个$N \times D$的矩阵，它的第$n$行为向量$x_n − \bar{x}$，类似地，定义$\Omega$为一个$M \times N$的矩阵，它的第$n$行是向量$\mathbb{E}[z_n]$。这样PCA的EM算法的E步骤（12.54）就变成了

$$\Omega = (W_{old}^TW_{old})^{-1}W_{old}^T\tilde{X}^T \tag{12.58}$$

M步骤（12.56）的形式为

$$W_{new} = \tilde{X}^T\Omega^T(\Omega\Omega^T)^{-1} \tag{12.59}$$

与之前一样，可以使用一种在线的方式执行。这些方程有一个很简单的意义，如下所述。根据我们之前的讨论，我们看到E步骤涉及到数据点在当前估计的主子空间上的正交投影。对应地，M步骤表示对主子空间的重新估计，使得平方重建误差最小，其中投影固定。

我们可以给出这个EM算法的一个简单的物理类比，这对于$D = 2$和$M = 1$的情形很容易进行可视化。考虑二维空间中的一组数据点，令一维主子空间用一个固体的杆表示。现在使用一个遵守胡克定律(存储的能量正比于弹簧长度的平方)的弹簧将每个数据点与杆相连。在E步骤中，我们保持杆固定，让附着的点沿着杆上下滑动，使得能量最小。这使得每个数据点独立的到达对应的数据点在杆上的正交投影的位置。在M步骤中，我们令附着点固定，然后松开杆，让杆达到能量最小的位置。然后E步骤和M步骤不断重复，直到满足一个收敛准则，如图12.12所示。