PCA数学原理再探讨

PCA（Principal Component Analysis） 是一种常见的数据分析方式，常用于高维数据的降维，可用于提取数据的主要特征分量。我之前也写过一篇简单的 PCA 的数学理解，日前看到一篇文章讲解的很棒，就整理转载过来，作为以后复习的总结。

理解 PCA，我觉得从其数学推导入手是非常好的途径。PCA 的数学推导可以从最大可分型和最近重构性两方面进行，前者的优化条件为划分后方差最大，后者的优化条件为点到划分平面距离最小，这里我将从最大可分性的角度进行证明。

特征表示和降维

在机器学习任务中，一般用一个特征向量描述一个样本，每个样本都有很多维特征，在机器学习中，一条样本由一个列向量表示。但是对于样本而言，某些特征之间有线性关系甚至直接关系。

举个例子，假如某学籍数据有两列 M 和 F，其中 M 列的取值是如何此学生为男性取值 1，为女性取值 0；而 F 列是学生为女性取值 1，男性取值 0。此时如果我们统计全部学籍数据，会发现对于任何一条记录来说，当 M 为 1 时 F 必定为 0，反之当 M 为 0 时 F 必定为 1。在这种情况下，我们将 M 或 F 去掉实际上没有任何信息的损失，因为只要保留一列就可以完全还原另一列。

当然上面是一个极端的情况，在现实中也许不会出现，不过类似的情况还是很常见的。例如上面淘宝店铺的数据，从经验我们可以知道，“浏览量”和“访客数”往往具有较强的相关关系，而“下单数”和“成交数”也具有较强的相关关系。这里我们非正式的使用“相关关系”这个词，可以直观理解为“当某一天这个店铺的浏览量较高（或较低）时，我们应该很大程度上认为这天的访客数也较高（或较低）”。这种情况表明，如果我们删除浏览量或访客数其中一个指标，我们应该期待并不会丢失太多信息。因此我们可以删除一个，以降低机器学习算法的复杂度。

上面这些让我们能够隐约感觉到特征之间的相关性，这些特征就可以去重。另一方面，对于样本个数比较少但是特征非常丰富的训练样本，少量样本在高维特征空间中就变成了稀疏分布，模型拟合就比较发散，而且计算量增加，通过降维把高维空间的特征映射到低维空间中，也就是降低特征的维度，可以减少计算量，加快模型拟合。以上就是降维的朴素思想描述和动机。

向量内积

文章首先介绍了关于向量内积的知识，这里不再对这些基础进行总结了。只需要记住：设向量 B 的模为 1，向量 A 与 B 的内积值等于 A 向 B 所在直线投影的矢量长度。这就是内积的一种几何解释，也是我们得到的第一个重要结论。

基

向量的基（也称基向量）用于描述一个多维空间。在多维空间中，如果我们要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。关于多维空间的基，这里也不啰嗦了，线性代数基础知识，但是需要记住一个重要结论，也是我一直比较推崇的，理解关于矩阵相乘的物理解释：两个矩阵相乘的意义是将右边矩阵中的 每一列 向量 $a_i$ 变换到左边矩阵中以 每一行 向量为基所表示的空间中去。

这里也解释了为什么一个样本数据要用一个列向量来表示，因为这个列向量通过矩阵相乘的形式被变换到其他特征空间（特征维度可能不一样）中，就可以得到这个样本在另一个特征空间的描述了。

最大可分性

描述一个多维空间中的某一向量，选择不同的基可以对同样一组数据给出不同的表示，如果基的数量少于向量本身的维数，则可以达到降维的效果。

这里说的向量本身的维度其实也是数据的一种描述方式，比如如果一条数据的向量维度是13，则表示这条数据用13个基来描述，第 i 个特征值是该数据在第 i 个基上的投影长度。我们用 PCA 降维就是找到另一种用更少的基描述数据的方式。

但是我们还没回答一个最关键的问题：如何选择基才是最优的。或者说，如果我们有一组 N 维向量，现在要将其降到 K 维（K 小于 N），那么我们应该如何选择 K 个基才能最大程度保留原有的信息？一种直观的看法是：希望投影后的投影值尽可能分散，因为如果重叠就会有样本消失。当然这个也可以从熵的角度进行理解，熵越大所含信息越多，这就是最大可分性的出发点。

方差

方差用于描述一组数组的分散程度，变量的方差可以看做是每个元素与变量均值的差的平方和的均值，即：

$$Var(a) = \frac {1}{m} \sum_{i=1}^{m}(a_i-\mu)^2$$

为了方便处理，我们将每个变量的均值都化为 0 ，因此方差可以直接用每个元素的平方和除以元素个数表示：

$$Var(a) = \frac {1}{m} \sum_{i=1}^{m}(a_i)^2$$

于是上面的如何判断分散问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

协方差

在一维空间中我们可以用方差来表示数据的分散程度。而对于高维数据，我们用协方差进行约束，协方差可以表示两个变量的相关性。为了让两个变量尽可能表示更多的原始信息，我们希望它们之间不存在线性相关性，因为相关性意味着两个变量不是完全独立，必然存在重复表示的信息。

