PCA算法(Principal Component Analysis)揭秘

经典PCA算法

PCA算法的应用包括降维、有损数据压缩、特征抽取、数据可视化等。目前PCA算法有两个通用定义，能殊途同归，得到相同的结果。一方面，我们可以用正交投影来定义PCA，即将数据投影到更低维的线性子空间，也被称为主子空间(principal subspace)，需要满足最终投射数据的方差是最大化的。另一方面，一个等价的定义是寻找最小化平均投影代价的线性投影，即通过计算数据点和其投影的平均二次方距离来实现。正交投影的过程如图1所示：

(图1)

其中紫色的直线表示低维子空间，红色的数据点正交投影到该子空间，得到绿色的点，其中需要满足方差最大化原则。而PCA另一个可选的定义是，最小化二次投影误差之和，如蓝色的线段所示。

PCA的最大方差提法

考虑数据观测集 $\left \{ x_{n}\right \}$ ，其中 $n=1,...,N$ ，且 $x_{n}$ 是维度为D的欧式变量。我们的目标是投影这些数据到某个子空间，其中后者维度 $M< D$ ，且满足投影数据的方差最大化。此时我们假设M的值已知，后面我们会思考如何通过数据集来得到M的恰当值。

由浅入深，首先考虑投影到一维空间(M=1)，我们用D维向量 $u_{1}$ 表示该空间的方向，为方便起见(同时也不失一般性)，我们应该选择一个单位向量，即满足 $u_{1}^{T}u_{1}=1$ ，注意我们只对 $u_{1}$ 的方向感兴趣，而不是对它的绝对值(或强度)。此时，每个数据点 $x_{n}$ 可以被投影到一个标量值 $u_{1}^{T}x_{n}$ ，而被投影的数据平均值为 $u_{1}^{T}\bar{x}$ ，其中样本采样均值 $x\bar{}$ 由以下基本公式给出：

(式1)

而投影数据的方差可根据方差公式得到：

(式2)

其中上式等号右边的S是数据的方差矩阵，定义为：

(式3)

现在考虑如何最大化投影方差，即 $u_{1}^{T}Su_{1}$ ，其中把 $u_{1}$ 看成变量，即选择什么样的子空间 $u_{1}$ 能达到这个目标。很明显，这是一个带约束的最优化问题，因为 $||u_{1}||$ 不能趋于无穷大。另外，一个归一化约束条件是 $u_{1}^{T}u_{1}=1$ ，为了强制满足该条件，我们借助拉格朗日乘子法，其中系数用 $\lambda_{1}$ 表示，这样就得到一个无约束的最大化求解问题：

(式4)

通过设置上式对 $u_{1}$ 的导数为0，我们可得到一个稳定点满足：

(式5)

这意味着我们所求的 $u_{1}$ 必须是方差矩阵S的特征向量。若对上式等号两边左乘一个 $u_{1}^{T}$ ，并利用 $u_{1}^{T}u_{1}=1$ 的事实，我们将得到方差的公式：

(式6)

这表明当我们设置 $u_{1}$ 为具有最大特征值 $\lambda_{1}$ 的特征向量时，方差S就会达到最大值。这个特征向量也叫做第一个主成分。

我们可以增量的方式额外定义更多的主成分，只需要选择每个新的方向为最大化投影方差的方向，并满足与已有的方向正交即可。现在考虑M维投影空间的通用情况，最优线性投影(方差最大化)被定义为M个特征向量 $u_{1},u_{2},...,u_{M}$ ，而协方差矩阵S对应了M个最大的特征值 $\lambda _{1},\lambda _{2},...,\lambda _{M}$ 。

简单总结：PCA包含评估数据集的协方差矩阵S和均值 $x\bar{}$ ，并试图找到S矩阵的M个特征向量，对应M个最大的特征值。对于一个DxD矩阵，计算它所有特征向量分解的代价为 $O(D^{3})$ ，如果我们试图将数据投影到前M个主成分，我们只需要找到前M个特征值和特征向量，这能够用更高效的方式实现，比如power method，它的时间复杂度为 $O(MD^{2})$ ，或者也可以用EM算法。

PCA的最小误差提法

现在讨论另一种可选的PCA解释法，即基于投影误差的最小化。首先我们引入完备的D维正交基向量集合 $\left \{ u_{i} \right \}$ ，其中 $i=1,...,D$ ，它满足：

$u_{i}^{T}u_{j}=\delta _{ij}$ (式7)

由于这组基向量是完备的，因此任何数据点能够用它的线性组合精确表示：

(式8)

