一文彻底搞懂PCA

时间 2021-01-14 标签 pca降维数据挖掘机器学习

一文彻底搞懂PCA

算法及实例

简介

PCA（Principal Component Analysis），即主成分分析，是一种基于线性代数的维归约技术，常用于数据挖掘中对数据的预处理，主要用于解决“维灾难”问题。
在阅读本文前，先介绍下文章会涉及到的一些数学理论，这些理论都是PCA技术所依赖的工具或理论基础。

向量相关，如转置、内积、正交
矩阵相关，如矩阵的乘法、对称矩阵及其对角化
数理统计相关，如方差、协方差、相关性

PCA利用线性代数技术，将数据由高维空间投影到低维空间。通俗来讲，是从原数据中，找到新的属性（主成分），这些属性是原属性的线性组合，是相互正交的（文中会解释为何是正交的）。

约定几个符号

n表示数据的特征数量，一个n维向量表示一条具有n个特征的数据
m表示数据的数量
m*n表示一个具有n个特征的、m行数据的矩阵
右上角大写的T表示转置

PCA的应用场景

在实际工作中，我们面对的数据具有大量的特征，但现在我们不妨以少数几个特征为例，足以展开研究。假如我们有一组房价的数据如下

面积（平方英尺）	面积（平方米）	卧室数量	浴室数量	阳台数量
1291	120	4	2	3
947	88	2	1	1

假如我们现在面临的“维灾难”问题是，我们的计算资源仅足够计算3个特征，而现在有5个特征。
必须承认，前两个表示面积的特征是线性相关的，即1 平方米=10.7639104 平方英尺。对于这种线性相关的特征，我们可以任取其一（不妨去掉第一个特征），于是5个特征还有4个。
再看后面三个特征，不难相信，通常越大的房屋面积，也将具有越多的卧室、浴室、阳台数量，但显然后三者和房屋面积不具有直接的线性关系。也就是说，我们知道这几个特征是相关，但又不是线性相关，无法简单的替换或删除。
我们应该如何做，能够从4个特征中保留3个，而不丢失信息？此时，就轮到PCA上场了。

PCA的数学原理

现在，为了由浅入深的理解PCA，不妨抛开PCA，先了解以下几个内容：

数据的向量表示

上文中的（第一行）数据，在去掉第一个特征后，通常用列向量表示为：
$(120,4,2,3)^T$
这里使用了转置，主要是因为横着书写一个行向量比较方便，而行向量的转置就是列向量。
现在不妨设两个向量A、B。我们知道，两个（列）向量的内积
$[A,B]=A^TB=(a_1,a_2,...,a_n)^T*(b_1,b_2,...,b_n)=a_1b_1+a_2b_2+...+a_nb_n$
内积运算的结果为一个实数，其计算方式虽然简单，但意义并不明显。我们知道n维向量可以等价表示为n维空间中的一条从原点发射的有向线段，为了便于理解，我们不妨以二维的向量为例：
$A=(x_1,y_1)，B=(x_2,y_2)$

则[A,B]的几何直观解释是A在B上的投影乘以B的长度，几何公式为
$|A|*|B|*cos(a)$ ，其中|A|表示A的模，即标量长度；代数公式为 $||A|| * ||B||*cos(a)$ ，其中||A||为A的范数。

现在，不妨我们让向量B的模始终为1，就有
$A⋅B=|A|*|B|cos(a)==|A|cos(a)$
也就是说，设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度，后文中我们将会继续用到这个工具。

向量的基

通常我们将一个二维向量描述成从原点出发的有向线段，并且直接使用线段终点的点坐标来表示。

现在，我们不妨换一种理解方式，首先要引入基的概念，或者说是单位向量。
在二维直角坐标系中，我们引入两个模为1的向量i（1，0）、j（0，1），这样一对特殊的向量也被成为单位向量，也就是二维空间中的一组基。

此时的向量，则理解成为两个方向的基的缩放，即
$(3,-2)=3*i+(-2)*j$

不难证明，对任意向量(x,y)，都有 $(x,y)=x*i+y*j$ 。而x、y，正是向量分别在两个基上的缩放值。
值得注意的是，例子中给出的基是正交的（即内积为0，或直观说相互垂直）单位向量，但可以成为一组基的唯一要求就是线性无关，任意长度、非正交的基也是可以的（如下图）。不过因为正交基有较好的性质，所以一般使用的基都是正交单位向量。

上图中的v和w同样更够表达二维空间中的任意向量，也可以作为一组基。因此，当我们将一组标量(x,y)作用到不同的基上，显然会得到不一样的向量。

基变换的矩阵表示

下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(3,2)变换为新基上的坐标，就是用(3,2)与第一个基做内积运算，作为第一个新的坐标分量，然后用(3,2)与第二个基做内积运算，作为第二个新坐标的分量。实际上，我们可以用矩阵相乘的形式简洁的表示这个变换：

