数据降维之PCA透彻讲解

机器学习

---

0. 引言

　　
　　在做任何机器学习算法训练之前，无论该算法模型是分类、回归还是聚类，我们都需要从原始数据中提取数据特征，然后将提取的特征组合成特征向量的形式输入到模型中进行训练。但是现实中往往会存在这么一个问题：原始数据特征的维度特别高，不利于模型的训练。这时候我们需要对原始数据的特征向量进行降维处理，然后将降维后的数据输入到模型中进行训练，得出最终的训练模型。

　　PCA算法就是在这种需求之下产生的。它是专门用来对高维数据进行降维而设计的，通过将高维数据降维后得到的低维数能加快模型的训练速度，并且低维度的特征具有更好的可视化性质。但是有一点我们必须明确：降维或多或少地会带来一定的信息损失，因此我们在真正PCA算法的应用过程中一般需要设置一个损失阈值的方法来有效的控制信息的损失。

1. PCA算法的数学基础

　　PCA的理论基础建立在严格的线性代数矩阵分析之上。因此，为了让广大的算法爱好者都能够更好的理解PCA算法的工作原理。我们需要先温习一下与PCA算法相关的线性代数基础知识。

1.1 数据的向量表示及其基本运算

　　假设向量 $A=(a_1,a_2,...,a_n)^T$ 和 $B=(b_1,b_2,...,b_n)^T$，则两者的内积表示为如下计算：

$$A \cdot B =a_1b_1 + a_2b_2 + \cdots + a_nb_n = \sum_{i=1}^n{a_i b_i}$$ 内积运算将两个向量映射为一个实数。其计算方式非常容易理解，但是其意义并不明显。下面我们分析内积的几何意义。假设A和B是两个n维向量，我们知道n维向量可以等价表示为n维空间中的一条从原点发射的有向线段，为了简单起见我们假设A和B均为二维向量，则 $A=(x_1,y_1)$，$B=(x_2,y_2)$。则在二维平面上A和B可以用两条发自原点的有向线段表示，见下图：
　　

现在我们从A点向B所在直线引一条垂线。我们知道垂线与B的交点叫做A在B上的投影，再设A与B的夹角是 $\alpha$，则投影的矢量长度为$|A|=\sqrt{x_1^2+y_1^2}$，其中$|A|=\sqrt{x_1^2+y_1^2}$是向量A的模，也就是A线段的标量长度。

注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于0，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。

到这里还是看不出内积和这东西有什么关系，不过如果我们将内积表示为另一种我们熟悉的形式：

$$A\cdot B=|A|\cdot|B|\cdot cos(\alpha)$$

现在事情似乎是有点眉目了：A与B的内积等于A到B的投影长度乘以B的模。再进一步，如果我们假设B的模为1，即让|B|=1|B|=1，那么就变成了：

$$A\cdot B=|A|\cdot cos(\alpha)$$ 也就是说，设向量B的模为1，则A与B的内积值等于A向B所在直线投影的矢量长度！这就是内积的一种几何解释，也是我们得到的第一个重要结论。在后面的推导中，将反复使用这个结论。

1.2 向量空间的基

　　下面我们继续在二维空间内讨论向量。上文说过，一个二维向量可以对应二维笛卡尔直角坐标系中从原点出发的一个有向线段。例如下面这个向量：

　　在代数表示方面，我们经常用线段终点的点坐标表示向量，例如上面的向量可以表示为(3,2)，这是我们再熟悉不过的向量表示。

　　不过我们常常忽略，只有一个(3,2)本身是不能够精确表示一个向量的。我们仔细看一下，这里的3实际表示的是向量在x轴上的投影值是3，在y轴上的投影值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。那么一个向量(3,2)实际是说在x轴投影为3而y轴的投影为2。注意投影是一个矢量，所以可以为负。

更正式的说，向量(x,y)实际上表示线性组合：

$$x(1,0)^\mathsf{T}+y(0,1)^\mathsf{T}$$ 不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基。

　　所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。

　　我们之所以默认选择(1,0)和(0,1)为基，当然是比较方便，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。

　　例如，(1,1)和(-1,1)也可以成为一组基。一般来说，我们希望基的模是1，因为从内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基而直接获得其在新基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基可以变为$(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$和$(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$。

　　现在，我们想获得(3,2)在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，我们只要分别计算(3,2)和两个基的内积，不难得到新的坐标为$(\frac{5}{\sqrt2},-\frac{1}{\sqrt{2}})$。下图给出了新的基以及(3,2)在新基上坐标值的示意图：

　　另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。

1.3 基变换的矩阵表示

　　下面我们找一种简便的方式来表示基变换。还是拿上面的例子，想一下，将(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}$$

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

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

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

2. PCA算法的数学原理

　　上面我们讨论了选择不同的基可以对同样一组数据给出不同的表示，而且如果基的数量少于向量本身的维数，则可以达到降维的效果。但是我们还没有回答一个最最关键的问题：如何选择基才是最优的。或者说，如果我们有一组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轴都不是最好的投影选择。我们直观目测，如果向通过第一象限和第三象限的斜线投影，则五个点在投影后还是可以区分的。

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

2.1 方差

　　上文说到，我们希望投影后投影值尽可能分散，而这种分散程度，可以用数学上的方差来表述。此处，一个字段的方差可以看做是每个元素与字段均值的差的平方和的均值，即：

$$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}$$ 于是上面的问题被形式化表述为：寻找一个一维基，使得所有数据变换为这个基上的坐标表示后，方差值最大。

2.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个方差）。

2.3 协方差矩阵

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

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

假设我们只有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=\frac{1}{m}XX^\mathsf{T}$，则C是一个对称矩阵，其对角线分别个各个字段的方差，而第i行j列和j行i列元素相同，表示i和j两个字段的协方差。

2.4 协方差矩阵对角化

　　根据上述推导，我们发现要达到优化目前，等价于将协方差矩阵对角化：即除对角线外的其它元素化为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^\mathsf{T}$ 是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。

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

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

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

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

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

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

$$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}$$ 其中$Λ$为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。