其中系数 $\alpha _{ni}$ 对于不同的数据点是不同的。这个相当于对原有坐标系统进行旋转，新系统用 $u_{i}$ 表示，并且原有的D个成分 $\left \{ x_{n1},..., x_{nD}\right \}$ 被等价的 $\left \{ \alpha _{n1},..., \alpha _{nD}\right \}$ 取代。根据基向量的定义和正交属性，每项系数 $\alpha _{nj}$ 是数据点与基向量 $u_{j}$ 的内积，即 $\alpha _{nj}=x_{n}^{T}u_{j}$ ，因此结合式8，可得：

(式9)

我们目标是利用 $M<D$ 个变量近似表示这个数据点，方式是通过投影到更低维的子空间。M维线性子空间不失一般性可以用前M个基向量表示，具体而言，每个数据点 $x_{n}$ 可近似表示为：

(式10)

其中 $\left \{ z_{ni} \right \}$ 这组系数依赖特定的数据点(即它是变化的)，而 $\left \{ b_{i} \right \}$ 则对所有数据点保持一致(即它是常量)。我们可自由选择不同的 $\left \{ u_{i} \right \}$ 、 $\left \{ z_{ni} \right \}$ 以及 $\left \{ b_{i} \right \}$ 来最小化降维带来的失真。作为失真的衡量，我们可用原始数据 $x_{n}$ 和近似数据 $x\tilde{}_{n}$ 的平均二次距离来表达，即我们要最小化如下目标：

(式11)

首先将 $\left \{ z_{ni} \right \}$ 看成变量，把式10代入式11，令式11对 $z_{ni}$ 的偏导为0，并利用正交条件，我们得到：

(式12)

其中 $j=1,...,M$ 。类似的，令式11对 $b_{i}$ 的偏导为0，同样利用正交属性，我们得到：

(式13)

其中 $j=M+1,...,D$ 。现在计算原数据和近似数据的差值，即 $x_{n}-x\tilde{}_{n}$ ，通过将式9减去式10，并结合式13，我们得到：

(式14)

从中我们看出 $x_{n}$ 和 $x\tilde{}_{n}$ 两个向量之间的错位存在于主子空间的正交中，因为式14等式右边是 $\left \{ u_{i} \right \}$ 的一个线性组合(其中i=M+1,...,D)。正如图1所示，投影点 $x\tilde{}_{n}$ 必须在主子空间中，虽然可以在该空间内自由移动，而最小误差则在正交投影时取得。

现在我们得到失真度量 $J$ 将 $\left \{ u_{i} \right \}$ 看成唯一变量的表达式：

(式15)

我们仍然需要对 $J$ 求最小化，当然它也是一个带约束的最小化问题，否则只要令 $u_{i}=0$ ，就能得到一个无意义的最小化结果。这个约束来自正交条件，而最终问题的解将用协方差矩阵的特征向量的扩展来表达。在正式求解之前，首先让我们获取一个直觉，考虑一个2维数据空间(D=2)，以及一个1维主子空间的场景(M=1)。我们需要寻找一个方向 $u_{2}$ 来最小化 $J=u_{2}^{T}Su_{2}$ ，当然要满足一个归一化约束条件 $u_{2}^{T}u_{2}=1$ 。利用拉格朗日乘子法，引入系数 $\lambda _{2}$ ，强制满足该约束，我们得到一个最小化目标函数 $J\tilde{}$ ：

(式16)

令式16对 $u_{2}$ 的偏导为0，我们得到 $Su_{2}=\lambda _{2}u_{2}$ ，这表明 $u_{2}$ 是 $S$ 的特征向量，而特征值为 $\lambda _{2}$ 。因此任意一个特征向量将对应一个失真度量的稳定点。为了找到 $J$ 的最小值，我们重新代入 $u_{2}$ 的解到式16，得到： $J=\lambda _{2}$ 。于是我们得到 $J$ 的最小值，只需要令 $u_{2}$ 为对应于最小的两个特征值的特征向量。反过来说，我们应该让主子空间与最大特征值对应的特征向量对齐。这个结果与我们的直觉一致，为了最小化平均投影二次距离，我们应该选择让主成分子空间经过数据点的均值，并与方差最大化的方向对齐。当所有特征值相等时，任何主方向的选择都会导致相等的 $J$ 结果。

对于任意的 $D$ 和 $M$ (满足 $M<D$ )，最小化 $J$ 的通解可通过选取 $\left \{ u_{i} \right \}$ 为协方差矩阵的特征向量得到，即有：

(式17)

其中 $i=1,...,D$ ，且特征向量 $\left \{ u_{i} \right \}$ 是正交向量集。此时对应的失真度量 $J$ 可简化为：

(式18)

