主成分分析 PCA 理解

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

先举个例子，方便理解：

1）对于一个训练集，20个sample(i=1,2,3,...,20)，特征Xi是100维[Xi1,Xi2,Xi3,...Xij,...,Xi100](j=1,2,..,100)，那么它可以建立一个20*100的样本矩阵M。

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

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

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

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

3）然后求出这个协方差矩阵Z20x20的特征值和特征向量，一般情况下应该有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
=∑biλ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时，Bn=[1 0;0 0],Bnx=[2 0]'≈[2 0.03].

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

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

%一个修改后的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);