机器学习算法笔记系列之深入理解主成分分析PCA

雪柳花明 2017-08-22

展开全文

Author: shizhixin
Blog: http://blog.csdn.net/shizhixin
Weibo：http://weibo.com/zhixinshi
Email: zstarstone@163.com
Date: 2016-04-14
Note: 本笔记是本人根据网上和书本资料总结，并对所有作者表示感谢，特别感谢致谢文献中提到的作者！本笔记仅供学习交流！

1 概述

真实的训练数据总是存在各种各样的问题：

1、比如拿到一个汽车的样本，里面既有以“千米/每小时”度量的最大速度特征，也有“英里/小时”的最大速度特征，显然这两个特征有一个多余。

2、拿到一个数学系的本科生期末考试成绩单，里面有三列，一列是对数学的兴趣程度，一列是复习时间，还有一列是考试成绩。我们知道要学好数学，需要有浓厚的兴趣，所以第二项与第一项强相关，第三项和第二项也是强相关。那是不是可以合并第一项和第二项呢？

3、拿到一个样本，特征非常多，而样例特别少，这样用回归去直接拟合非常困难，容易过度拟合。比如北京的房价：假设房子的特征是（大小、位置、朝向、是否学区房、建造年代、是否二手、层数、所在层数），搞了这么多特征，结果只有不到十个房子的样例。要拟合房子特征->房价的这么多特征，就会造成过度拟合。

4、这个与第二个有点类似，假设在IR中我们建立的文档-词项矩阵中，有两个词项为“learn”和“study”，在传统的向量空间模型中，认为两者独立。然而从语义的角度来讲，两者是相似的，而且两者出现频率也类似，是不是可以合成为一个特征呢？

5、在信号传输过程中，由于信道不是理想的，信道另一端收到的信号会有噪音扰动，那么怎么滤去这些噪音呢？

而这里的特征很多是和类标签有关的，但里面存在噪声或者冗余。在这种情况下，需要一种特征降维的方法来减少特征数，减少噪音和冗余，减少过度拟合的可能性。

下面探讨一种称作主成分分析（PCA）的方法来解决部分上述问题。PCA的思想是将n维特征映射到d维上(d<n)，这d维是全新的正交特征。这d维特征称为主元，是重新构造出来的d维特征，而不是简单地从n维特征中去除其余n-d维特征。

2 引言

内积与投影

导读：大部分人估计都知道PCA是将数据点向新的方差最大的单位向量做投影，但是什么是到向量的投影，它和内积又有什么关系呢？

两个向量的内积定义为：

(a1,a2,⋯,an)T⋅(b1,b2,⋯,bn)T=a1b1+a2b2+⋯+anbn

两个向量的内积是一个实数，从数学计算的角度很容易计算出两个向量的内积，但是怎么从几何意义上理解内积呢？
为了简单，看二维向量的内积。假设有两个向量 A=(x1,y1)T 和 B=(x2,y2)T ，从 A 点对 OB 向量做一条垂线，垂线与 OB 的交点C即为 A 在 B 上的投影，再设 A 与 B 的夹角是 α，那么点 A 到 B 的投影长度即为 |A|cos(α)。

图1：内积和投影

向量的内积还有另外一种表达方式：

A⋅B=|A||B|cos(α)

即 A 与 B 的内积等于 A 到 B 的投影长度乘以 B 的模。再进一步，如果我们假设 B 为单位向量，其模为1，即|B|=1，那么就变成了：

A⋅B=|A|cos(α)

也就是说，假设向量 B 为单位向量（模为1），则A与单位向量B的内积值等于A向单位向量B所在直线的投影长度，即图 1 中 OC 的长度。

基与基的坐标

如图 1 所示，众所周知，向量 B 用点坐标表示为（3,2）T。其中的 3 本质上是向量在 x 轴上的投影对应的值是3， 2是向量在 y 轴上的投影对应的值是2。也就是说我们其实隐式引入了一个定义：以x轴和y轴上正方向长度为1的向量为标准。更正式的说，向量(x,y)T实际上表示向量在x轴和y轴的线性组合：

