里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

1 PCA与人脸识别及其理论基础

龚 勋（2007-4-29）

1.1 问题描述[1]

对于一幅图像可以看作一个由像素值组成的矩阵，也可以扩展开，看成一个矢量，如一幅N*N象素的图像可以视为长度为N2的矢量，这样就认为这幅图像是位于N2维空间中的一个点，这种图像的矢量表示就是原始的图像空间，但是这个空间仅是可以表示或者检测图像的许多个空间中的一个。不管子空间的具体形式如何，这种方法用于图像识别的基本思想都是一样的，首先选择一个合适的子空间，图像将被投影到这个子空间上，然后利用对图像的这种投影间的某种度量来确定图像间的相似度，最常见的就是各种距离度量。

1.1.1 K-L变换[1]

PCA方法是由Turk和Pentlad提出来的，它的基础就是Karhunen-Loeve变换(简称KL变换)，是一种常用的正交变换。下面我们首先对K-L变换作一个简单介绍:

假设X为n维的随机变量，X可以用n个基向量的加权和来表示:

X=∑αiφi

i=1n

式中: αi是加权系数，φi是基向量，此式还可以用矩阵的形式表示:

X=(φ1,φ2, ,φn)(α1,α2, ,αn)T=Φα

取基向量为正交向量，即

1 i=j

ΦTΦj= ΦTΦj=I

0 i≠j

则系数向量为：

α=ΦTX

综上所述，K-L展开式的系数可用下列步骤求出:

T

步骤一 求随即向量X的自相关矩阵R=E XX ，由于没有类别信息的样本集的µ均值向

T

量，常常没有意义，所以也可以把数据的协方差矩阵∑=E (x µ)(x µ) 作为

K_L坐标系的产生矩阵，这里µ是总体均值向量。

步骤二 求出自相关矩阵或协方差矩阵R的本征值λi和本征向量φi，Φ=(φ1,φi, ,φn)

步骤三 展开式系数即为α=ΦTX

K_L变换的实质是建立了一个新的坐标系，将一个物体主轴沿特征矢量对齐的旋转 变换，这个变换解除了原有数据向量的各个分量之间相关性，从而有可能去掉那些带 有较少信息的坐标系以达到降低特征空间维数的目的。

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

1.1.2 利用PCA进行人脸识别

完整的PCA人脸识别的应用包括几个步骤：人脸图像预处理；读入人脸库，训练形成特征子空间；把训练图像和测试图像投影到上一步骤中得到的子空间上；选择一定 的距离函数进行识别。下面详细描述整个过程（源码见’faceRec.m’）。 1. 读入人脸库

归一化人脸库后，将库中的每人选择一定数量的图像构成训练集，其余构成测试集。设归一化后的图像是n*m,按列相连就构成N=n*m维矢量，可视为N维空间中的一个点，可以通过K-L变换用一个低维子空间描述这个图像。 2. 计算K- L变换的生成矩阵

所有训练样本的协方差矩阵为（以下三个等价）：

M

=1. C(xkixkT)/M mximxT∑A

k=1

T

（1） 2. CA=(AiA)/M

M

3.CA= ∑(xi mx)(xi mx)T i=1

A={φ1,φ2,...,φM}, φi=xi mx,mx是平均人脸, M 训练人脸数，协方差矩阵CA 是

一个N*N的矩阵, N 是xi的维数。

为了方便计算特征值和特征向量，一般选用第2个公式。根据K - L 变换原理，我们所求的新坐标系即由矩阵AiAT的非零特征值所对应的特征向量组成。直接求N*N大小矩阵CA的特征值和正交归一特征向量是很困难的, 根据奇异值分解原理（见段落1.2.5和1.2.6），可以通过求解ATiA的特征值和特征向量来获得ATiA的特征值和特征向量，。

在计算得到CA的所有非零特征值[λ0,λ1, ,λr 1]（从大到小排序，1≤r<M）及其对应的单位正交特征向量[u0,u1, ,ur 1]后，可以得到特征空间U=[u0,u1, ,ur 1]∈ N*r，从而可以计算一张图片X在特征空间上的投影系数（也可以理解为X在空间U中的坐标）：

