7.1-PCA

概述

  PCA有两种通俗易懂的解:
1. 是最大化投影后数据的方差(让数据更分散)；
2. 是最小化投影造成的损失。

  这两个思路最后都能推导出同样的结果。
  下图应该是对PCA第二种解释展示得最好的一张图片了.

原理

  在理解特征提取与处理时，涉及高维特征向量的问题往往容易陷入维度灾难。随着数据集维度的增加，算法学习需要的样本数量呈指数级增加。有些应用中，遇到这样的大数据是非常不利的，而且从大数据集中学习需要更多的内存和处理能力。另外，随着维度的增加，数据的稀疏性会越来越高。在高维向量空间中探索同样的数据集比在同样稀疏的数据集中探索更加困难。
  PCA将数据投射到一个低维子空间实现降维。例如，二维数据集降维就是把点投射成一条线，数据集的每个样本都可以用一个值表示，不需要两个值。三维数据集可以降成二维，就是把变量映射成一个平面。一般情况下，n维数据集可以通过映射降成k维子空间，其中k≤n。
假如你是一本养花工具宣传册的摄影师，你正在拍摄一个水壶。水壶是三维的，但是照片是二维的，为了更全面的把水壶展示给客户，你需要从不同角度拍几张图片。下图是你从四个方向拍的照片：

  第一张图里水壶的背面可以看到，但是看不到前面。第二张图是拍前面，可以看到壶嘴，这张图可以提供了第一张图缺失的信息，但是壶把看不到了。从第三张俯视图里无法看出壶的高度。第四张图是你真正想要的，水壶的高度，顶部，壶嘴和壶把都清晰可见。
PCA的设计理念与此类似，它可以将高维数据集映射到低维空间的同时，尽可能的保留更多信息。PCA旋转数据集与其主成分对齐，将最多的信息保留到第一主成分中。假设我们有下图所示的数据集：

  数据集看起来像一个从原点到右上角延伸的细长扁平的椭圆。要降低整个数据集的维度，我们必须把点映射成一条线。下图中的两条线都是数据集可以映射的，映射到哪条线样本变化最大？

  显然，样本映射到黑色虚线的变化比映射到红色点线的变化要大的多。实际上，这条黑色虚线就是第一主成分。第二主成分必须与第一主成分正交，也就是说第二主成分必须是在统计学上独立的，会出现在与第一主成分垂直的方向，如下图所示：

  后面的每个主成分也会尽量多的保留剩下的信息，唯一的要求就是每一个主成分需要和前面的主成分正交。
现在假设数据集是三维的，散点图看起来像是沿着一个轴旋转的圆盘。

  这些点可以通过旋转和变换使圆盘完全变成二维的。现在这些点看着像一个椭圆，第三维上基本没有信息，可以被忽略。
当数据集不同维度上的方差分布不均匀的时候，PCA最有用。（如果是一个球壳行数据集，PCA不能有效的发挥作用，因为各个方向上的方差都相等；没有丢失大量的信息维度一个都不能忽略）。

理论依据

  首先看一个例子：这里给出5个同学的成绩

header 1	学生1	学生2	学生3	学生4	学生5
语文成绩	50	60	70	80	90
数学成绩	80	82	84	86	88

  此时如果将二维特征【语文成绩，数学成绩】压缩成1维，仅仅用数学成绩表达似乎不是最佳方案，即使沿着投影方向2，只用语文成绩表达，其分散程度似乎也比不上投影方向1的降维结果。
  如何找到最佳的（方差最大，即最分散。这也容易理解，毕竟分散的数据更容易区分开来）降维方式呢？显然最佳方案不总是用语文/数学成绩单一替代，这便需要本文的主角——PCA算法。简而言之：PCA算法其表现形式是降维，同时也是一种特征融合算法。
  推导之前先介绍几个操作的几何意义：