x(1,0)T+y(0,1)T

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

图2：什么是基

所以，要准确描述向量，首先要确定一组基，然后给出在基所在的各个直线上的投影值，这些投影值就是当前基所确定的坐标。只不过我们经常省略第一步，而默认以(1,0T)和(0,1)T为基。我们之所以默认选择(1,0T)和(0,1)T为基，因为它们分别是x和y轴正方向上的单位向量，因此就使得二维平面上点坐标和向量一一对应，非常方便。但实际上任何两个线性无关的二维向量都可以成为一组基，所谓线性无关在二维平面内可以直观认为是两个不在一条直线上的向量。
例如，(1,1)T和(−1,1)T也可以成为一组基。一般来说，我们希望基的模是1，因为从上节中内积的意义可以看到，如果基的模是1，那么就可以方便的用向量点乘基（坐标与基的内积）而直接获得其在基上的坐标了！实际上，对应任何一个向量我们总可以找到其同方向上模为1的向量，只要让两个分量分别除以模就好了。例如，上面的基(1,1)T和(−1,1)T可以变为(12√,12√)T和(−12√,12√)T。

现在，我们想获得(3,2)T在新基上的坐标，即在两个方向上的投影矢量值，那么根据内积的几何意义，我们只要分别计算(3,2)T和两个基的内积，不难得到新的坐标为(52√,−12√)T。图3给出了新的基以及(3,2)T在新基上坐标值的示意图（建议左上斜四十五度看）：

图3：基变换

另外这里要注意的是，我们列举的例子中基是正交的（即内积为0，或直观说相互垂直），但可以成为一组基的唯一要求就是线性无关，非正交的基也是可以的。不过因为正交基有较好的性质，所以一般使用的基都是正交的。（如果新基是正交矩阵，这里即从旧基到新基的过渡矩阵为正交矩阵C，那么非常方便的可以求得在新基下的坐标 η=C−1ξ=CTξ，正交方阵是欧氏空间中标准正交基到标准正交基的过渡矩阵。）

基变换的矩阵表示

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

(1/2√1/2√−1/2√1/2√)T(32)=(5/2√−1/2√)

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

(1/2√1/2√−1/2√1/2√)T(112233)=(2/2√04/2√06/2√0)(1)

于是一组向量的基变换被干净的表示为矩阵的相乘。

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

数学表示为：

⎛⎝⎜⎜⎜⎜⎜pT1pT2⋮pTR⎞⎠⎟⎟⎟⎟⎟(a1a2⋯am)=⎛⎝⎜⎜⎜⎜⎜pT1a1pT2a1⋮pTRa1pT1a2pT2a2⋮pTRa2⋯⋯⋱⋯pT1ampT2am⋮pTRam⎞⎠⎟⎟⎟⎟⎟

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

特别要注意的是，这里R可以小于N，而R决定了变换后数据的维数。也就是说，我们可以将一N维数据变换到更低维度的空间中去，变换后的维度取决于基的数量。因此这种矩阵相乘的表示也可以表示降维变换，而PCA的降维也是这个过程，正如式（1）所示，在新基中原来二维的三个点可以之用一维的三个点表示，因为在新基中第二维的坐标都为0。

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

这里只是简单的介绍基和基的变换，想要更深入理解基相关的内容建议学习线性代数和矩阵论相关的知识，包括线性空间、线性无关组、正交矩阵、线性空间的基与坐标，基变换与坐标变换，过渡矩阵等内容。矩阵论的书籍推荐：程云鹏：《矩阵论》

期望方差协方差

给定一个含有m个样本的集合X={X1,…,Xm}，
均值：

X¯=∑mi=1Xim

标准差：

s=∑mi=1(Xi−X¯)2m−1−−−−−−−−−−√

方差：

s2=∑mi=1(Xi−X¯)2m−1

