PCA主成分分析原理简述

原贴: http://www.fengchang.cc/post/89

之前读书看过,没有总结,忘得差不多了,最近读书又看了一遍,写总结如下。

 

还是老方法,先建立直觉,提供一个比较好的建立直觉的资料:A TUTORIAL ON PRINCIPAL COMPONENT ANALYSIS 我的圈评

 

我就引用这里面的例子来说明。

PCA主成分分析原理简述

上图是一个弹簧实验的示意图。假设我们通过三个摄影机 A B C来记录连续若干时刻的球的位置,目标是求出球的运动方程。首先弹簧本身会沿着x方向伸缩,但是由于现实世界的误差的原因,在y方向上依然可能会有轻微的扰动。在已经确定好x, y, z的坐标轴的情况下,这个问题倒也不难,难就难在,在现实世界的很多问题中,并不存在这样的设定好的坐标轴给你参考,所以基本上你收集到的信息可能是比较杂乱的,比如上图中,我们假设每个摄像机记录两个维度的信息A记录(Ax, Ay), B记录(Bx, By), C记录(Cx, Cy), 注意这里x, y 只表示两个维度,不代表在哪个坐标轴,那么结果可以认为就是一个摄像机记录一个观察角度下球的向量,最终的一条记录数据长这样(Ax,Ay,Bx,By,Cx,Cy),这种数据我们可以想象的一个缺陷是,它有6个维度,但是对实际弹簧球的运动来说,实在太多了,且没有重点,因为弹簧只在一个方向上运动为主,那么如果我们拿这样的数据集去做机器学习中的训练,这个训练数据其实质量是相当差的。

 

又好比我收集学生的如下6个维度的信息:年龄,性别,智商,肤色,是否直发,眼睛颜色,来预测其人物理考试的成绩会落在那个分数段,可想而知,智商对结果的影响一定是最大的,而是否直发这个维度的信息可能根本就无关紧要。那么PCA要做的一件事就是通过“降维”来处理训练数据。

 

问题来了,怎么降呢?如果是上述物理考试成绩预测的问题,直接把是否直发这样的维度删掉,就是一种处理方案,这样就把6维训练数据,降到5维了。但是这么粗暴的方法,并不适合所有的问题,比如上述弹簧球的问题就不适合,因为你并不能明确知道到底哪个维度的数据是无关紧要的,甚至实际上可能某一个维度的信息,实际上是由另外两个维度的信息叠加合成的(可以参考牛顿物理学中基本的运动的合成),这个时候问题就比较复杂了。

 

那么PCA是怎么做的呢?简单来说就是线性代数当中的基变换。

 

姿势前提:基:假设线性子空间的基B={v1,v2,...,vk}, 向量 a = v1c1+v2c2+...+vkck,那么(c1,c2,...,ck)即为a基于B的坐标。

 

这里不得不先说一下,理解PCA的前提是基本的线性代数知识,涉及到方阵,矩阵对角化,特征值,特征向量,这几个概念理解得越好,对PCA的理解就越水到渠成。我本人只是在大学数学学过线性代数,没有多精通,所以也是先复习了一下这几个概念,所以自认为对PCA的理解还没有达到多深的水平。

 

假设我们的原始训练数据是一个m乘n矩阵,其中m是观察的原始维度,n是训练样本的个数, 所以相当于每一列是一个观察样本,每一行是一个维度的所有观察集合,那么针对上面弹簧实验的问题,假如我们在记录了200个数据时刻的数据,那么训练数据就是6*200的矩阵, 我就记录为X吧,这样的矩阵,如果左乘一个矩阵P,例如P的尺寸为5*6, 那么结果就变成一个5*200的矩阵,先不管目的为何,这样左乘,实际上就是对原始数据做了一个基变换,且从6维降到了5维。

 

然后我们希望达到一个目标,如果在基变换后的训练数据,各个维度之间的相关度都很低,那就相当的好,为什么呢?因为在你胡乱观察数据的时候,收集来的信息可能是有很大冗余的,比如你既收集了一个人的出生年月,又收集了他的年龄,又收集了他的属相,然后存为3个维度的信息,这很明显就有冗余,因为这三者是相关的,从出生年月这个维度,完全可以推导出另外两个维度。这种维度之间的冗余性,就是PCA要解决的问题,通过线性变换,达到一个效果,让维度之间尽可能独立,减少冗余性。至于为什么线性变换可以达到这个效果,不在本文的讨论范围内了, 我也理解得不深。但这就是我们的优化目标,那么下面就需要用一种正式的数学公式来衡量这个冗余性,PCA是这样做的:

PCA主成分分析原理简述