表明 $J$ 只是与主子空间正交的所有特征向量对应的特征值的简单求和。而 $J$ 的最小值可通过选取最小的 $D-M$ 个特征值来得到。反过来说，主子空间中的特征向量对应了 $M$ 个最大的特征值。

到目前，尽管我们讨论的是 $M<D$ 的情况，PCA的结论对于 $M=D$ 的情况同样成立，后者相当于没有维度衰减，而只是对坐标轴进行了简单旋转，来对齐主成分。

与PCA技术类似的一个线性降维技术叫做CCA(即canonical correlation analysis)，二者区别是，PCA只处理单个随机变量，而CCA能够考虑两个甚至更多个变量，本文不深入讨论。

PCA的应用

数据压缩

PCA在数据压缩方面的应用可用一个数字的案例(图2用数字3)说明，由于协方差矩阵的每个特征向量来自原有D维向量空间，我们可将特征向量表示为与原数据相同尺寸的图像。如图2所示，图像数字3的前5个特征向量以及对应的特征值分别用5张图像表示，从左到右依次是平均向量 $x\bar{}$ 和PCA的4个特征向量 $u_{1},...,u_{4}$ ，其中对应的特征值分别用 $\lambda _{1},...,\lambda _{4}$ 表示。

(图2)

该例子完整的特征值光谱(降序表示)如图3a所示。

(图3)

与图2对应，i越大， $\lambda _{i}$ 越小。与此同时，根据式18，失真度量 $J$ 是 $M+1$ 到 $D$ 的特征值之和， $J$ 与 $M$ 的关系也是 $M$ 越大， $J$ 越小，其具体的函数曲线如图3b所示。

如果将式12和式13代入式10，我们得到PCA中数据向量的近似表示：

其中式20中 $x\bar{}$ 定义为：

(式21)

很容易验证，只需把式21代入式20，由于 $x\bar{}$ 是 $(x\bar{}^{T}u_{i})u_{i}$ 从i=1到D的求和，它与该项的从i=1到M的求和的差即为i=M+1到D的求和，而另外一项保持不变，因此可得式19。

式19和式20体现了对原始数据的压缩，因为对每个点，我们用M维向量(即包含成分 $x_{n}^{T}u_{i}-x\bar{}^{T}u_{i}$ )来代替D维向量 $x_{n}$ ，而M越小，压缩程度越大。数字样例的PCA重构如图4所示。

(图4)

M是主成分的个数，图4表明：M越大，图像越清晰，即重构越准确；M越小，图像越模糊，即压缩程度越高。而当M=D时是最完美的，此时不存在压缩。

数据预处理

PCA另一个应用是数据预处理。此时不是为了降维，而是转换数据集来标准化它的某些属性，这在将来的模式识别算法中很重要。一般来说，当原始变量具有不同的单位，或有显著的变异性时，就需要该操作。比如我们在对一些数据集执行K-means算法之前，会先对变量进行单独的线性重缩放，来保证每个变量具有零均值和单位方差。标准化之后的协方差矩阵的每个元素可表示为：

(式22)

其中 $\sigma _{i}$ 是 $x_{i}$ 的方差。该矩阵也被称为相关性矩阵，它有一个属性，即当 $x_{i}$ 和 $x_{j}$ 是完美相关的时候， $\rho _{ij}=1$ ，反之若二者完全无关， $\rho _{ij}=0$ 。而当我们使用PCA时，我们将得到更优的归一化结果，即零均值和单位方差，这将对所有不同变量解除相关性。为了实现该目标，我们首先将特征向量公式(式17)写成如下形式：

$SU=UL$ (式23)

其中 $L$ 是DxD的对角阵，其每个对角元素是 $\lambda _{i}$ 。 $U$ 是DxD的正交矩阵，每个列是 $u_{i}$ 。我们对每个数据点 $x_{n}$ 执行如下转换：

(式24)

其中 $x\bar{}$ 是式1定义的样本均值，而 $L^{-1/2}$ 是 $L^{1/2}$ 的逆矩阵。显然我们得到的 $\left \{ y_{n} \right \}$ 具有零均值，而它的协方差是单位矩阵，因为有：

(式25)

上式第一个等号直接将式24代入定义。第二个等号利用了方差矩阵 $S$ 的定义(即式3)。第三个等号利用式23，将等式两边左乘 $U^{T}$ ，得到 $U^{T}SU=U^{T}UL$ ，因为 $U$ 是正交矩阵，因此 $U^{T}SU=L$ 。思考：第四个等号如何得到，你知道吗？

这项技术在预处理中叫做数据的白化或球形化，对著名的Old Faithful数据集的预处理效果如图5所示。其中左侧是原始数据，中间是单变量标准归一化后的结果，红色的垂直交叉线是主轴，其范围是 $\pm \lambda _{i}^{1/2}$ ，右侧是白化后的结果，它具有零均值和单位协方差。