很显然，均值描述的是样本集合的平均值，而标准差给我们描述的则是样本集合的各个样本点到均值的距离之平均。比如两个考生四次测试的成绩分别为[80，100，80，60] 和 [78，80，82，80]，两个考生的平均成绩均是80，但是我们愿意选择哪个考生呢？也就是选择哪个发挥更加稳定的考生呢？即每次考试偏离平均值的程度，那就用标准差来描述。计算两者的标准差，前者是14.1，后者是1.41，显然后者较为集中，故其标准差小一些，标准差描述的就是这种“散布度”。之所以除以m-1而不是除以m，是因为这样能使我们以较小的样本集更好的逼近总体的标准差，即统计上所谓的“无偏估计”。而方差则仅仅是标准差的平方。

协方差
方差是描述一维数据样本本身相对于均值的偏离程度，比如上面所说的考生某门考试成绩的稳定程度的描述。但是，如果遇到含有多维数据的数据集，比如要统计多个学科的考试成绩，面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，比如，某个考生数学成绩很好（数学成绩平均值高，且其方差小），他的物理成绩也很好（物理成绩平均值高，且其方差小），这些可以用均值和方差来描述，但是我们想要了解考生的数学成绩很好是不是和他的物理成绩很好相关呢？或者他的数学成绩很好和他的历史成绩差负相关呢？（一般来说，数学成绩好，逻辑思维强，对应的物理成绩也会相对好些，个人觉得哈）。那么要描述这个考生数学成绩和物理成绩是否有正相关，或者描述其数学成绩和历史成绩是否负相关，那么需要协方差来度量两个随机变量X,Y关系的统计量，我们可以仿照方差的定义：

var(X)=∑mi=1(Xi−X¯)(Xi−X¯)m−1

随机变量X,Y的协方差的定义为

cov(X,Y)=∑mi=1(Xi−X¯)(Yi−Y¯)m−1

协方差的结果有什么意义呢？如果结果为正值，则说明两者是正相关的(从协方差可以引出“相关系数”的定义)，也就是说如果考生的数学成绩和物理成绩的协方差为正，那么代表他的数学成绩和物理成绩是相关的，如果数学成绩和历史成绩的协方差为负，那么他的数学成绩和物理成绩是负相关的，如果协方差为0，那么代表考生的数学成绩和物理成绩没有相关性。

从协方差的定义上我们也可以看出一些显而易见的性质，如：
1. cov(X,X)=var(X)
2. cov(X,Y)=cov(Y,X)

协方差矩阵

上面提到的X,Y是一个二维问题，而协方差也只能处理二维问题，那维数多了自然就需要计算多个协方差，分别计算两两之间的协方差，比如n维的数据集就需要计算C2n=n!(n−2)!∗2个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

Cn×n=⎛⎝⎜⎜⎜⎜⎜c11c21⋮cn1c12c22⋮cn2⋯⋯⋱⋯c1nc2n⋮cnn⎞⎠⎟⎟⎟⎟⎟其中ci,j=cov(Dimi,Dimj))

这个定义还是很容易理解的，我们可以举一个简单的三维的例子，假设数据集有{X,Y,Z}三个维度，则协方差矩阵为：

C3×3=⎛⎝⎜cov(X,X)cov(Y,X)cov(Z,X)cov(X,Y)cov(Y,Y)cov(Z,Y)cov(X,Z)cov(Y,Z)cov(Z,Z)⎞⎠⎟

可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。C可以作为刻画不同维度（分量）之间相关性的一个评判量。若不同分量之间的相关性越小，则C的非对角线元素的值就越小，特别地，若不同分量彼此不相关，那么C就变成了一个对角阵。

特别提醒的是，协方差矩阵是针对数据的维度定义的，是代表维度两两之间的相关性，换句话说也是特征两两之间的相关性，切记！

更有意思的是，两个随机变量(X,Y)，根据协方差的定义：

cov(X,Y)=∑mi=1(Xi−X¯)(Yi−Y¯)m−1=(X−X¯)⋅(Y−Y¯)m−1

可以看到，在均值为0的情况下，两个维度向量的协方差简洁的表示为其内积除以元素数m−1。
那么推广到求协方差矩阵，假设有一个样本数据集X，X用一个n∗m的矩阵表示，其中n表示每个样本的特征维度，m表示样本点个数，那我们的协方差矩阵应该是一个n * n的对称矩阵，从上述协方差和内积的关系，不难得出协方差矩阵，具体推导过程也可以参考文章关于协方差矩阵的理解：

