主成成分分析pca算法 原理解析

目前，pca算法已经广泛应用于各方面，就拿图像处理，经常做的一件事就是当提取的图像特征维度比较高时，为了简化计算量以及储存空间，需要对这些高维数据进行一定程度上的降维，并尽量保证数据的不失真。

 

先举个例子，方便理解：

    1）对于一个训练集，100个sample(i=1,2,3,...,100)，特征X_i是20维.[X_i1,X_i2,X_i3,...X_ij,...,X_i20](j=1,2,..,20)，那么它可以建立一个20*100的样本矩阵M。

    2）紧接着我们开始求这个样本的协方差矩阵，得到一个20*20的协方差矩阵，计算过程如下：

            •先求解出Xi的平均Xav=(∑x_i)/20;

            •对每一个X_i,计算Xi-X_av,即Mi(第 i 行)变为 M_i-X_av,记为Mn;

            •则容易得到协方差矩阵Z为Mn*Mn'( ' 表示转置 ) 。

    3）然后求出这个协方差矩阵Z_20x20的特征值和特征向量，一般情况下应该有20个特征值和特征向量，现在根据特征值的大小，取出较大的特征值以及其所对应的特征向量，（假设提取的特征值为较大的5个特征值），那么这5个特征向量就会构成一个20*5的矩阵V，这个矩阵就是我们要求的特征矩阵。

    4）用Mn'去乘以V，得到一个base矩阵（*），大小为100x5。

    5）任取一个样本1x100,乘上这个100*5的特征矩阵，就得到了一个1*5的新的样本，显然每个sample的维数下降了，然后再用这个1x5向量去比较相似性。

 

注：

›上述3）过程中特征值的选取在不确定具体要降到多少维的情况下，一般还可以根据n个特征值之和大于总和的90%进行选取。

›上面的（*）处base矩阵的求解不唯一，也可以自行修正。

 

大致说完了PCA降维的过程，现在会有人要问为什么只提取特征值较大的几个特征值就可以近似代替原样本矩阵呢。

好了，话不多说。下面就讲讲矩阵的特征值和特征向量的数学意义：

 

为简单起见，以二维矩阵A=[1 0;0 -1]（矩阵的秩为2）为例，以平面上一点（x,y）经过A变换后变为（x',y'）若这两点在一条直线在，那么可以理解为矩阵A的作用恰好使得向量[x y]' 只是在原有方向上变换了长度而已，即Ax=λx (x为一列向量).对于A矩阵，容易得到A的两个特征值及相应的特征向量 λ1=1 ,e1=[1 0]' , λ2=-1 ,e2=[0 -1]' ,二维平面上任意一点（x,y）=b1*(1,0)+b2*(0,-1)（b1,b2均为实常数）; 那么A[x y]'=A*(b1*e1+b2*e2)=b1*λ1+b2*λ2 =∑b_iλ_i ;

 

把这个公式推广到高维空间，在计算（x',y'）对于λ值比较小的特征维可以忽略.

B=[1 0;0 0.01] ,其中B的两个特征值及相应的特征向量 λ1=1 ,e1=[1 0]' , λ2=0.01 ,e2=[0 1]'

那么x=[2 3]' 经过B变换为 Bx=[2 0.03]';

如果我们认为λ2远小于λ1，忽略掉λ2时，B_n=[1 0;0 0],B_nx=[2 0]'≈[2 0.03].

 

通俗点讲，pca算法就是去寻找那些在该维度上方差比较大的维，同时忽略比较平均的维度。假如上面所说的X特征向量的第一个元素都为1，那么这一列数据是可以忽略的，因为他并不能起到区分的作用，相反我们是要寻找那些在某一维度分布比较广的维，并提取出来。

打个比方，平面区域一个斜75度的椭圆，且长轴远大于短轴，那么椭圆上的点在短轴上的分布明显弱于长轴，当短轴远小于长轴时，近似为一条直线，便失去了短轴这个维度。

下面贴上一个matlab版的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

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

%一个修改后的PCA进行人脸识别的Matlab代码

clear;

% calc xmean,sigma and its eigen decomposition

allsamples=[];%所有训练图像

