表1 原始数据矩阵
图1 坐标系旋转20°。虚线为原坐标系，实线为新坐标系
图2 坐标系旋转 60°。虚线为原坐标系，实线为新坐标系
表 3 根据第一主成分重构的数据
比较表1，可见重构后的数据很接近原始数据。也就是说，在允许的误差范围内，可以用 L (L<m) 维主成分重构m维变量，从而达到降维的目的。
PCA 是一种非常有效的人脸识别技术。人脸识别需要存储和传输大量数字图像，如何从海量的人脸数据中提取特征信息，去除不必要的信息或噪声是人脸识别面对的主要挑战。人脸库图像数据存储在一个n行m列的二维数据矩阵X中，利用PCA算法确定出L个主成分。这些主成分上方差之和包含了人脸的主要特征，而其它主成分上的方差之和只占总方差的很小一部分，可以忽略不计，这样就能大大节省图像存储空间和提高数据传输速度。下面通过一个计算实例来介绍一下如何用 PCA 进行人脸识别， 附录中的 MATLAB 代码供读者参考。
该例子中所采用的400张人脸照片取自人脸库https://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html共有20个人，每人10张照片（见图 3）。
图 3 人脸数据库
每张照片的分辨率为 112 x 92 个像素。首先将每张图像从二维数组转换成一维向量，依次放入一个 n 行m 列的数据矩阵 X 中。对照上节的例子，本例中 X 的行数 n = 112 x 92 = 10304，列数 m = 400。根据文献 [3]， 人脸识别算法主要有以下步骤：
1、计算平均值
2、根据方程
3、计算特征脸eigenface [4]
图4 特征脸 eigenfaces（第1~第16主成分）
4.根据特征值λ计算L值，即保留人脸主要特征所需要的最少的主成分个数。若特征值为 λ1, λ2, …, λm，则前 L各主成分方差与总方差之比为
令 r > 95%，就能够确定L的值，本例中算得L=190，也就是说，只要用190张特征脸就可以重构400个已知人脸图像了。
5.计算权重
图 5 利用特征脸重构的人脸图像
7、最后通过计算权重
图7 人脸识别结果。左为测试样本，右为被识别的已知样本
结语
机器人脸识别是一项重要的身份确认技术，近几十年来发展非常迅速，广泛应用在银行，安检，交通等场所。人脸识别是运用一定的算法，通过计算机图像处理，快速地将待测人脸和人脸数据库中的已知人脸进行比较，从而确定待测人脸和已知人脸是否属于同一个人。随着神经网络，深度学习等技术的发展，人脸识别技术日臻完善。主成分分析法是最初出现的，也是最基本的人脸识别算法，本文简单地介绍了如何应用主成分分析法进行人脸特征提取，即所谓特征脸的概念，通过具体范例演示了如何 快速进行人脸识别，并和大家分享了MATLAB程序代码，希望对这一领域感兴趣的读者有所帮助。
附录：人脸识别算例 MATLAB 代码
n = 400;
rows = 112; cols = 92; % height and width of image
S = uint8(zeros((rows*cols),n)); % Initialization of data matrix
for k = 1:40
cd(strcat('s',num2str(k)));
for j = 1:10
I = imread(strcat(num2str(j),'.pgm'));
S(:,(k-1)*10+j) = reshape(I,size(I,1)*size(I,2),1);
end
cd ..
end
ri=round(400*rand(1,1)); % Randomly pick an index.
T = S(:,ri); % T contains the image we later on will use to test the algorithm
X = S(:,[1:ri-1 ri+1:end]); % X contains the rest of the 399 images.
X = double(X);
m = mean(X,2);
centredX = X-m; % subtract mean
L = centredX'*centredX; % covariance
[V, lambda] = eig(L); % obtain eigenvalue & eigenvector
[lambda_sorted, ind] = sort(diag(lambda),'descend'); % sort eigenvalues in descending order
V_sorted = V(:, ind);
U = centredX*V_sorted; % eigenfaces
figure(1)
for i = 1:16
subplot(4,4,i)
imagesc(reshape(U(:,i),[rows, cols]))
axis image
axis off
colormap gray
end
% evaluate the number of principal components needed to represent
% 95% Total variance.
Evalsum = sum(lambda_sorted);
csum = 0;
for i = 1:n
csum = csum + lambda_sorted(i);
tv = csum / Evalsum;
if tv > 0.95
N = i;
break
end
end
% reconstruct the training images
weights = V_sorted(:,1:N);
U = centredX*weights;
Xhat = U*weights'+m;
figure(2)
for i = 1:16
subplot(4,4,i)
imagesc(reshape(Xhat(:,(i-1)*10+1),[rows, cols]))
axis image
axis off
colormap gray
end
% Recognition
T = double(T);
centredT = T-m;% Subtract the mean
s = centredT'*U(:,1:N);% calculate the weight of the test image
u = centredT*s;% project the test image into the pca space
That = u*s'+m;% reconstruct the test image
figure(3)
subplot(121)
imshow(reshape(uint8(T), [rows, cols]))
title('Test image','Fontsize',14)
subplot(122)
imagesc(reshape(That, [rows, cols]))
axis image
axis off
colormap gray
title('Reconstructed test image','Fontsize',14)
figure(4);
subplot(121)
imshow(reshape(uint8(T), [rows, cols]))
title('Looking for ...','FontWeight','bold','Fontsize',16,'color','red');
subplot(122)
z = zeros(1,size(X,2));
for i = 1:size(X,2)
z(i) = norm(weights(i,:)-s,2);
if rem(i,20) == 0
imshow(reshape(uint8(X(:,i)),112,92))
end
drawnow;
end
subplot(122)
[~, ind] = min(z);
imshow(reshape(uint8(X(:,ind)),112,92))
title('Found!','FontWeight','bold','Fontsize',16,'color','red')
参考文献：
1. Pearson, K., On lines and planes of closest fit to systems of points in space. Philosophical Magazine, 1901. 2(7-12): p. 559-572.
2. Hotelling, H., Analysis of a complex of statistical variables into principal components. Journal of Educational Psychology, 1933. 24: p. 417-441.
3. Turk, M. and A. Pentland, EIGENFACES FOR RECOGNITION. Journal of Cognitive Neuroscience, 1991. 3(1): p. 71-86.
4. Sirovich, L. and M. Kirby, LOW-DIMENSIONAL PROCEDURE FOR THE CHARACTERIZATION OF HUMAN FACES. Journal of the Optical Society of America a-Optics Image Science and Vision, 1987. 4(3): p. 519-524.