(图5)

数据可视化

数据可视化也是PCA的应用之一。在图6的oil ﬂow数据集案例中，每个数据点被投影到2维( $M=2$ )的主子空间中，此时 $x_{n}$ 在笛卡尔坐标中用 $x_{n}^{T}u_{1}$ 和 $x_{n}^{T}u_{2}$ 表示，其中 $u_{1}$ 和 $u_{2}$ 是分别对应最大和次大特征值的特征向量。

(图6)

其中绿色、红色、蓝色分别对应3种不同的流配置数据。

PCA在高维数据的应用

在有些应用中，数据点的个数少于数据空间的维度，比如当只有几百个图像时，每个图像可能对应几百万维度的空间(这里指三色的像素值)，即 $N<D$ ，而它的线性子空间的维度最大只有 $N-1$ ，因此在应用PCA时只有极少数点的值大于 $N-1$ 。事实上在执行PCA时，我们会发现至少有 $D-N+1$ 个特征值为0，意味着在对应的特征向量方向上数据集具有零方差。一般来说，计算DxD矩阵的特征向量算法的时间复杂度为 $O(D^{3})$ ，因此直接应用PCA会导致计算不可行。

我们可用如下方式解决，首先定义 $X$ 为NxD的中心标准化矩阵，即它的每一行为 $(x_{n}-x\bar{})^{T}$ ，注意此时式3定义的协方差矩阵可写为 $S=N^{-1}X^{T}X$ ，而对应的特征向量等式(利用式17)变为：

(式26)

上式两边同时左乘一个 $X$ ，得到：

(式27)

假设把上式的 $Xu_{i}$ 记为 $v_{i}$ ，即 $Xu_{i}=v_{i}$ ，我们得到：

(式28)

这是NxN方阵 $N^{-1}XX^{T}$ 的特征向量等式。我们发现该方阵具有和原协方差矩阵一样的 $N-1$ 个特征值，当然前者具有额外的 $D-N+1$ 个零特征值。因此我们可以在 $O(N^{3})$ 的复杂度内计算特征向量，而不再需要 $O(D^{3})$ 。为了确定特征向量，我们对式28两边分别左乘一个 $X^{T}$ ，得到：

(式29)

从中可看出， $X^{T}v_{i}$ 就是矩阵 $S$ 的特征向量，其对应的特征值为 $\lambda _{i}$ 。注意：该特征向量尚未归一化。为了近似归一化，我们对 $X^{T}v_{i}$ (即 $u_{i}$ )进行缩放，以使 $||u_{i}||=1$ ，只需对 $u_{i}$ 执行如下除法转换：

(式30)

简单来说，在高维数据中应用PCA，我们首先计算 $XX^{T}$ ，并得到对应的特征值和特征向量，然后利用式30计算在原数据空间中的特征向量。

概率化PCA(Probabilistic PCA)

经典PCA将数据投影到更低维的线性子空间，事实上，PCA也可以用概率隐变量模型的极大似然估计方法来解释。这种PCA提法，也叫做概率化PCA(Probabilistic PCA，或简称PPCA)。它比传统PCA有如下优势：

PPCA表示高斯分布的约束形式，其中尽管自由参数的个数被约束，但允许模型捕获数据集的主要相关性。
我们可为PPCA推导出EM算法，其计算高效，只有少量主特征向量需要被计算，从而避免了计算整个协方差矩阵。
概率模型和EM算法的结合使得我们可以处理丢失数据。
多个PPCA模型可被有序混合，并用EM算法来训练。
PPCA是贝叶斯PCA的基础，在后者中，主子空间的维度可被自动发现。
PPCA中似然函数的存在允许它和其他概率密度模型的直接比较，与之不同，传统PCA会给那些接近主子空间的数据点赋予较低的重构代价，即使它们和原始训练数据有较大的距离。
PPCA可用于建模基于类别的条件密度，因此可被用于分类问题。
PPCA可被用于生成式的提供来自某个分布的样本数据。

PPCA是线性高斯框架下的一个简单模型，其中所有的边际分布和条件分布都是高斯分布。为了解释PPCA，首先引入对应于主子空间的显式隐变量z。然后再定义先验高斯分布p(z)，以及高斯条件分布p(x|z)，其中x是基于隐变量z的观测变量。需要指出，p(z)是一个零均值、单位协方差的高斯分布，即：

(式31)

而条件分布p(x|z)的表达式为：

(式32)

其中x的均值是z的一般线性函数，其中系数和截距分别是DxM维的W矩阵和D维的μ向量。注意：这也是一个朴素贝叶斯模型。后面我们会发现，W矩阵的各列跨越了数据空间内对应于主子空间的线性子空间。另外需要强调，假设零均值、单位协方差并不失一般性，因为即使一般的高斯模型也会得到一样的结果。