for i=1:40

    for j=1:5

      a=imread(strcat('C:\Documents and Settings\Foreigners\桌面\ORL\s',num2str(i),'\',num2str(j),'.bmp'));

      % imshow(a);

      b=a(1:112*92); % b是行矢量 1×N，其中N＝10304，提取顺序是先列后行，即从上到下，从左到右

      b=double(b);

      allsamples=[allsamples; b];  % allsamples 是一个M * N 矩阵，allsamples 中每一行数据代表一张图片，其中M＝200

  end

end

samplemean=mean(allsamples); % 平均图片，1 × N

for i=1:200 xmean(i,:)=allsamples(i,:)-samplemean; % xmean是一个M × N矩阵，xmean每一行保存的数据是“每个图片数据-平均图片”

end;

 

sigma=xmean*xmean';   % M * M 阶矩阵

[v d]=eig(sigma);

d1=diag(d);

[d2 index]=sort(d1); %以升序排序

cols=size(v,2);% 特征向量矩阵的列数

for i=1:cols

    vsort(:,i) = v(:, index(cols-i+1) ); % vsort 是一个M*col(注:col一般等于M)阶矩阵，保存的是按降序排列的特征向量,每一列构成一个特征向量

    dsort(i)   = d1( index(cols-i+1) );  % dsort 保存的是按降序排列的特征值，是一维行向量

end  %完成降序排列

%以下选择90%的能量

dsum = sum(dsort);

    dsum_extract = 0;

    p = 0;

    while( dsum_extract/dsum < 0.90)

        p = p + 1;

        dsum_extract = sum(dsort(1:p));

    end

i=1;

% (训练阶段)计算特征脸形成的坐标系

while (i<=p && dsort(i)>0)

    base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i);   % base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化，详见《基于PCA的人脸识别算法研究》p31

    i = i + 1;

end

 

% add by wolfsky 就是下面两行代码，将训练样本对坐标系上进行投影,得到一个 M*p 阶矩阵allcoor

allcoor = allsamples * base;

accu = 0;

 

% 测试过程

for i=1:40

    for j=6:10 %读入40 x 5 副测试图像

        a=imread(strcat('C:\Documents and Settings\Foreigners\桌面\ORL\s',num2str(i),'\',num2str(j),'.bmp'));

        b=a(1:10304);

        b=double(b);

        tcoor= b * base; %计算坐标，是1×p阶矩阵

        for k=1:200

                mdist(k)=norm(tcoor-allcoor(k,:));

        end;

        %三阶近邻

         

        [dist,index2]=sort(mdist);

        class1=floor( index2(1)/5 )+1;

        class2=floor(index2(2)/5)+1;

        class3=floor(index2(3)/5)+1;

        %class=class1;%%blue_lg

        if class1~=class2 && class2~=class3

            class="class1";

        elseif class1==class2

            class="class1";

            %elseif class2==class3

          % class="class2";

        end;

        if class==i

            accu=accu+1;

        end;

    end;

end;

accuracy=accu/200 %输出识别率

%zuobiao=[1:100];

%plot(zuobiao,accuracy);

　　

转载： http://www.cnblogs.com/blue-lg/archive/2012/05/14/2499581.html

文章由两部分构成，第一部分主要讲解PCA算法的步骤，第二部分讲解PCA算法的原理。

 

那么首先进入第一部分

                                   --PCA算法的步骤

---

---

① 样本矩阵X的构成

    假设待观察变量有M个，其实相当于一个数据在M维各维度上的坐标，我们的目标是在保证比较数据之间相似性不失真的前提下，将描述数据的维度尽量减小至L维（L<M）。

    样本矩阵X在这里用x₁,x₂,...,x_N共N个数据（这些数据都是以列向量的形式出现）来表示，那么X=[x₁ x₂ ... x_N]_MxN显而易见。

 

② 计算样本X均值

    计算第m维（m=1,2,...,M）的均值如下：

        

             

 

③ 计算观察值与均值的偏差

    在每一维上，用当前值X[m,n]减去u[m]，用矩阵运算表示如下：

    

       

    明显，h是一行向量，u是一列向量。

 

④ 计算协方差矩阵

    根据协方差矩阵的计算定义wikihttp://en.wikipedia.org/wiki/Covariance_matrix

    我们认为b_i代表B的第i行，那么由协方差矩阵

                  

 

    推知 

    <>表示向量的内积，也就是每一项的积之和。

 

⑤ 计算协方差矩阵C的特征值和特征向量

    若XA=aA,其中A为一列向量，a为一实数。那么a就被称为矩阵X的特征值，而A则是特征值a对应的特征向量。    

    顺便扯一下，在matlab里面，求解上述A以及a，使用eig函数。如[D,V] = eig(C)，那么D就是n个列特征向量，而V是对角矩阵，对角线上的元素就是特征值。

 

⑥ 排序

    将特征值按从大到小的顺序排列，并根据特征值调整特征向量的排布。

    D'=sort(D);V'=sort(V);

 

⑦ 计算总能量并选取其中的较大值

    若V'为C的对角阵，那么总能量为对角线所有特征值之和S。

    由于在步骤⑥里面已对V进行了重新排序，所以当v'前几个特征值之和大于等于S的90%时，可以认为这几个特征值可以用来"表征"当前矩阵。

    假设这样的特征值有L个。

 

⑧ 计算基向量矩阵W

    实际上，W是V矩阵的前L列，所以W的大小就是 MxL 。

 

⑨ 计算z-分数（这一步可选可不选）

     Z[i,j]=B[i,j]/sqrt(D[i,i])

 

⑩ 计算降维后的新样本矩阵

 

                

     W^*表示W的转置的共轭矩阵，大小为 LxM , 而Z的大小为 MxN , 所以Y的大小为 LxN , 即降维为 N 个 L 维向量。

     

现在进入第二部分

                                   --PCA算法的原理

---

---

 首先让我们来了解一下特征值和特征向量的几何意义

 

根据上面步骤⑤的说法，很明显的发现一个特定的变换，它的特征向量是这样一种向量，它经过这种特定的变换后保持方向不变，只是进行长度上的伸缩而已(再想想特征向量的原始定义Ax=cx, cx是方阵A对向量x进行变换后的结果，但显然cx和x的方向相同)。

一个变换（矩阵）可由它的所有特征向量完全表示，而每一个向量所对应的特征值，就代表了矩阵在这一向量上的贡献率——说的通俗一点就是能量（power）。

请看，当x=c₁·x₁+c₂·x₂时，而x₁和x₂均为矩阵A的特征向量(λ1和λ2分别为他们的特征值)，那么Ax=c1·λ1·x1+c2·λ2·x2 ,只不过是在两基方向上的投影长度发生了变化。

 

我们知道，一个变换可由一个矩阵乘法表示，那么一个空间坐标系也可视作一个矩阵，而这个坐标系就可由这个矩阵的所有特征向量表示，用图来表示的话，可以想象就是一个空间张开的各个坐标角度，这一组向量可以完全表示一个矩阵表示的空间的“特征”，而他们的特征值就表示了各个角度上的能量（可以想象成从各个角度上伸出的长短，越长的轴就越可以代表这个空间，它的“特征”就越强，或者说显性，而短轴自然就成了隐性特征），因此，通过特征向量/值可以完全描述某一几何空间这一特点，使得特征向量与特征值在几何（特别是空间几何）及其应用中得以发挥。

                                       

假设2x2矩阵的两个特征向量如上图红色虚线和绿色线所示，那么当椭圆的长轴远大于短轴时，可以将2维的点降维成一维，及当前向量在长轴上的投影值。

 

大家可以看看这篇blog，写的蛮不错的。http://blog.csdn.net/hexbina/article/details/7525850

 

再举一个例子，如果大家再引入一个M维的向量p和已经降维过后的向量q比较的话，需要先减去②里面每一维的均值，然后除以每一维的特征值的平方根，再乘以W^*就OK了。

用公式表示为:

                                     

 

转载：http://www.cnblogs.com/blue-lg/archive/2012/07/23/2604995.html

1.问题描述  

在许多领域的研究与应用中，往往需要对反映事物的多个变量进行大量的观测，收集大量数据以便进行分析寻找规律。多变量大样本无疑会为研究和应用提供了丰富的信息，但也在一定程度上增加了数据采集的工作量，更重要的是在大多数情况下，许多变量之间可能存在相关性，从而增加了问题分析的复杂性，同时对分析带来不便。如果分别对每个指标进行分析，分析往往是孤立的，而不是综合的。盲目减少指标会损失很多信息，容易产生错误的结论。

 

2.过程

主成分分析法是一种数据转换的技术，当我们对一个物体进行衡量时，我们将其特征用向量（a1,a2,a3,...an）进行表示，每一维都有其对应的variance（表示在其均值附近离散的程度）；其所有维的variance之和，我们叫做总的variance；我们对物体进行衡量时，往往其特征值之间是correlated的，比如我们测量飞行员时，有两个指标一个是飞行技术（x1）,另一个是对飞行的喜好程度（x2），这两者之间是有关联的，即correlated的。我们进行PCA（主成分分析时），我们并没有改变维数，但是我们却做了如下变换，设新的特征为（x1,x2,x3...,xn）;

其中

1)x1的variance占总的variance比重最大；

2)除去x1,x2的variance占剩下的variance比重最大；

....

依次类推；

最后，我们转换之后得到的(x1,x2,...xn)之间都是incorrelated，我们做PCA时，仅取（x1，x2,....xk）,来表示我们测量的物体，其中，k要小于n。主成分的贡献率就是某主成分的方差在全部方差中的比值。这个值越大，表明该主成分综合X1，X2，…，XP信息的能力越强。如果前k个主成分的贡献率达到85%，表明取前k个主成分基本包含了全部测量指标所具有的信息，这样既减少了变量的个数又方便于对实际问题的分析和研究。

    注意，当（a1,a2,a3,...an）之间都是incorrelated时，我们就没有做PCA的必要了

 

数据点在上图所示的方向上进行投影后，数据仍然有着很大的variance,但在下图所示的方向上，投影后的数据的variance就很小。

我们所需要做的就是找到这一系列的向量，使得数据在其上的投影有着较大的variance。

 

3.数学描述

  为了能够找到这一系列的向量，我们对数据进行预处理

1) Alcohol 
2) Malic acid 
3) Ash 
4) Alcalinity of ash 
5) Magnesium 
6) Total phenols 
7) Flavanoids 
8) Nonflavanoid phenols 
9) Proanthocyanins 
10)Color intensity 
11)Hue 
12)OD280/OD315 of diluted wines 
13)Proline