　　关于如何证明以上结论以及如何求协方差矩阵 $C$ 的n个正交的单位特征向量和及其对应的特征值的过程我们会在将来单独的一篇文章《实对称矩阵单位正交特征向量及其对应的特征值求解过程剖析》来详细介绍。

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

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

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

3. PCA算法执行过程及其损失误差分析

3.1 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维后的数据

3.2 PCA算法的损失误差分析

　　由上述的分析可知，在PCA算法的执行过程中，计算出了协方差矩阵的n个特征向量。但是我们只取了最大的前k个特征值对应的特征向量，而将其它的n-k个特征向量直接去除：

$$\begin{align} X & = [X_1,X_2,\cdots,X_k,X_{k+1},\cdots,X_n]^T \\ 　& \approx[X_1,X_2,\cdots,X_k,X_{k+1}]^T \\ 　& = \bar X \end{align}$$

　　而这样直接将后边的n-k个特征去掉会导致一定的信息损失。信息的损失随着k的取值的增大而减少，反之信息的损失就越多。由于后边的特征值非常小，他们对应的信息量也就越少，因此将 特征值接近0的特征去除不会导致过多的损失。
　　
　　但是现在问题来了，如何衡量前k个特征所占重要性（也就是信息量）的大小呢？具体的计算方法为：计算前k个特征值在所有特征值中所占的比例。 由于协方差矩阵的特征值为方差，因而特征值之比即为方差的百分比，具体计算如下：

$$\alpha=\frac{\sum_{j=1}^k\lambda_j}{\sum_{j=1}^n\lambda_j}$$ 其中，$\lambda_j=\sum_{i=1}^m{({X_j}^{(i)})}^2$。
　　
　　一般情况下，我们只需要满足 $\alpha\ge 0.95$ 就可以了，当然具体情况具体分析。另外，PCA算法具有缩放不变性，即所有的特征分量被放大或缩小相同的倍数，PCA输出的特征向量不会发生变化[2]。

4. PCA算法应用实例

4.2 完整代码和输出结果

import matplotlib.pyplot as plt
from sklearn import datasets
from sklearn.decomposition import PCA
from sklearn.discriminant_analysis import LinearDiscriminantAnalysis
iris = datasets.load_iris()
X = iris.data
y = iris.target
target_names = iris.target_names
pca = PCA(n_components=2)
X_r = pca.fit(X).transform(X)
lda = LinearDiscriminantAnalysis(n_components=2)
X_r2 = lda.fit(X, y).transform(X)
# Percentage of variance explained for each components
print('explained variance ratio (first two components): %s'
      % str(pca.explained_variance_ratio_))
plt.figure()
colors = ['navy', 'turquoise', 'darkorange']
lw = 2
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
    plt.scatter(X_r[y == i, 0], X_r[y == i, 1], color=color, alpha=.8, lw=lw,
                label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('PCA of IRIS dataset')
plt.figure()
for color, i, target_name in zip(colors, [0, 1, 2], target_names):
    plt.scatter(X_r2[y == i, 0], X_r2[y == i, 1], alpha=.8, color=color,
                label=target_name)
plt.legend(loc='best', shadow=False, scatterpoints=1)
plt.title('LDA of IRIS dataset')
plt.show()

输出结果：explained variance ratio (first two components): [ 0.92461621 0.05301557]

5. 总结

　　
　　PCA算法非常巧妙地利用协方差矩阵来计算出样本集在不同方向上的分散程度，利用方差最大的方向作为样本集的主方向。其主要过程是：首先利用样本集及特征构建一个样本矩阵，然后利用样本矩阵计算得到协方差矩阵，再计算协方差矩阵的特征值和特征向量，保留特征值前k大的特征向量作为新的维度方向。再将原始样本数据转换到新的空间维度。
　　
　　PCA本质上是将方差最大的方向作为主要特征，并且在各个正交方向上将数据“离相关”，也就是让它们在不同正交方向上没有相关性。

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

---

参考文献：
　
　1. PCA主成分分析方法
　2. PCA的数学原理
　3. 主成分分析（PCA）原理详解

---

对机器学习和人工智能感兴趣，请扫码关注微信公众号！