Y=UT*X∈ r*1 （2）

3. 识别

利用公式（2），首先把所有训练图片进行投影，然后对于测试图片也进行同样的投影，采用判别函数对投影系数进行识别。

1.2 PCA的理论基础 1.2.1 投影[2]

设d维样本x1，x2， ,xn，以及一个d维基w，那么标量：

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

yi=wTxi

是相当于xi在基上的坐标值。如果w=1,yi就是把xi向方向为w的直线进行投影的结果，可以从图 1看到。推广之，如果有一组基（m个）组成的空间W=[w1,w2, ,wm]，那么可以得到xi在空间W上的坐标为：Y=WTx∈ m*1。

证明：wTxi=w x cosθ

又∵x cosθ=y, w=1

wTxi=y

图 1 投影图

进一步，表达式w=m+ae表示w是一条通过点m，方向为e的直线。

1.2.2 PCA的作用及其统计特性[3]

采用PCA对原始数据的处理，通常有三个方面的作用—降维、相关性去除、概率估计。下面分别进行介绍：

去除原始数据相关性

从统计学上讲，E{[X E(X)][Y E(Y)]}称为随机变量X与Y协方差，记为

Cov(X,Y

)。令ρXY=

，称为随机变量X与Y的相关系数。ρXY=1则X与

Y是相关的，ρXY=0，则X与Y是不相关的。

命题 1对于矩阵A来说，如果AA是一个对角阵，那么A中的向量是非相关的。

T

由PCA处理的人脸库数据的非相关性可以从两点进行说明。 （1） 基底的非相关性

特征空间基U=[u0,u1, ,ur 1]是非相关的，即UUT=I。

（2） 投影系数的非相关性

根由SVD可知A={φ1,φ2,...,φM}＝UΛVT, 其中φi=xi mx,mx是平均人脸。据公式（2）可以把A映射到特征空间上，得到： B=UT*A，其中B是非相

关的，可由下面得到证明： Y的协方差矩阵为：CB=

1112

BBT=UTAATU=Λ （3） MMM

由命题 1可知，B是非相关的。

统计参数（均值及方差）

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

均值即mx--平均人脸。

命题 2随机变量方差越大，包含的信息越多，当一个变量方差为0时，该变量为常数，不含任何信息。

用PCA计算主分量，就是寻找一组向量，使得原始数据A={φ1,φ2,...,φM}在这组向量上的投影值的方差尽可能大。

最大方差对应的向量就是第一主成份，

用PCA计算主分量就是求原始数据A={φ1,φ2

,...,φM}（其中φi=xi mx）协方差矩阵的特征向量U=[u0,u1, ,ur 1]，由公式（3）可知，P=uiTA=(p1,p2, ,pm)是A在ui上的投影值，其中P的方差就是ui对应的特征值λi 命题 3 所有原始数据在主分量ui上的投影值方差为λi。

降维

如果在原始空间表示一幅n*m大小的图片X，那么需要一个N＝n*m维矢量，但是当用公式(2)把它映射到特征空间后，只需要一个r*1维的向量就可。

另外，由命题2可知，可以根据方差的大小来判断特征向量的重要性。由ORL图片库的200个人脸计算得到的特征值呈图 2分布，可知特征向量重要性呈指数下降，据此可以只选用前面几个重要的特征向量来构建特征空间。

通过计算，前71个特征值占了90.17%，因此r可以取71而非200，从而达到进一步降维的作用。

图 2 特征值的分布

1.2.3 特征脸

U=[u0,u1, ,ur 1]中的每一个单位向量都构成一个特征脸，如图 3所示。由这些特征脸

所张成的空间称为特征脸子空间，需要注意对于正交基的选择的不同考虑，对应较大特(正交基)的特征向量则用于描述人脸的具体细节，或者从频域来看，

主分量表示了人脸的低频部分，而此分量则描述了人脸的高频部分（源码见’EigenFace.m’）。

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