样本数为130，在matlab下按照以上步骤，进行PCA，得到的特征值如下：

选取前k个特征值使得前k个主成分的贡献率达到85%，计算得到的结果为k=1，其对应的特征向量为u=

令X=X*u即可得到新的X，其中X原来维数为130×13，进行PCA后的维数为130×1。

今天看论文的时候又看到了协方差矩阵这个破东西，以前看模式分类的时候就特困扰，没想到现在还是搞不清楚，索性开始查协方差矩阵的资料，恶补之后决定马上记录下来，嘿嘿~本文我将用自认为循序渐进的方式谈谈协方差矩阵。

统计学的基本概念

学过概率统计的孩子都知道，统计里最基本的概念就是样本的均值，方差，或者再加个标准差。首先我们给你一个含有n个样本的集合，依次给出这些概念的公式描述，这些高中学过数学的孩子都应该知道吧，一带而过。

均值：
标准差：
方差：

很显然，均值描述的是样本集合的中间点，它告诉我们的信息是很有限的，而标准差给我们描述的则是样本集合的各个样本点到均值的距离之平均。以这两个集合为例，[0，8，12，20]和[8，9，11，12]，两个集合的均值都是10，但显然两个集合差别是很大的，计算两者的标准差，前者是8.3，后者是1.8，显然后者较为集中，故其标准差小一些，标准差描述的就是这种“散布度”。之所以除以n-1而不是除以n，是因为这样能使我们以较小的样本集更好的逼近总体的标准差，即统计上所谓的“无偏估计”。而方差则仅仅是标准差的平方。