协方差公式为：

$$Cov(a,b)=\frac {1} {m-1} \sum_{i=1}^{m}(a_i-\mu_a)(b_i-\mu_b)$$

因为均值为 0，因此：

$$Cov(a,b)=\frac {1} {m} \sum_{i=1}^{m}a_ib_i$$

当样本数较大时，不必在意其是 m 还是 m-1，为了方便计算，我们分母取 m。

当协方差为 0 时，表示两个变量完全独立。为了让协方差为 0，我们选择第二个基时只能在与第一个基正交的方向上进行选择，因此最终选择的两个方向一定是正交的。

至此，我们得到了降维问题的优化目标：将一组 N 维向量降为 K 维，其目标是选择 K 个单位正交基，使得原始数据变换到这组基上后，各变量两两间协方差为 0，而变量方差则尽可能大（在正交的约束下，取最大的 K 个方差）。

这里的变量不是指某条数据，而是一个表示特征的基向量。因此 PCA 降维的目标就是找 K 个正交基（向量内积为 0 则正交）。

协方差矩阵

针对我们给出的优化目标，接下来我们将从数学的角度来给出优化目标。

我们看到，最终要达到的目的与变量内方差及变量间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。

假设我们只有 a 和 b 两个变量，有 m 个样本，那么我们将它们组成矩阵 X：

$$X= \begin{pmatrix} a_1 & a_2 & \cdots & a_m \\ b_1 & b_2 & \cdots & b_m \end{pmatrix}$$

然后：

$$\frac{1}{m} X X^T= \begin{pmatrix} \frac{1}{m}\sum_{i=1}^{m}a_{i}^{2} & \frac{1}{m}\sum_{i=1}^{m}a_ib_i \\ \frac{1}{m}\sum_{i=1}^{m}a_ib_i & \frac{1}{m}\sum_{i=1}^{m}b_{i}^{2} \end{pmatrix} =\begin{pmatrix} Cov(a,a) & Cov(a,b) \\ Cov(b,a) & Cov(b,b) \\ \end{pmatrix}$$

我们可以看到这个矩阵对角线上的分别是两个变量的方差，而其它元素是 a 和 b 的协方差。两者被统一到了一个矩阵里。

根据上述公式，我们很容易被推广到一般情况：

设我们有 m 个 n 维数据记录，将其排列成矩阵 $X_{n,m}$，设 $C=\frac{1}{m}XX^T$ ，则 C 是一个对称矩阵，其对角线分别对应各个变量的方差，而第 i 行 j 列和 j 行 i 列元素相同，表示 i 和 j 两个维度特征的协方差。

矩阵对角化

根据我们的优化条件，我们需要将除对角线外的其它元素化为 0（协方差为 0），并且在对角线上将元素按大小从上到下排列（方差尽可能大），这样我们就达到了优化目的。这样说可能还不是很明晰，我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系。

设原始数据矩阵 X 对应的协方差矩阵为 C，而 P 是一组基按行组成的矩阵，设 Y=PX，则 Y 为 X 对 P 做基变换后的数据。设 Y 的协方差矩阵为 D，我们推导一下 D 与 C 的关系：

$$\begin{align} D &= \frac{1}{m}YY^T \\ &= \frac{1}{m}(PX)(PX)^T \\ &= \frac{1}{m}PXX^TP^T \\ &= P(\frac{1}{m}XX^T)P^T \\ &= P(C)P^T \\ \end{align}$$

这样我们就看清楚了，我们要找的 P 是能让原始协方差矩阵对角化的 P。换句话说，优化目标变成了寻找一个矩阵 P，满足 $PCP^T$ 是一个对角矩阵，并且对角元素按从大到小依次排列，那么 P 的前 K 行就是要寻找的基，用 P 的前 K 行组成的矩阵乘以 X 就使得 X 从 N 维降到了 K 维并满足上述优化条件。

至此，我们离 PCA 还有仅一步之遥，我们还需要完成对角化。

由上文知道，协方差矩阵 C 是一个是对称矩阵，在线性代数中实对称矩阵有一系列非常好的性质：

1. 实对称矩阵不同特征值对应的特征向量必然正交。
2. 设特征向量 $\lambda$ 重数为 r，则必然存在 r 个线性无关的特征向量对应于 $\lambda$ ，因此可以将这 r 个特征向量单位正交化。

由上面两条可知，一个 n 行 n 列的实对称矩阵一定可以找到 n 个单位正交特征向量，设这 n 个特征向量为：$e_1,e_2,\cdots,e_n$，我们将其按列组成矩阵：$E=(e_1,e_2,\cdots,e_n)$。

则对协方差矩阵 C，有如下结论：

$$E^TCE=\Lambda = \begin{pmatrix} \lambda_1 & & & \\ & \lambda_2 & & \\ & & \cdots & \\ & & & \lambda_n \end{pmatrix}$$

其中 $\Lambda$ 为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

到这里，我们发现我们已经找到了需要的矩阵 P：$P=E^T$。