1 2 10 50 70 average

图 3 特征脸，分别是第1，2，10，50，70分量，最后一张是平均脸。

1.2.4

图片重建

要进行图片X的重建，首先对X投影到特征空间上，得到系数Y=UT(X mX),然后选用一部分系数与特征向量进行原始图片的重建：X'=mX+U(1:t)*Y(1:t)，其中1:t表示取前t个元素。（见’reconstruct.m’）

在图 4中，其中前两张图片来自训练样本，第3张来自测试样本，可以看到对于训练样本，PCA系数可以对图片实现很好重建，而对于训练样本以外的图片重建效果很差。

Original1550100150199

图 4 人脸图像重建。第列张图片是输入原始图，其它列图片是重建结果，数字表示t的数目。

1.2.5 奇异值分解（SVD）[1]

设A是秩为r的m*n（m>>n）维矩阵，则存在两个正交矩阵和一个对角阵：

A=[a1,a2, ,ar]=UΛVT

其中U=[u0,u1, ,ur 1]，V=[v0,v1, ,vr 1]，Λ=diag(λ0,λ1, ,λr 1)，且UUT=I，VVT=I，

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

λi呈降序排列。其中λi2为AAT∈ m*m和ATA∈ n*n的非零特征值，ui和vi分别是AAT和ATA对应于λi2的特征向量。可得一个推论：

U=AVΛ 1

可以计算ATA的特征值λi2及相应的正交归一特征向量vi后，可由推论知AAT的正交归一特征向量

ui=

1

λi

Avi

注意，协方差矩阵CA=(AiAT)/M的特征值为：λi2/M。

1.2.6 利用小矩阵计算大矩阵特征向量

高阶矩阵的特征向量可以转化为求低阶矩阵的特征向量：

设：A是秩为r 的m*n（m>>n）维矩阵，CX=AAT∈ m*m，是一个矩阵，现在要求CX的特征值及特征向量，可通过先求小矩阵ATA∈ n*n的特征向量[v0,v1, ,vr 1]和特征值

[λ0,λ1, ,λr 1]，两者之间有以下关系：

ATA vi=λi vi

左乘A

→AAT(A vi)=λi(A vi)

显然，CX=AAT的特征向量是A vi（注意没有单位化），[λ0,λ1, ,λr 1]亦为其特征值。

结论：1.2.5与1.2.6的方法计算协方差矩阵的特征向量，特征值的结果是一致的，只是要注意1.2.5中的特征值要除以M，1.2.6中的特征向量要单位化。

1.2.7 图片归一化

图片标准化通常是一个整体概念，要求把图片归一到均值为0，方差为1下情况下。 这个概念类似于一般正态分布向标准正态分布的转化： 命题 4 若X～N(µ,σ2)，则Z=

X µ

σ

所以要对一组图片中的一张Xi进行归一化（标准化），只需要减去均值，除以方差就可以了。

T

均值mX=∑Xi，方差为D=E (Xm)(Xm) XX

iM

～N(0,1)

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

1.3 参考文献

[1] 邓楠, 基于主成份分析的人脸识别, 西北大学硕士学位论文，2006.06

[2] R.O. Duda, P.E. Hart, and D.G. Stork, Pattern Classification, seconded. John Wiley & Sons,

2001.

[3] Sami Romdhani .Face Image Analysis using a Multiple Feature Fitting Strategy. PhD Thesis,

University of Basel, Switzerland, January 2005.

[4] Sami Romdhani .FACE RECOGNITION USING PRINCIPAL COMPONENTS ANALYSIS.

1.4 附录—matlab源码 1.4.1 人脸识别

% FaceRec.m

% PCA 人脸识别修订版，识别率88％

% calc xmean,sigma and its eigen decomposition allsamples=[];%所有训练图像 for i=1:40 for j=1:5