为什么需要协方差？

上面几个统计量看似已经描述的差不多了，但我们应该注意到，标准差和方差一般是用来描述一维数据的，但现实生活我们常常遇到含有多维数据的数据集，最简单的大家上学时免不了要统计多个学科的考试成绩。面对这样的数据集，我们当然可以按照每一维独立的计算其方差，但是通常我们还想了解更多，比如，一个男孩子的猥琐程度跟他受女孩子欢迎程度是否存在一些联系啊，嘿嘿~协方差就是这样一种用来度量两个随机变量关系的统计量，我们可以仿照方差的定义：

来度量各个维度偏离其均值的程度，标准差可以这么来定义：

协方差的结果有什么意义呢？如果结果为正值，则说明两者是正相关的(从协方差可以引出“相关系数”的定义)，也就是说一个人越猥琐就越受女孩子欢迎，嘿嘿，那必须的~结果为负值就说明负相关的，越猥琐女孩子越讨厌，可能吗？如果为0，也是就是统计上说的“相互独立”。

从协方差的定义上我们也可以看出一些显而易见的性质，如：

协方差多了就是协方差矩阵

上一节提到的猥琐和受欢迎的问题是典型二维问题，而协方差也只能处理二维问题，那维数多了自然就需要计算多个协方差，比如n维的数据集就需要计算个协方差，那自然而然的我们会想到使用矩阵来组织这些数据。给出协方差矩阵的定义：

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