我们可以通用视角看待PPCA，首先得到隐变量，然后采样得到观测变量，比如D维观测变量x可看成M维隐变量z的线性转换，即：

$x=Wz+\mu +\varepsilon$ (式33)

其中z是M维高斯隐变量， $\varepsilon$ 是D维零均值高斯噪音，其协方差为 $\delta ^{2}I$ 。PPCA的工作过程如图7所示。左图是隐变量z的概率密度，中图表明观测值x可按如下方式得到：首先根据p(z)对隐变量z做出一个 $z\hat{}$ ，然后利用各向同性高斯分布作出x的值(如红色圈圈所示)，其均值为 $wz\hat{}+\mu$ ，协方差为 $\delta ^{2}I$ 。右图绿色椭圆为p(x)边际密度函数的轮廓线。整个过程是从隐空间到数据空间的映射，后面讲到的贝叶斯PCA将会是相反的过程。

(图7)

为了利用MLE确定 $W$ 、 $\mu$ 、 $\delta ^{2}$ 的值，我们需要写下p(x)的表达式，根据 $p(x)=\int p(x,z)dz$ ，以及 $p(x,z)=p(x|z)p(z)$ ，可得：

$p(x)=\int p(x|z)p(z)dz$ (式34)

作为线性高斯模型，p(x)也满足高斯分布，即有：

(式35)

其中 $C$ 是DxD协方差矩阵，定义为：

(式36)

上式也可通过对式33求均值和协方差直接得到，具体是：

其中用到一个事实是：z和 $\varepsilon$ 是独立的随机变量，因此是不相关的。直觉上，我们可理解p(x)为拿到一个各向同性高斯“喷雾罐”，并在主子空间移动喷涂高斯墨水(其密度由 $\delta ^{2}$ 决定)，以先验分布确定权重。累积的墨水密度最终得到一个饼状的边际分布p(x)。

虽然p(x)由 $W$ 、 $\mu$ 、 $\delta ^{2}$ 三组参数决定，但是这里面存在冗余，因为存在隐空间坐标轴的旋转。为了说明，首先考虑一个矩阵 $W\tilde{}=WR$ ，其中 $R$ 是正交矩阵。利用正交属性 $RR^{T}=I$ ，我们发现 $W\tilde{}W\tilde{}^{T}$ 就是 $WW^{T}$ ，意味着 $W\tilde{}W\tilde{}^{T}$ 存在于协方差矩阵 $C$ 中。具体计算如下：

(式39)

因此该量独立于 $R$ 。可以说这里存在一整个 $W\tilde{}$ 矩阵族群，每个都能得到一样的分布。这种不变性可从隐空间旋转的角度来理解，后面会回到对独立参数个数的讨论中来。

当评估预测分布时，我们需要计算 $C^{-1}$ ，即DxD矩阵的逆。利用矩阵逆的求解，可得 $C^{-1}$ ：

(式40)

其中 $M$ 是MxM的矩阵，具体为：

(式41)

由于我们直接对 $M$ 求逆，而不是直接对 $C$ 求逆，计算 $C^{-1}$ 的代价从 $O(D^{3})$ 缩减到 $O(M^{3})$ 。

除了计算p(x)，我们还需计算后验分布p(z|x)，根据线性高斯模型，可直接得：

(式42)

从上式可看出：该后验分布的均值依赖 $x$ ，而后验协方差独立于 $x$ 。

最大似然PCA

现在考虑如何用最大似然估计(MLE)来决定模型参数，对于观测数据集 $X=\left \{ x_{n} \right \}$ ，已知式35，其对应的对数似然函数可直接得到，为：

(式43)

令上式对 $\mu$ 的偏导为0，得到我们期望的结果 $\mu =x\bar{}$ ，其中 $x\bar{}$ 是式1定义的样本均值。回代之后，我们得到上式的新形式：

(式44)

其中 $S$ 是式3定义的协方差矩阵。由于最大似然函数是 $\mu$ 的二次函数，其解具有唯一的最大值。令式43对 $W$ 、 $\delta ^{2}$ 求偏导并求最大值则更复杂，但同样有闭合的精确解。可以证明，所有对数似然函数的稳定点可写为：

(式45)

其中 $U_{M}$ 是DxM矩阵，其各列是协方差矩阵 $S$ 的特征向量集合的大小为M的任意子集。 $L_{M}$ 是MxM对角矩阵，其元素来自特征值 $\lambda _{i}$ 。 $R$ 是任意MxM正交矩阵。