C=1m−1SST，其中S为X按列去均值的矩阵

比如X为三维特征（a,b,c），m 个样本所组成的矩阵，其中心化后为 S 设为：

S=⎛⎝⎜a1b1c1a2b2c2⋯⋯⋯ambmcm⎞⎠⎟ST=⎛⎝⎜⎜⎜a1a2⋯amb1b2⋯bmc1c2⋯cm⎞⎠⎟⎟⎟

那么协方差矩阵

C=1m−1⎛⎝⎜⎜∑mi=1a2i∑mi=1(biai)∑mi=1(ciai)∑mi=1(aibi)∑mi=1b2i∑mi=1(cibi)∑mi=1(aici)∑mi=1(bici)∑mi=1c2i⎞⎠⎟⎟=1m−1SST

下面我们举例来求协方差矩阵，包括两种求解方式（其实本质上是一样的），一种是按照定义求，另一种直接用矩阵的乘积去求，最后和Python中的Numpy内部的cov函数进行比较计算是否正确。

#第一种方式
#根据定义$cov(X,Y)=\frac{\sum_{i=1}^m (X_{i}-\bar{X})(Y_{i}-\bar{Y})}{m-1}$求得协方差矩阵的每个元素
#注意其求和的过程便是求内积的过程，从整个矩阵的角度来看，也是矩阵相乘的过程
def cov_define(dataSet):
    num_dim = dataSet.shape[0]
    num_samples = dataSet.shape[1]
    cov_mat = np.zeros((num_dim, num_dim))
    mean_data = np.mean(dataSet,1)
    move_mean_data = (dataSet.transpose() - mean_data).transpose()
    for i in range(0, num_dim):
        for j in range(0, num_dim):
            cov_mat[i,j] = np.sum(move_mean_data[i,:] * move_mean_data[j,:]) / (num_samples - 1)
    return  cov_mat
1
2
3
4
5
6
7
8
9
10
11
12
13
14
1
2
3
4
5
6
7
8
9
10
11
12
13
14

#第二种方式
#cov = 1/(n-1) * X * X^T，其中X为m*n的矩阵,m为样本数，n为维度
def cov_tran(dataSet):
    mean_data = np.mean(dataSet,1)
    move_mean_data = (dataSet.transpose() - mean_data).transpose()
    my_cov = (np.dot(move_mean_data,move_mean_data.transpose())) / (dataSet.shape[1] - 1)
    return my_cov1
2
3
4
5
6
7
1
2
3
4
5
6
7

#主程序模块
#get cov by cov = 1/(n-1) * X * X^T
sampleDataSet = np.random.randint(0,10,(3,5)) #X为5个3维的样本点
covdata_tran = cov_tran(sampleDataSet)

#get cov by the function cov of python.numpy
mean_data = np.mean(sampleDataSet,1)
move_mean_sample = (sampleDataSet.transpose() - mean_data).transpose()
covdata_np = np.cov(move_mean_sample,rowvar=1)#If rowvar is non-zero (default), then each row represents a variable, with observations in the columns.
covdata_def = cov_define(sampleDataSet)

print 'sampleDataSet = \n',sampleDataSet
print 'covdata_tran = \n',covdata_tran
print 'covdata_np = \n',covdata_np
print 'covdata_def = \n',covdata_def1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15

#输出结果如下：
sampleDataSet = 
[[6 3 6 6 2]
 [2 0 5 8 2]
 [7 3 3 8 0]]
covdata_tran = 
[[  3.8    3.95   5.1 ]
 [  3.95   9.8    5.4 ]
 [  5.1    5.4   10.7 ]]
covdata_np = 
[[  3.8    3.95   5.1 ]
 [  3.95   9.8    5.4 ]
 [  5.1    5.4   10.7 ]]
covdata_def = 
[[  3.8    3.95   5.1 ]
 [  3.95   9.8    5.4 ]
 [  5.1    5.4   10.7 ]]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