可见，协方差矩阵是一个对称的矩阵，而且对角线是各个维度上的方差。

Matlab协方差实战

上面涉及的内容都比较容易，协方差矩阵似乎也很简单，但实战起来就很容易让人迷茫了。必须要明确一点，协方差矩阵计算的是不同维度之间的协方差，而不是不同样本之间的。这个我将结合下面的例子说明，以下的演示将使用Matlab，为了说明计算原理，不直接调用Matlab的cov函数(蓝色部分为Matlab代码)。

首先，随机产生一个10*3维的整数矩阵作为样本集，10为样本的个数，3为样本的维数。

1MySample = fix(rand(10,3)*50)

根据公式，计算协方差需要计算均值，那是按行计算均值还是按列呢，我一开始就老是困扰这个问题。前面我们也特别强调了，协方差矩阵是计算不同维度间的协方差，要时刻牢记这一点。样本矩阵的每行是一个样本，每列为一个维度，所以我们要按列计算均值。为了描述方便，我们先将三个维度的数据分别赋值：

1dim1 = MySample(:,1);2dim2 = MySample(:,2);3dim3 = MySample(:,3);

计算dim1与dim2，dim1与dim3，dim2与dim3的协方差：

 1sum( (dim1-mean(dim1)) .* (dim2-mean(dim2)) ) / ( size(MySample,1)-1 ) % 得到  74.53332sum( (dim1-mean(dim1)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -10.08893sum( (dim2-mean(dim2)) .* (dim3-mean(dim3)) ) / ( size(MySample,1)-1 ) % 得到  -106.4000

搞清楚了这个后面就容易多了，协方差矩阵的对角线就是各个维度上的方差，下面我们依次计算：

1std(dim1)^2 % 得到   108.32222std(dim2)^2 % 得到   260.62223std(dim3)^2 % 得到   94.1778

这样，我们就得到了计算协方差矩阵所需要的所有数据，调用Matlab自带的cov函数进行验证：

1cov(MySample)

把我们计算的数据对号入座，是不是一摸一样？

Update：今天突然发现，原来协方差矩阵还可以这样计算，先让样本矩阵中心化，即每一维度减去该维度的均值，使每一维度上的均值为0，然后直接用新的到的样本矩阵乘上它的转置，然后除以(N-1)即可。其实这种方法也是由前面的公式通道而来，只不过理解起来不是很直观，但在抽象的公式推导时还是很常用的！同样给出Matlab代码实现：

1X = MySample - repmat(mean(MySample),10,1);    % 中心化样本矩阵，使各维度均值为02C = (X'*X)./(size(X,1)-1);

总结

理解协方差矩阵的关键就在于牢记它计算的是不同维度之间的协方差，而不是不同样本之间，拿到一个样本矩阵，我们最先要明确的就是一行是一个样本还是一个维度，心中明确这个整个计算过程就会顺流而下，这么一来就不会迷茫了~

P.S.写论文要选Latex，在wordpress里编辑公式还得用Latex，用Latex还真对得起咱学计算机这张脸~