$\begin{pmatrix}1/\sqrt{2}&1/\sqrt{2} \\-1/\sqrt{2}&1/\sqrt{2}\end{pmatrix}\begin{pmatrix}3\\2\end{pmatrix}=\begin{pmatrix} 5/\sqrt{2}\\-1/\sqrt{2}\end{pmatrix}$

太漂亮了！其中矩阵的两行分别为两个基，乘以原向量，其结果刚好为新基的坐标。可以稍微推广一下，如果我们有m个二维向量，只要将二维向量按列排成一个两行m列矩阵，然后用“基矩阵”乘以这个矩阵，就得到了所有这些向量在新基下的值。例如(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：

$\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} &1/\sqrt{2}\end{pmatrix}\begin{pmatrix} 1 & 2 & 3 \\ 1 & 2 & 3\end{pmatrix} =\begin{pmatrix} 2/\sqrt{2} & 4/\sqrt{2} & 6/\sqrt{2} \\0&0&0\end{pmatrix}$
于是一组向量的基变换被干净的表示为矩阵的相乘。

一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。

数学表示为：
$\begin{pmatrix} p_1 \\ p_2 \\ \vdots \\ p_R \end{pmatrix} \begin{pmatrix} a_1 & a_2 & \cdots & a_M\end{pmatrix}=\begin{pmatrix} p_1a_1 & p_1a_2 & \cdots & p_1a_M \\ p_2a_1 & p_2a_2 & \cdots & p_2a_M \\ \vdots & \vdots & \ddots & \vdots \\ p_Ra_1 & p_Ra_2 & \cdots & p_Ra_M \end{pmatrix}$

其中pi是一个行向量，表示第i个基，aj是一个列向量，表示第j个原始数据记录。

特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换。

最后，上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。更抽象的说，一个矩阵可以表示一种线性变换。很多同学在学线性代数时对矩阵相乘的方法感到奇怪，但是如果明白了矩阵相乘的物理意义，其合理性就一目了然了。

协方差矩阵及优化目标

上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组N维向量，现在要将其降到K维（K小于N），那么我们应该如何选择K个基才能最大程度保留原有的信息？

要完全数学化这个问题非常繁杂，这里我们用一种非形式化的直观方法来看这个问题。

为了避免过于抽象的讨论，我们仍以一个具体的例子展开。假设我们的数据由五条记录组成，将它们表示成矩阵形式：

$\begin{pmatrix} 1 & 1 & 2 & 4 & 2 \\ 1 & 3 & 3 & 4 & 4 \end{pmatrix}$

其中每一列为一条数据记录，而一行为一个字段。为了后续处理方便，我们首先将每个字段内所有值都减去字段均值，其结果是将每个字段都变为均值为0（这样做的道理和好处后面会看到）。

我们看上面的数据，第一个字段均值为2，第二个字段均值为3，所以变换后：

$\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}$

我们可以看下五条数据在平面直角坐标系内的样子：

现在问题来了：如果我们必须使用一维来表示这些数据，又希望尽量保留原始的信息，你要如何选择？

通过上一节对基变换的讨论我们知道，这个问题实际上是要在二维平面中选择一个方向，将所有数据都投影到这个方向所在直线上，用投影值表示原始记录。这是一个实际的二维降到一维的问题。

那么如何选择这个方向（或者说基）才能尽量保留最多的原始信息呢？一种直观的看法是：希望投影后的投影值尽可能分散。

以上图为例，可以看出如果向x轴投影，那么最左边的两个点会重叠在一起，中间的两个点也会重叠在一起，于是本身四个各不相同的二维点投影后只剩下两个不同的值了，这是一种严重的信息丢失，同理，如果向y轴投影最上面的两个点和分布在x轴上的两个点也会重叠。所以看来x和y轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。

下面，我们用数学方法表述这个问题。

方差

上文说到，我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：
$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}$
于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

协方差

对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。

如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。

数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

$Cov(a,b)=\frac{1}{m}\sum_{i=1}^m{a_ib_i}$
可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。

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

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

协方差矩阵

上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。

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

假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：

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

然后我们用X乘以X的转置，并乘上系数1/m：

$\frac{1}{m}XX^\mathsf{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}$

奇迹出现了！这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。

根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：

设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设C=1mXX𝖳，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差。

协方差矩阵对角化

根据上述推导，我们发现要达到优化目前，等价于将协方差矩阵对角化：即除对角线外的其它元素化为0，并且在对角线上将元素按大小从上到下排列，这样我们就达到了优化目的。这样说可能还不是很明晰，我们进一步看下原矩阵与基变换后矩阵协方差矩阵的关系：

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

$\begin{array}{l l l} D & = & \frac{1}{m}YY^\mathsf{T} \\ & = & \frac{1}{m}(PX)(PX)^\mathsf{T} \\ & = & \frac{1}{m}PXX^\mathsf{T}P^\mathsf{T} \\ & = & P(\frac{1}{m}XX^\mathsf{T})P^\mathsf{T} \\ & = & PCP^\mathsf{T} \end{array}$