更进一步，有学者指出似然函数的最大值在当M个特征向量对应的特征值是最大的前M个解时达到(其余解为鞍点)。假设特征向量按照特征值降序排列，那么M个主特征向量是 $u_{1},...,u_{M}$ 。此时 $W$ 的列定义了标准PCA的主子空间。 $\delta ^{2}$ 对应的最大似然解为：

(式46)

也就是说 $\delta _{ML}^{2}$ 是被丢弃维度的平均方差。由于 $R$ 是正交的，它可理解为MxM隐空间的旋转矩阵。当把 $W$ 的解代入 $C$ 的表达式，并使用正交属性 $RR^{T}=I$ ，我们发现 $C$ 与 $R$ 独立。这意味着隐空间在旋转过程中的预测密度不变。特别当 $R=I$ 时， $W$ 的列是受方差参数 $\lambda _{i}-\delta ^{2}$ 控制缩放的主特征向量。一旦我们认识到对于独立高斯分布（在本例中为潜在空间分布和噪声模型）的卷积，方差是相加的，这些缩放因子的解释就很清楚了。因此在 $u_{i}$ 方向的方差 $\lambda _{i}$ 由 $\lambda _{i}-\delta ^{2}$ 部分从单位方差隐空间分布到数据空间的投影(通过相应的 $W$ 矩阵的列)加上方差 $\delta ^{2}$ 的各项同性部分(它通过噪音模型被加到所有的方向上)之和构成。

考虑协方差矩阵(式36)的形式，由单位向量 $v$ 指定的方向上预测分布的方差为 $v^{T}Cv$ ，其中 $v^{T}v=1$ 。假设 $v$ 与主子空间正交，即它由一些特征向量的线性组合构成。因此有 $v^{T}U=0$ ，以及 $v^{T}Cv=\delta ^{2}$ 。模型预测的噪音方差正交于主子空间，如式46所示，它是被丢弃特征值的平均值。假设 $v=u_{i}$ ，以及 $u_{i}$ 是主子空间其中一个被保留的特征向量。因此有 $v^{T}Cv=(\lambda _{i}-\delta ^{2})+\delta ^{2}=\lambda _{i}$ ，这意味着该模型正确捕获了主轴上的数据方差，并且其余方向上的方差近似用单个均值 $\delta ^{2}$ 表示。

一种构造最大似然密度模型的方式是找到协方差矩阵的特征值和特征向量，并基于以上结果评估 $W$ 、 $\delta ^{2}$ 。为了方便，可以选择 $R=I$ ，但是如果最大似然解是用数值优化的方式找到的话，比如用共轭梯度法，或EM算法，那么 $R$ 的结果值本质上是任意的。这表明 $W$ 的各列不需要正交。如果需要正交基向量，那么矩阵 $W$ 也可以被恰当的后处理。另外，EM算法也可被修改来得到正交主方向，并以对应特征值的降序排序。

隐空间的旋转不变性代表一种统计上的非不可识别性，类似于混合模型中的离散隐变量。我们考虑 $M=D$ 的场景，此时没有降维，因此有 $U_{M}=U$ ，以及 $L_{M}=L$ 。利用正交属性 $UU^{T}=I$ 和 $RR^{T}=I$ ，那么协方差矩阵 $C$ 就变成：

(式47)

上式利用了 $C$ 的定义(式36)以及 $W$ 的定义(式45)， $S$ 的属性(式23)以及正交属性。 $C$ 就是 $S$ ，这意味着什么？说明在标准最大似然解中，协方差矩阵就来自样本协方差。

与传统PCA的从D维数据空间投影到M维子空间不同，PPCA从隐空间投影到数据空间(利用式33)。对于可视化和数据压缩应用，可用贝叶斯理论来逆转该映射。

由式42可得z对x的条件期望为：

(式48)

其中 $M$ 由式41给出。这将投影到数据空间的以下点：

(式49)

该式和标准线性回归一致，也是线性高斯模型的最大似然函数的结果。类似的，从式42得到后验协方差为 $\delta ^{2}M^{-1}$ ，也是独立于x。

令 $\delta ^{2}\rightarrow 0$ 时，后验均值趋于如下表达式：

(式50)