更形象的协方差矩阵的实例可以参考协方差矩阵的实例与意义和协方差矩阵的几何解释。

3 原理与流程

前面长篇大论，也该言归正传了，但是各位一定要相信，前面所有的都是有帮助的，理解了前面的预备知识，对于PCA的原理的理解就可以说是水到渠成了。

维基百科中主成分分析（Principal Component Analysis(PCA)）的数学定义是：一个正交化线性变换，把数据变换到一个新的坐标系统中，使得这一数据的任何投影的第一大方差在第一个坐标（称为第一主成分）上，第二大方差在第二个坐标（第二主成分）上，依次类推。

PCA 是最常用的线性降维方法，它的目标是通过某种线性投影，将高维的数据映射到低维的空间中表示，并期望在所投影的维度上数据的方差最大，以此使用较少的数据维度，同时保留住较多的原数据点的特性。通俗的理解，如果把所有的点都映射到一起，那么几乎所有的信息（如点和点之间的距离关系）都丢失了，而如果映射后方差尽可能的大，那么数据点则会分散开来，以此来保留更多的信息。

下面我们来看个例子：

X=(1113234424)

中心化之后的矩阵为：

S=(−1−2−10002101)

其在平面上的坐标如图所示：

图4：PCA的样本点

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

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

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

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

PCA数学表达

下面，我们用数学方法表述这个问题。
假设样本集 X={x1,x2,...,xm}，样本点xi在新空间中超平面上的投影是WTxi，若所有样本点的投影能尽可能分开，则应该使投影后样本点的方差最大化，投影后样本点的方差可以表示为：∑iWTxixTiW，于是优化目标可写为：

maxW tr(WTXXTW) s.t. WTW=I(2)

对于式（2）使用拉格朗日乘子法可得

XXTW=λW

于是，只需对协方差矩阵XXT进行特征值分解，将求得的特征值排序：λ1≥λ2≥...≥λn，再取前d个特征值对应的特征向量构成W=(w1,w2,...,wd)，就是主成分分析的解。
此外，对PCA 优化问题比较详细的推导的可以参见博文PCA－Principal Components Analysis数学原理和主成分分析-最大方差解释。

PCA 算法步骤

输入：样本集 X={x1,x2,...,xm}；

过程：
1. 对所有样本进行中心化：xi←xi−1m∑mixi;
2. 计算样本的协方差矩阵XXT
3. 对协方差矩阵XXT做特征分解；
4. 取最大的d个特征值所对应的特征向量w1,w2,...,wd。

输出：投影矩阵W=(w1,w2,...,wd)

4 算法实例

上文提到的，有5个二维的样本集X={(11),(13),(23),(44),(24)}:

X=(1113234424)

PCA过程：
1. 对所有样本进行中心化，得到矩阵S

S=(−1−2−10002101)

计算协方差矩阵C=1m−1S∗ST

C=14(−1−2−10002101)⎛⎝⎜⎜⎜⎜⎜⎜−1−1020−20011⎞⎠⎟⎟⎟⎟⎟⎟=(1.5111.5)
矩阵分解，求特征值和特征向量，得到：

λ1=0.5, λ2=2.5

对应的特征向量：

w1=(−1/2√1/2√), w2=(1/2√1/2√)

选取λ2=2.5对应的特征向量w2=(1/2√ 1/2√)
得到投影矩阵为W=(w2)，即可以将样本点X投影到W平面，其投影是Y，

Y=WTX=(1/2√1/2√)(−1−2−10002101)=(−3/2√−1/2√03/2√−1/2√)

输出结果

sampleDataSet =
[[1 1 2 4 2]
 [1 3 3 4 4]]
move_mean_sample =
[[-1. -1.  0.  2.  0.]
 [-2.  0.  0.  1.  1.]]
np_cov =
[[ 1.5  1. ]
 [ 1.   1.5]]
eigV=
[ 0.5  2.5]
eigVector=
[[-0.70710678  0.70710678]
 [ 0.70710678  0.70710678]]
pca_data=
[[-2.12132034 -0.70710678  0.          2.12132034  0.70710678]]1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16