现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足PCP𝖳是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

至此，我们离“发明”PCA还有仅一步之遥！

现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。

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

1）实对称矩阵不同特征值对应的特征向量必然正交。

2）设特征向量λ重数为r，则必然存在r个线性无关的特征向量对应于λ，因此可以将这r个特征向量单位正交化。

由上面两条可知，一个n行n列的实对称矩阵一定可以找到n个单位正交特征向量，设这n个特征向量为 $e1,e2,⋯,en$ ，我们将其按列组成矩阵：

$E=\begin{pmatrix} e_1 & e_2 & \cdots & e_n \end{pmatrix}$
则对协方差矩阵C有如下结论：

$E^\mathsf{T}CE=\Lambda=\begin{pmatrix} \lambda_1 & & & \\ & \lambda_2 & & \\ & & \ddots & \\ & & & \lambda_n \end{pmatrix}$

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

以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。

到这里，我们发现我们已经找到了需要的矩阵P：

$P=E^\mathsf{T}$

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

至此我们完成了整个PCA的数学原理讨论。在下面的一节，我们将给出PCA的一个实例。

算法及实例

为了巩固上面的理论，我们在这一节给出一个具体的PCA实例。

PCA算法
总结一下PCA的算法步骤：

设有m条n维数据。

1）将原始数据按列组成n行m列矩阵X
2）将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值
3）求出协方差矩阵 $C=\frac{1}{m}XX^\mathsf{T}$
4）求出协方差矩阵的特征值及对应的特征向量
5）将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P

6） $Y=PX$ 即为降维到k维后的数据

实例

这里以上文提到的

$\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}$
为例，我们用PCA方法将这组二维数据其降到一维。

因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：

$C=\frac{1}{5}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}\begin{pmatrix} -1 & -2 \\ -1 & 0 \\ 0 & 0 \\ 2 & 1 \\ 0 & 1 \end{pmatrix}=\begin{pmatrix} \frac{6}{5} & \frac{4}{5} \\ \frac{4}{5} & \frac{6}{5} \end{pmatrix}$
然后求其特征值和特征向量，具体求解方法不再详述，可以参考相关资料。求解后特征值为：

$\lambda_1=2,\lambda_2=2/5$
其对应的特征向量分别是：

$c_1\begin{pmatrix} 1 \\ 1 \end{pmatrix},c_2\begin{pmatrix} -1 \\ 1 \end{pmatrix}$

其中对应的特征向量分别是一个通解， $c1$ 和 $c2$ 可取任意实数。那么标准化后的特征向量为：
$\begin{pmatrix} 1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix},\begin{pmatrix} -1/\sqrt{2} \\ 1/\sqrt{2} \end{pmatrix}$
因此我们的矩阵P是：

$P=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}$

可以验证协方差矩阵C的对角化：

$PCP^\mathsf{T}=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \\ -1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} 6/5 & 4/5 \\ 4/5 & 6/5 \end{pmatrix}\begin{pmatrix} 1/\sqrt{2} & -1/\sqrt{2} \\ 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}=\begin{pmatrix} 2 & 0 \\ 0 & 2/5 \end{pmatrix}$

最后我们用P的第一行乘以数据矩阵，就得到了降维后的表示：

$Y=\begin{pmatrix} 1/\sqrt{2} & 1/\sqrt{2} \end{pmatrix}\begin{pmatrix} -1 & -1 & 0 & 2 & 0 \\ -2 & 0 & 0 & 1 & 1 \end{pmatrix}=\begin{pmatrix} -3/\sqrt{2} & -1/\sqrt{2} & 0 & 3/\sqrt{2} & -1/\sqrt{2} \end{pmatrix}$

降维投影结果如下图：

进一步讨论

根据上面对PCA的数学原理的解释，我们可以了解到一些PCA的能力和限制。PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。

因此，PCA也存在一些限制，例如它可以很好的解除线性相关，但是对于高阶相关性就没有办法了，对于存在高阶相关性的数据，可以考虑Kernel PCA，通过Kernel函数将非线性相关转为线性相关，关于这点就不展开讨论了。另外，PCA假设数据各主特征是分布在正交方向上，如果在非正交方向上存在几个方差较大的方向，PCA的效果就大打折扣了。

最后需要说明的是，PCA是一种无参数技术，也就是说面对同样的数据，如果不考虑清洗，谁来做结果都一样，没有主观参数的介入，所以PCA便于通用实现，但是本身无法个性化的优化。

希望这篇文章能帮助朋友们了解PCA的数学理论基础和实现原理，借此了解PCA的适用场景和限制，从而更好的使用这个算法。

记得关注我的微信公众号

注：部分内容引用其他文章或视频截图，侵删。