它表示数据点到隐空间的正交投影，此时恢复到标准PCA模型。最终我们发现PPCA在定义多变量高斯分布时的一个重要角色，就是允许模型捕获主相关性的同时，控制模型自由度，即独立参数的个数。一般高斯分布的协方差矩阵有D(D+1)/2个独立参数，当然均值也有D个额外参数，因此参数个数与D成二次方关系，在高维场景代价较大。如果限制协方差矩阵为对角阵，它就只有D个独立参数，此时就是线性关系。此时它把变量看成独立的，因此无法表达相关性。PPCA提供了一种优雅的折中，一方面能捕获M个最显著的相关性，另一方面确保参数总数与D保持线性关系。矩阵C依赖参数W(维度DxM)，而 $\delta ^{2}$ 包含参数DM+1个。但要指出，这里的参数存在冗余，由于坐标轴的旋转。正交矩阵R体现了旋转的尺寸为MxM。该矩阵的第一列一共有M-1个独立参数，这是由于列向量必须归一化到单位长度。第二列有M-2个独立参数，除了归一化，还必须与前一列正交。以此类推。总结一下，R一共有M(M-1)/2个独立参数。因此协方差矩阵C的自由度为：

$DM+1-M(M-1)/2$ (式51)

显然该模型当M固定时，独立参数个数与D成线性相关。当取M=D-1时，恢复到全协方差高斯的标准结果。此时D-1个线性独立方向的方差被W的列控制，其余方向的方差为 $\delta ^{2}$ 。当M=0时，模型为各向同性协方差的情况。

PCA的EM算法

如前文所述，PPCA可用连续隐空间z的边际分布来表达，其中每个数据点 $x_{n}$ 对应了一个 $z_{n}$ 。此时可用EM算法来寻找模型参数的极大似然估计值。在高维空间中，使用迭代式的EM算法比传统的直接计算协方差矩阵具有一定优势。EM算法还可以进一步扩展到FA模型(factor analysis model)，但是后者没有闭合解。

PPCA的EM算法可用EM通用框架来推导。首先写出完全数据对数似然函数，并求它的期望。对该期望求最大值，得到新的参数值。因为数据点被假设独立，完全数据对数似然函数为如下形式：

(式52)

其中Z矩阵的第n行为 $z_{n}$ 。我们已知 $\mu$ 精确的最大似然解为样本均值 $x\bar{}$ (由式1给出)，在该阶段对 $\mu$ 进行替换也比较方便。利用式31和式32，引入高斯概率分布，我们得到式52的期望，为：

(式53)

在E步中，我们用旧的参数值计算出式53所需的以下期望值：

上面两式基于式42定义的后验分布，以及标准结果 $E[z_{n}z_{n}^{T}]=cov[z_{n}]+E[z_{n}]E[z_{n}]^{T}$ ，其中 $M$ 由式41定义。

在M步中，令式53对 $W$ 和 $\delta ^{2}$ 求最大化，其中对 $\delta ^{2}$ 的最大化是直接的，而对 $W$ 的最大化将用到矩阵的迹的一个求偏导性质，即：

(待证明)

为此我们得到M步的等式(即求 $W$ 和 $\delta ^{2}$ 新值的迭代式)：

PPCA的EM算法首先初始化参数，然后在E步中迭代计算潜在空间后验分布的统计值(利用式54、式55)，并在M步修改参数值(利用式56、式57)。EM算法的优势是计算高效，不像传统PCA基于协方差矩阵的特征向量分解，EM算法是迭代式的，每一轮都比传统的高效，尤其是高维空间。协方差矩阵的自分解需要计算量 $O(D^{3})$ ，而一般我们只对前M个特征向量和特征值感兴趣，因此只需要 $O(MD^{2})$ 的计算量。但是协方差矩阵本身就需要 $O(ND^{2})$ 的计算量，其中N是数据点个数。EM算法没有显式构造协方差矩阵，其中计算量最大的步骤是求和过程，为 $O(NDM)$ ，对于大的 $D$ ，即 $M<<D$ ，这个比起 $O(ND^{2})$ 节省了不少计算量。

EM算法还可用在线(on-line)方式实现，一边读入D维数据，一边处理，然后丢弃。注意到E步每个数据点的计算可被单独执行，M步中累计求和也可以增量的方式进行。该方法在D和N很大时具有明显优势。

EM算法另一个优雅的特征是当 $\delta ^{2}\rightarrow 0$ ，仍能使用EM算法。由式55可知，E步唯一需要计算的量是 $E[z_{n}]$ ，而且M步也被简化了，因为 $M=W^{T}W$ 。为了说明算法的简单性，我们定义 $X\tilde{}$ 为NxD的矩阵，其中第n行是向量 $x_{n}-x\bar{}$ ，同时定义 $\Omega$ 为DxM的矩阵，其第n行为向量 $E[z_{n}]$ 。那么E步的公式(式54)变为：

(式58)

M步公式(式56)也变为：

(式59)