1. 内积和投影

  下面先来看一个高中就学过的向量运算：内积。两个维数相同的向量的内积被定义为：

(a1,a2,⋯,an)T⋅(b1,b2,⋯,bn)T=a1b1+a2b2+⋯+anbn $(a_1,a_2,\cdots ,a_n{)}^T\cdot (b_1,b_2,\cdots ,b_n{)}^T=a_1{b}_1+a_2{b}_2+\cdots +a_n{b}_n$

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

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

注意这里我们专门区分了矢量长度和标量长度，标量长度总是大于等于0，值就是线段的长度；而矢量长度可能为负，其绝对值是线段长度，而符号取决于其方向与标准方向相同或相反。
  到这里还是看不出内积和这东西有什么关系，不过如果我们将内积表示为另一种我们熟悉的形式：

A⋅B=|A||B|cos(a) $A\cdot B=|A||B|cos(a)$

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

A⋅B=|A|cos(a) $A\cdot B=|A|cos(a)$

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

基的概念

给出下面这个向量：

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

x(1,0)T+y(0,1)T $x(1,0{)}^T+y(0,1{)}^T$

  不难证明所有二维向量都可以表示为这样的线性组合。此处(1,0)和(0,1)叫做二维空间中的一组基。

  所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，就可以了。只不过我们经常省略第一步，而默认以(1,0)和(0,1)为基。基就是骨架，基的个数即为撑起给定维度空间的最小向量个数，通常设定基德模为1。

基变换的矩阵表示

  但基不是唯一的，例如上一个例子中选(1,0)和(0,1)为基，同样，上面的基可以变为 (12√,12√) $(\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$和 (−12√,12√) $(-\frac{1}{\sqrt{2}},\frac{1}{\sqrt{2}})$.
  现在我们想获得(3,2)在新基上的坐标，即在两个方向上的投影矢量值，就需要利用到基变换的矩阵表示。
  根据下图，我们只要分别计算(3,2)和两个基的内积，不难得到新的坐标为 (52√,−12√) $(\frac{5}{\sqrt{2}},-\frac{1}{\sqrt{2}})$.

  同样我们可以利用基变换的矩阵表示进行求解：

  再举一例：(1,1)，(2,2)，(3,3)，想变换到刚才那组基上，则可以这样表示：

  于是一组向量的基变换被干净的表示为矩阵的相乘。
  推广成一般性结论：
  一般的，如果我们有M个N维向量，想将其变换为由R个N维向量表示的新空间中，那么首先将R个基按行组成矩阵A，然后将向量按列组成矩阵B，那么两矩阵的乘积AB就是变换结果，其中AB的第m列为A中第m列变换后的结果。
  对应数学表达：

  其中 pi $p_i$是一个行向量，表示第i个基， aj $a_j$是一个列向量，表示第j个原始数据记录。
  上述分析同时给矩阵相乘找到了一种物理解释：两个矩阵相乘的意义是将右边矩阵中的每一列列向量变换到左边矩阵中每一行行向量为基所表示的空间中去。这种操作便是降维。

最优基的选取方法

方差

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

Var(a)=1m∑mi=1(ai−μ)2 $Var(a)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}(a_i-\mu {)}^2$

  由于上面我们已经将每个字段的均值都化为0了，因此方差可以直接用每个元素的平方和除以元素个数表示:

Var(a)=1m∑mi=1a2i $Var(a)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i^2$

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

协方差

  对于上面二维降成一维的问题来说，找到那个使得方差最大的方向就可以了。不过对于更高维，还有一个问题需要解决。考虑三维降到二维问题。与之前相同，首先我们希望找到一个方向使得投影后方差最大，这样就完成了第一个方向的选择，继而我们选择第二个投影方向。
  如果我们还是单纯只选择方差最大的方向，很明显，这个方向与第一个方向应该是“几乎重合在一起”，显然这样的维度是没有用的，因此，应该有其他约束条件。从直观上说，让两个字段尽可能表示更多的原始信息，我们是不希望它们之间存在（线性）相关性的，因为相关性意味着两个字段不是完全独立，必然存在重复表示的信息。
  数学上可以用两个字段的协方差表示其相关性，由于已经让每个字段均值为0，则：