假设我们的训练数据叫X(就以我们上面说的6*200的训练数据为例),那么通过上图中的公式定义了一个叫covariance matrix的矩阵,这个矩阵里面的每一个元素,都代表着两个维度之间的协方差,协方差这个东西,就是定义两个变量相关性的的一个量化指标,这个指标的定义跟什么均值,方差之类的指标是一个level的,我这里就不再赘述了。总之上图就可以定义维度与维度之间的协方差(如果大家对比协方差的标准公式定义可能会有一点困惑,就是均值哪去了?其实是这样的,这里的数据都做了预处理,每个原始数据都已经减去了均值,这样处理后的均值都变成0了,所以看起来均值就不见了)。

 

接下来这个协方差矩阵就是一个有代表性的东西,它的每个元素都代表了两个维度之间的相关性,假设第i行第j列的元素a<i,j>,就代表了第i维数据跟第j维数据的相关性。那么还记得我们的优化目标吗?就是要通过线性变换,让不同维度之间的相关性尽可能小,对应到这个协方差矩阵,就是希望该矩阵对角线以外的相关性都是0,只有对角线上的元素非0。这里就要用到线性代数了,首先根据定义,这个协方差矩阵一定是个方阵,方阵的话,如果能对角化就说明,我们总是有办法找到这样一个线性变换,让它变成对角矩阵,这不就达到我们的优化目标了吗?

插入一个矩阵可对角化的定义:如果存在一个可逆矩阵 P 使得 P −1AP 是对角矩阵,则它就被称为可对角化

 

再梳理一下,可以这样描述PCA做的事,对原始训练数据矩阵(记为X),要寻找一个线性变换矩阵P,使得PX对应的协方差矩阵为对角阵。找到这样一个P之后,PX就是我们要变换后的结果。

 

这应该是能一句话说明的最简单描述了,但是这里P的维度其实是要注意的,因为理论上可能有多种,比如对于上面弹簧的训练数据(6*200的训练数据),P可能是5*6,也可能是4*6或者3*6的尺寸,总之就是k*6,这个k就是我们降维后想要达到的维度,通常是人为指定的,至于如何选这个k,就是另一个话题了,以后有空再说。

 

下面说一下在指定了k之后,如何解这个P,基本上就是所谓矩阵对角化那一套了,我这里不说底层细节,因为其实展开了说挺难的,涉及到的都是线性代数的运算知识,包括特征值,大体就是把矩阵的特征值全找出来,然后挑出k个最大的特征值,找出对应的k个特征向量,这k个特征向量拼起来就是矩阵P。

 

下面给一个python sklearn的代码实现,用的是Iris数据集,长这样,也不多,就150条记录。

PCA主成分分析原理简述

比较简单,而且没有底层的实现细节,因为sklearn已经封装好了,不过通过debug的过程,可以看看矩阵的维度是怎么变化的,加深对PCA该过程的理解。

 

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

#!/usr/bin/python

# -*- coding: utf-8 -*-

 

"""

=========================================================

PCA example with Iris Data-set

=========================================================

 

Principal Component Analysis applied to the Iris dataset.

 

See `here <https://en.wikipedia.org/wiki/Iris_flower_data_set>`_ for more

information on this dataset.

 

"""

print(__doc__)

 

 

# Code source: Gaël Varoquaux

# License: BSD 3 clause

 

import numpy as np

import matplotlib.pyplot as plt

from mpl_toolkits.mplot3d import Axes3D

 

 

from sklearn import decomposition

from sklearn import datasets

 

np.random.seed(5)

 

centers = [[11], [-1-1], [1-1]]

iris = datasets.load_iris()

= iris.data

= iris.target

 

fig = plt.figure(1, figsize=(43))

#clear figure

plt.clf()

ax = Axes3D(fig, rect=[00, .951], elev=48, azim=134)

 

# clear axis

plt.cla()

pca = decomposition.PCA(n_components=3)

pca.fit(X)

= pca.transform(X)

 

for name, label in [('Setosa'0), ('Versicolour'1), ('Virginica'2)]:

    ax.text3D(X[y == label, 0].mean(),

              X[y == label, 1].mean() + 1.5,

              X[y == label, 2].mean(), name,

              horizontalalignment='center',

              bbox=dict(alpha=.5, edgecolor='w', facecolor='w'))

# Reorder the labels to have colors matching the cluster results

= np.choose(y, [120]).astype(np.float)

ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y, cmap=plt.cm.nipy_spectral,

           edgecolor='k')

 

ax.w_xaxis.set_ticklabels([])

ax.w_yaxis.set_ticklabels([])

ax.w_zaxis.set_ticklabels([])

 

plt.show()

 

刚才说了,对线性代数以上若干概念的理解直接关乎对PCA理解的深度,附知乎帖子一篇: 如何理解矩阵特征值

 

代码原贴:PCA example