需要指出这些都能用在线方式实现。并且这些公式理解也很简单，在E步包含数据点到主子空间的正交投影，M步包含对主子空间的重新估计，为的是最小化二次重构误差。我们可以给出EM的物理类比，特别当D=2以及M=1时，就能简单可视化。在二维数据点集合中，让一维主子空间用实线表示，让数据点附着到线段上，在E步中保持线段固定，允许数据点上下滑动。这会导致数据点在正交投影方向有独立的位置。在M步中保持点固定，让线段滑动到最优的位置。E步和M步不断重复直到收敛条件的满足，如图8所示。

(图8)

贝叶斯PCA

到目前的讨论中，都假定PCA主子空间的维度M已知，现实中我们必须选择合适的值，为了可视化，一般选取M=2，但在其他应用中，合适的M值并不清楚。一个方式是画出特征值光谱(如图3a所示)，并查看是否天然形成了两个具有巨大差异的类别，其中一组值更小，另一组值更大，若成立，则意味着M的天然取值。现实中，这样的差异并不常见。

由于PPCA有一个定义良好的似然函数，我们可使用交叉验证来决定维度的值，即通过在验证集中选择最大的对数似然。但这种方法计算代价较高，尤其当涉及PCA的概率混合模型时，此时我们需要为每个成分确定恰当的维度。

既然我们有了PPCA，自然而然我们想到用贝叶斯的方式去建模。为了做这个事，我们需要分别对模型参数 $\mu$ 、 $W$ 、 $\delta ^{2}$ 进行边际化。这可以通过使用变分框架来近似分析上难以处理的边际化来完成。变分下界得到的边际似然值将与不同的M值以及给出所选最大边际似然的值进行比较。

这里介绍一个更简单的方法，基于证据近似(evidence approximation)，当数据点个数较多，对应的后验分布的峰值较高的时候，该方法比较合适。它包含对W指定的先验选择，并允许主子空间中过剩维度被修剪掉。这对应了一种技术叫自动相关性确定(automatic relevance determination或ARD)。具体来说，我们基于W每列定义独立高斯先验分布，它表示主子空间的向量。每个高斯有个独立的方差，它由超参数 $\alpha _{i}$ 精确控制，具体如下：

(式60)

其中 $w_{i}$ 是 $W$ 矩阵的第i列。 $\alpha _{i}$ 可通过最大化边际概率函数迭代得到。该优化的结果是，一些 $\alpha _{i}$ 可能趋于无穷大，而对应的参数向量 $w _{i}$ 趋于零，得到一个稀疏解。主子空间的有效维度由有限 $\alpha _{i}$ 的个数决定，相应的向量 $w _{i}$ 可被认为与数据分布建模“相关”。这种贝叶斯方法自动在"提升与数据的拟合"和"降低模型复杂度"之间做了权衡利弊，前者通过使用更多的 $w _{i}$ 向量和对应的特征值 $\lambda _{i}$ ，后者通过抑制部分 $w _{i}$ 向量。

$\alpha _{i}$ 在训练中通过最大化对数边际概率被重新计算，该边际似然函数为：

(式61)

其中 $p(X|W,\mu ,\delta ^{2})$ 由式43定义。为了简化，也将 $\mu$ 、 $\delta ^{2}$ 看成待估计的参数，而不是把它们当成先验参数。该积分不好处理，我们可使用拉普拉斯近似。我们假设该后验分布的峰值很陡峭，那么 $\alpha _{i}$ 的估值就有如下形式：

(式62)

注意其中 $w _{i}$ 向量的维度是D。这些估值被EM算法迭代更新，来决定 $W$ 和 $\delta ^{2}$ 。E步再次由式54、式55确定，M步也再次由式57确定。唯一的变化是M步的W，它变为：

(式63)

其中 $A$ 是 $\alpha _{i}$ 构成的对角阵，即 $A=diag(\alpha _{i})$ 。而 $\mu$ 和前文一样，是样本均值。如果我们取 $M=D-1$ ，并且所有的 $\alpha _{i}$ 是有限值，该模型就变为全协方差高斯；而如果所有的 $\alpha _{i}$ 变为无限大，那么模型就变为各向同性高斯，此时该模型将包含主子空间有效维度的所有允许值。也可以考虑更小的M值，这虽然可以节省计算，但也会限制子空间的最大维度。该算法和标准PPCA的比较结果如图9所示。

(图9)

图中 $W$ 矩阵每个元素用一个方形表示，其中白色是正值，黑色是负值，其面积正比于元素的大小。合成数据集包括D=10维的300个数据点，这些数据点是从高斯分布中采样的，对于这10维中的数据集，在3个方向上标准差为 1.0，其余7个方向上的标准差为 0.5，也就是说M=3个方向的方差大于剩下7个方向。左图是最大似然PPCA的结果，右图是贝叶斯PCA的结果。我们看到贝叶斯模型能够通过抑制6个剩余自由度从而发现恰当的维度。