Cov(a,b)=1m∑mi=1aibi $Cov(a,b)=\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i$

  可以看到，在字段均值为0的情况下，两个字段的协方差简洁的表示为其内积除以元素数m。
  当协方差为0时，表示两个字段完全独立。为了让协方差为0，我们选择第二个基时只能在与第一个基正交的方向上选择。因此最终选择的两个方向一定是正交的。
  至此，我们得到了降维问题的优化目标：将一组N维向量降为K维（K大于0，小于N），其目标是选择K个单位（模为1）正交基，使得原始数据变换到这组基上后，各字段两两间协方差为0，而字段的方差则尽可能大（在正交的约束下，取最大的K个方差）。

协方差矩阵

  上面我们导出了优化目标，但是这个目标似乎不能直接作为操作指南（或者说算法），因为它只说要什么，但根本没有说怎么做。所以我们要继续在数学上研究计算方案。
  我们看到，最终要达到的目的与字段内方差及字段间协方差有密切关系。因此我们希望能将两者统一表示，仔细观察发现，两者均可以表示为内积的形式，而内积又与矩阵相乘密切相关。于是我们来了灵感：
  假设我们只有a和b两个字段，那么我们将它们按行组成矩阵X：

X=(a1b1a2b2⋯⋯ambm) $X=\left(\begin{array}{cccc}{a}_1& a_2& \cdots & a_m\\ b_1& b_2& \cdots & b_m\end{array}\right)$

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

1mXXT=(1m∑mi=1a2i1m∑mi=1aibi1m∑mi=1aibi1m∑mi=1b2i) $\frac{1}{m}X{X}^T=\left(\begin{array}{cc}\frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i^2& \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i\\ \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{a}_i{b}_i& \frac{1}{m}\underset{i=1}{\overset{m}{\sum }}{b}_i^2\end{array}\right)$

  奇迹出现了！这个矩阵对角线上的两个元素分别是两个字段的方差，而其它元素是a和b的协方差。两者被统一到了一个矩阵的。
  根据矩阵相乘的运算法则，这个结论很容易被推广到一般情况：
  设我们有m个n维数据记录，将其按列排成n乘m的矩阵X，设 C=1mXXT $C=\frac{1}{m}X{X}^T$，则C是一个对称阵，其对角线分别为各个字段的方差，而第i行和j列的元素和第j行i列的元素相同，表示i和j两个字段的协方差。

协方差矩阵对角化

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

D=====1mYYT1m(PX)(PX)T1mPXXTPTP(1mXXT)PTPCPT $\begin{array}{lll}D& =& \frac{1}{m}Y{Y}^T\\ & =& \frac{1}{m}(PX)(PX{)}^T\\ & =& \frac{1}{m}PX{X}^T{P}^T\\ & =& P(\frac{1}{m}X{X}^T)P^T\\ & =& PC{P}^T\end{array}$

  现在事情很明白了！我们要找的P不是别的，而是能让原始协方差矩阵对角化的P。换句话说，优化目标变成了寻找一个矩阵P，满足 PCPT $PC{P}^T$是一个对角矩阵，并且对角元素按从大到小依次排列，那么P的前K行就是要寻找的基，用P的前K行组成的矩阵乘以X就使得X从N维降到了K维并满足上述优化条件。
  至此，我们离“发明”PCA还有仅一步之遥！
  现在所有焦点都聚焦在了协方差矩阵对角化问题上，有时，我们真应该感谢数学家的先行，因为矩阵对角化在线性代数领域已经属于被玩烂了的东西，所以这在数学上根本不是问题。
  由上文知道，协方差矩阵C是一个是对称矩阵，在线性代数上，实对称矩阵有一系列非常好的性质：