a=imread(strcat('e:\ORL\s',num2str(i),'\',num2str(j),'.jpg')); % 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);

% 按特征值大小以降序排列 dsort = flipud(d1); vsort = fliplr(v);

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

%以下选择90%的能量 dsum = sum(dsort); dsum_extract = 0; p = 0; while( dsum_extract/dsum < 0.9) p = p + 1; dsum_extract = sum(dsort(1:p)); end i=1;

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

base = xmean' * vsort(:,1:p) * diag(dsort(1:p).^(-1/2));

% base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化（使其方差为1） % 详见《基于PCA的人脸识别算法研究》p31

% xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程 %while (i<=p && dsort(i)>0)

% base(:,i) = dsort(i)^(-1/2) * xmean' * vsort(:,i); % base是N×p阶矩阵，除以dsort(i)^(1/2)是对人脸图像的标准化（使其方差为1） % 详见《基于PCA的人脸识别算法研究》p31 % i = i + 1; % xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程 %end

% 以下两行add by gongxun 将训练样本对坐标系上进行投影,得到一个 M*p 阶矩阵allcoor allcoor = allsamples * base; % allcoor 里面是每张训练人脸图片在M*p子空间中的一个点，即在子空间中的组合系数， accu = 0; % 下面的人脸识别过程中就是利用这些组合系数来进行识别

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

a=imread(strcat('e:\ORL\s',num2str(i),'\',num2str(j),'.jpg')); 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)-1)/5 )+1; class2=floor((index2(2)-1)/5)+1; class3=floor((index2(3)-1)/5)+1;

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

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 %输出识别率

1.4.2 特征人脸

% eigface.m

function [] = eigface()

% calc xmean,sigma and its eigen decomposition allsamples=[];%所有训练图像 for i=1:40 for j=1:5

a=imread(strcat('e:\ORL\s',num2str(i),'\',num2str(j),'.jpg')); % 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);

% 按特征值大小以降序排列

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

dsort = flipud(d1); vsort = fliplr(v);

%以下选择90%的能量 dsum = sum(dsort); dsum_extract = 0; p = 0; while( dsum_extract/dsum < 0.9) p = p + 1; dsum_extract = sum(dsort(1:p)); end

p = 199;

% (训练阶段)计算特征脸形成的坐标系 %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; % xmean' * vsort(:,i)是小矩阵的特征向量向大矩阵特征向量转换的过程 %end

base = xmean' * vsort(:,1:p) * diag(dsort(1:p).^(-1/2));

% 生成特征脸 for (k=1:p),

temp = reshape(base(:,k), 112,92); newpath = ['d:\test\' int2str(k) '.jpg']; imwrite(mat2gray(temp), newpath); end

avg = reshape(samplemean, 112,92);

imwrite(mat2gray(avg), 'd:\test\average.jpg');

% 将模型保存

save('e:\ORL\model.mat', 'base', 'samplemean');

1.4.3 人脸重建

% Reconstruct.m

function [] = reconstruct()

load e:\ORL\model.mat;

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程

% 计算新图片在特征子空间中的系数 img = 'D:\test2\10.jpg' a=imread(img);

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

c = b * base; % c 是图片a在子空间中的系数, 是1*p行矢量

% 根据特征系数及特征脸重建图

% 前15个 t = 15;

temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean';

imwrite(mat2gray(reshape(temp, 112,92)),'d:\test2\t1.jpg');

% 前50个 t = 50;

temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean';

imwrite(mat2gray(reshape(temp, 112,92)),'d:\test2\t2.jpg');

% 前100个 t = 100;

temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean';

imwrite(mat2gray(reshape(temp, 112,92)),'d:\test2\t3.jpg');

% 前150个 t = 150;

temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean';

imwrite(mat2gray(reshape(temp, 112,92)),'d:\test2\t4.jpg');

% 前199个 t = 199;

temp = base(:,1:t) * c(1:t)'; temp = temp + samplemean';

imwrite(mat2gray(reshape(temp, 112,92)),'d:\test2\t5.jpg');

里面讲述了PCA用于人脸识别,附带的源码能让我们更好的了解PCA的算法过程