P 是协方差矩阵的特征向量单位化后按行排列出的矩阵，其中每一行都是 C 的一个特征向量。如果设 P 按照 $\Lambda$ 中特征值的从大到小，将特征向量从上到下排列，则用 P 的前 K 行组成的矩阵乘以原始数据矩阵 X，就得到了我们需要的降维后的数据矩阵 Y。

原文还介绍基于拉格朗日乘子法的 PCA 数学推导，我这里就不展开了。

PCA 求解步骤

总结一下 PCA 的算法步骤：

设有 m 条 n 维数据：

1. 将原始数据按列组成 n 行 m 列矩阵 X；
2. 将 X 的每一行进行零均值化，即减去这一行的均值；
3. 求出协方差矩阵 $C=\frac {1}{m}XX^T$ ；
4. 求出协方差矩阵的特征值及对应的特征向量；
5. 将特征向量按对应特征值大小从上到下按行排列成矩阵，取前 K 行组成矩阵 P；
6. 求得矩阵 $Y=PX$，即为降维到 K 维后的数据。

零均值化

当对训练集进行 PCA 降维时，也需要对验证集、测试集执行同样的降维。而对验证集、测试集执行零均值化操作时，均值必须从训练集计算而来，不能使用验证集或者测试集的中心向量。

其原因也很简单，因为我们的训练集时可观测到的数据，测试集不可观测所以不会知道其均值，而验证集再大部分情况下是在处理完数据后再从训练集中分离出来，一般不会单独处理。如果真的是单独处理了，不能独自求均值的原因是和测试集一样。

另外我们也需要保证一致性，我们拿训练集训练出来的模型用来预测测试集的前提假设就是两者是独立同分布的，如果不能保证一致性的话，会出现 Variance Shift 的问题。

PCA 的一些性质

1. 缓解维度灾难：PCA 算法通过舍去一部分信息之后能使得样本的采样密度增大（因为维数降低了），这是缓解维度灾难的重要手段；
2. 降噪：当数据受到噪声影响时，最小特征值对应的特征向量往往与噪声有关，将它们舍弃能在一定程度上起到降噪的效果；
3. 过拟合：PCA 保留了主要信息，但这个主要信息只是针对训练集的，而且这个主要信息未必是重要信息。有可能舍弃了一些看似无用的信息，但是这些看似无用的信息恰好是重要信息，只是在训练集上没有很大的表现，所以 PCA 也可能加剧了过拟合；
4. 特征独立：PCA 不仅将数据压缩到低维，它也使得降维之后的数据各特征相互独立；

PCA 和 SVD 的对比

这是两个不同的数学定义。我们先给结论：特征值和特征向量是针对方阵才有的，而对任意形状的矩阵都可以做奇异值分解。

PCA

方阵的特征值分解，对于一个方针 A，总可以写成：$A=Q \Lambda Q^{-1}$。其中，Q 是这个矩阵 A 的特征向量组成的矩阵，$\Lambda$ 是一个对角矩阵，每一个对角线元素就是一个特征值，里面的特征值是由小排列的，这些特征值所对应的特征向量就是描述这个矩阵变化方向（从主要的变化到次要的变化排列）。也就是说矩阵 A 的信息可以由其特征值和特征向量表示，这也是矩阵特征值和特征向量的重要物理意义之一。

SVD

矩阵的奇异值分解其实就是对于矩阵 A 的协方差矩阵 $A^TA$ 和 $AA^T$ 做特征值分解推导出来的：

$$A_{m,n}=U_{m,m} \Lambda_{m,n} V_{n,n}^{T} \approx U_{m,k} \Lambda_{k,k} V_{k,n}^{T}$$

其中：U 和 V 都是正交矩阵，有 $U^TU=I_m$，$V^TV=I_n$。这里的约等于是因为 $\Lambda$ 中有 n 个奇异值，但是由于排在后面的很多接近 0，所以我们可以仅保留比较大的 k 个奇异值。

$$A^TA=(U \Lambda V^T)^T U \Lambda V^T=V \Lambda^TU^TU \Lambda V^T=V \Lambda^2V^T$$ $$AA^T=U \Lambda V^T(U \Lambda V^T)^T =U \Lambda V^TV \Lambda^TU^T=U \Lambda^2U^T$$

所以，V U 两个矩阵分别是 $A^TA$ 和 $AA^T$ 的特征向量，中间的矩阵对角线的元素是 $A^TA$ 和 $AA^T$ 的特征值。我们也很容易看出 A 的奇异值和 $A^TA$ 的特征值之间的关系。

PCA 需要对协方差矩阵 $C=\frac {1}{m} XX^T$ 进行特征值分解；SVD 也是对 $A^TA$ 进行特征值分解。如果取 $A=\frac{X^T}{\sqrt{m}}$ ，则二者基本等价。所以 PCA 问题可以转换成 SVD 求解。实际上 Sklearn 的 PCA 就是用 SVD 进行求解的，原因有以下几点：

1. 当样本维度很高时，协方差矩阵计算太慢；
2. 方针特征值分解计算效率不高；
3. SVD 出了特征值分解这种求解方式外，还有更高效更准球的迭代求解方式，避免了 $A^TA$ 的计算；
4. 其实 PCA 与 SVD 的右奇异向量的压缩效果相同。