1. 实对称矩阵不同特征值对应的特征向量必然正交。
2. 设特征向量 λ $\lambda$重数为r，则必然存在r个线性无关的特征向量对应于 λ $\lambda$，因此可以将这r个特征向量单位正交化.

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

E=(e1e2⋯en) $E=\left(\begin{array}{cccc}{e}_1& e_2& \cdots & e_n\end{array}\right)$

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

ETCE=Λ=⎛⎝⎜⎜⎜⎜λ1λ2⋱λn⎞⎠⎟⎟⎟⎟ $E^{T}CE=\mathrm{\Lambda }=\left(\begin{array}{cccc}{\lambda }_1& & & \\ & {\lambda }_2& & \\ & & \ddots & \\ & & & {\lambda }_n\end{array}\right)$

  其中 Λ $\mathrm{\Lambda }$为对角矩阵，其对角元素为各特征向量对应的特征值（可能有重复）。
  以上结论不再给出严格的数学证明，对证明感兴趣的朋友可以参考线性代数书籍关于“实对称矩阵对角化”的内容。
  到这里，我们发现我们已经找到了需要的矩阵P：

P=ET $P=E^T$

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

算法

  总结一下PCA的算法步骤：
  设有m条n维数据。
1. 将原始数据按列组成n行m列矩阵X.
2. 将X的每一行（代表一个属性字段）进行零均值化，即减去这一行的均值.
3. 求出协方差矩阵 C=1mXXT $C=\frac{1}{m}X{X}^T$
4. 求出协方差矩阵的特征值及对应的特征向量.
5. 将特征向量按对应特征值大小从上到下按行排列成矩阵，取前k行组成矩阵P.
6. Y=PX $Y=PX$即为降维到k维后的数据.

实例

A-示例1：二维数据降维

  这里以上文提到的

(−1−2−10002101) $\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)$

  为例，我们用PCA方法将这组二维数据其降到一维。
  因为这个矩阵的每行已经是零均值，这里我们直接求协方差矩阵：

C=15(−1−2−10002101)⎛⎝⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜⎜−1−1020−20011⎞⎠⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟⎟=⎛⎝⎜⎜⎜65454565⎞⎠⎟⎟⎟ $C=\frac{1}{5}\left(\begin{array}{ccccc}-1& -1& 0& 2& 0\\ -2& 0& 0& 1& 1\end{array}\right)\left(\begin{array}{cc}-1& -2\\ \\ -1& 0\\ \\ 0& 0\\ \\ 2& 1\\ \\ 0& 1\end{array}\right)=\left(\begin{array}{cc}\frac{6}{5}& \frac{4}{5}\\ \\ \frac{4}{5}& \frac{6}{5}\end{array}\right)$

  然后求其特征值和特征向量，具体求解方法不再详述，可以参考相关资料。求解后特征值为：

λ1=2,λ2=2/5 ${\lambda }_1=2,{\lambda }_2=2/5$

  其对应的特征向量分别是：

c1(11),c2⎛⎝⎜⎜−1111),c2⎛⎝⎜⎜−11⎞⎠⎟⎟ $c_1\left(\begin{array}{c}1\\ 1\end{array}\right),c_2\left(\begin{array}{c}-1\\ \\ 1\end{array}\right)$

  其中对应的特征向量分别是一个通解，c1和c2可取任意实数。那么标准化后的特征向量为：

(1/⎛⎝⎜⎜−11⎞⎠⎟⎟ $c_1\left(\begin{array}{c}1\\ 1\end{array}\right),c_2\left(\begin{array}{c}-1\\ \\ 1\end{array}\right)$

  其中对应的特征向量分别是一个通解，c1和c2可取任意实数。那么标准化后的特征向量为：

(1/2‾√1⎞⎠⎟⎟