文章首发于我的个人博客
1、PCA介绍
主成分分析(PCA)是多元统计分析同用来分析数据的一种方法。通过正交变换将一组可能存在相关性的变量转换为一组线性不想关的向量,转换后的这组变量叫做主成分。而这个转换的过程中,可以丢弃很多相关的成分或者对描述这个物体不重要的成分。从而达到对原始数据降维,提取重要特征的目的。
其方法主要是通过对协方差矩阵进行特征分解,以得出数据的主成分(即特征向量)与它们的权值(即特征值)。PCA是最简单的以特征量分析多元统计分布的方法。其结果可以理解为对原数据中的方差做出解释:哪一个方向上的数据值对方差的影响最大?换而言之,PCA提供了一种降低数据维度的有效办法;如果分析者在原数据中除掉最小的特征值所对应的成分,那么所得的低维度数据必定是最优化的(也即,这样降低维度必定是失去讯息最少的方法)。主成分分析在分析复杂数据时尤为有用,比如人脸识别。
2、PCA人脸识别
2.1 第一步
使用ORL数据库。人脸图像总样本为(Q=400),共有P(P=40)个人,每人人均L(L=10)幅图像,每张图像大小为 112 x 92 = 10304。
写出训练样本矩阵
A = ( X 1 , X 2 , . . . , X 40 ) A = (X_1,X_2,...,X_{40})A=(X1,X2,...,X40)
- 其中Xi为由第i个人的图像的每一列向量堆叠而成的一列M*N维列向量,即把矩阵向量化。需要注意的是,由于每个人都有L张图像,所以Xi对应的是每个人的平均脸,即:
X i = 1 L ∑ j = 1 L X i j X_i = \frac{1}{L}\sum_{j=1}^{L}X_{ij}Xi=L1j=1∑LXij
- 其中X i j X_{ij}Xij为第i个人第j个样本。
2.2 计算平均脸
- 计算所有训练图像的平均脸
m = 1 P ∑ i = 1 P X i m = \frac{1}{P}\sum_{i=1}^{P}X_im=P1i=1∑PXi
2.3 计算差值脸
- 计算每张人脸与平均脸的差值(去均值)。
d i = X i − m d_i = X_i - mdi=Xi−m
2.4 构造协方差矩阵
C = 1 P ∑ i = 1 P d i d i T = 1 P Y Y T C = \frac{1}{P}\sum_{i=1}^{P}d_id_i^T = \frac{1}{P}YY^TC=P1i=1∑PdidiT=P1YYT
2.5 求解协方差矩阵的特征值和特征向量,构造特征脸
如果协方差矩阵的维数为MN*NM,当维数较大时,计算量较大,所以采用奇异值分解(SVD)方法,通过求解Y T Y Y^TYYTY的特征值λ i \lambda _iλi和特征向量v i v_ivi来获得Y Y T YY^TYYT的特征值和特征向量。
将协方差矩阵的特征值按从大到小排序:λ 1 > λ 2 . . . > λ t > λ P \lambda_1 > \lambda_2 ... > \lambda_t > \lambda_Pλ1>λ2...>λt>λP。根据特征值的贡献率选取前t个最大的特征值及其对应的特征向量v i v_ivi。
其中,贡献率是指选取的特征值的和占所有特征值的比,即:
ψ = ∑ i = 1 t λ i ∑ i = 1 P λ i > a \psi = \frac{\sum_{i=1}^{t}\lambda _i }{\sum_{i=1}^{P}\lambda _i} > aψ=∑i=1Pλi∑i=1tλi>a
- 这里选取a = 0.9,得到前t个特征向量集作为特征子空间。求出原协方差矩阵的特征向量:
u i = 1 λ i A v i u_i = \frac{1}{\sqrt{\lambda_i}}Av_iui=λi1Avi
- 则特征脸空间为:
w = ( u i , u 2 , . . . , u t ) w = (u_i, u_2, ..., u_t)w=(ui,u2,...,ut)
2.6 将训练集投影到特征子空间
- 将训练集每个人脸对应的去均值脸投影到“特征脸”空间,即
Ω i = w T d i ( i = 1 , 2 , . . . 40 ) \Omega _i = w^Td_i(i=1,2,...40)Ωi=wTdi(i=1,2,...40)
2.7 识别阶段
- 将待识别的人脸图像$\Gamma $与平均脸的差值脸投影到特征空间,得到其特征向量表示:
Ω T = w T ( Γ − m ) \Omega ^T = w^T(\Gamma -m)ΩT=wT(Γ−m)
- 使用欧氏距离来识别人脸,将待识别人脸对应特征子空间的向量与训练集中每个人脸的特征脸计算欧氏距离,选择向量距离最小值对应的类别,即识别属于某一类别的人脸:
θ = 1 2 m a x { ∣ ∣ Ω i − Ω j ∣ ∣ } \theta = \frac{1}{2}max\{||\Omega_i-\Omega_j||\}θ=21max{∣∣Ωi−Ωj∣∣}
3、通过矩阵维度变换说明算法的流程
这部分,我将根据以上PCA的步骤以矩阵维度变换的形式,直观说明PCA的流程。
- 首先,我们得到的训练样本矩阵A的维度为(10304,40),所有人脸的均值向量m的维度为(10304,1)。由此得到的去均值人脸d i d_idi数据的维度为(10304,40)。
- 其次,构造的协方差矩阵Y Y T YY^TYYT,其维度为:(10304,40)x (40,10304)=(10304,10304).但是,这个矩阵维度太高,计算量极大。因此,我们通过计算Y T Y Y^TYYTY求特征值和特征向量,维度为:(40,10304)x (10304,40) = (40,40)。最终求得特征值维度为(40,1),特征向量维度为(40,40)。
- 选取能量值占所有特征值90%的前t个特征向量组成特征子空间,其维度为(40,t)。求得特征脸空间w ww,维度为:(10304,40)x (40,t)= (10304,t)。
- 将每一幅人脸与平均脸得差值脸适量投影到特征子空间w T d i w^Td_iwTdi:(t,10304)x (10304,40)= (t,40).
- 将待识别人脸与平均脸得差值脸投影到特征空间:(t,10304)x (10304,1)= (t,1).
- 计算(t,40)和(t,1)对应的40个欧氏距离,最小值对应的类别即为该类别人脸。
4、实验结果
4.1 特征脸

4.2 人脸重建
t=16是前t个特征值的和占所有特征值和的前90%取值。
4.3 人脸识别
我实现的PCA人脸识别代码准确率达到71.25%。
5、代码
void PCA_achieve_2() {
//ORL训练集路径
std::string train_image_path = "../images/ORL_faces_train_3/s";
int per_class_num = 8; //每一个类别的样本数
int class_nums = 40; //类别数
//定义一个(10304 x 40*per_class_num)的数组存放所有训练图像
cv::Mat all_train_image(10304, 40*per_class_num, CV_32FC1, cv::Scalar(0));
std::cout << all_train_image.rows << std::endl;
//均值脸(10304 x 1)
cv::Mat all_train_image_mean;
//每一个人脸数据集合
std::vector<cv::Mat> per_train_image(40);
int n_col = 0;
for (int i = 1; i <= 40; i++) {
std::string train_image_sub_path = train_image_path + std::to_string(i);
cv::Mat per_image(112, 92, CV_32FC1, cv::Scalar(0));
for (int j = 1; j <= per_class_num; j++) {
std::string image_path = train_image_sub_path + "/" + std::to_string(j) + ".png";
std::cout << image_path << std::endl;
auto image = cv::imread(image_path, 0);
image.convertTo(image, CV_32FC1, 1.0 / 255);
//cv::Rect()只能获取对象的拷贝,而不能获取对象的引用。operator=()重载应该是引用,所以用+=即可完成赋值。
all_train_image(cv::Rect(n_col++, 0, 1, 10304)) += image.reshape(0, 10304);
}
}
//求所有图像的均值,得到均值脸
cv::reduce(all_train_image, all_train_image_mean, 1, cv::REDUCE_AVG);
//从all_train中取出每一个person的数据,计算均值.
for (int i = 0; i < 40; i++) {
cv::reduce(all_train_image(cv::Rect(i * train_num, 0, train_num, 10304)), per_train_image[i], 1, cv::REDUCE_AVG);
}
//去中心化
cv::Mat decenter_image(10304, 40, 5, cv::Scalar(0));
for (int i = 0; i < 40; i++) {
per_train_image[i] -= all_train_image_mean;
}
for (int i = 0; i < 40; i++) {
decenter_image(cv::Rect(i, 0, 1, 10304)) += per_train_image[i];
}
//decenter_image:(10304 x 40)
//协方差矩阵,然而直接使用(10304 x 40)的维度太高,直接计算特征值与特征向量太耗时
//使用其转置计算。得到的特征值和特征向量相同
cv::Mat temp_L;
cv::mulTransposed(decenter_image, temp_L, true);
std::cout << "dim:" << temp_L.size();
cv::Mat eigen_values;
cv::Mat eigen_vectors;
cv::eigen(temp_L, eigen_values, eigen_vectors);
std::cout << "eigen_values:" << eigen_values;
std::cout << "eigen_values_size:" << eigen_values.size();
std::cout << "eigen_values_type:" << eigen_values.type();
std::cout << "eigen_vector_size:" << eigen_vectors.size();
std::cout << "eigen_vector_type:" << eigen_vectors.type();
double eigen_value_sum = 0;
for (int j = 0; j < eigen_values.rows; j++) {
eigen_value_sum += eigen_values.at<float>(j, 0);
std::cout << eigen_values.at<float>(j, 0) << std::endl;
}
float lamb = 0;
int t = 0;
std::cout << "eigen_value_sum:" << eigen_value_sum << std::endl;
for (; t < eigen_values.rows; t++) {
lamb += eigen_values.at<float>(t, 0);
std::cout << lamb / eigen_value_sum << std::endl;
if (lamb / eigen_value_sum >= 0.9) {
break;
}
}
std::cout << "get_t:" << t << std::endl;
//L_eigen_vec : (40 x t)
cv::Mat L_eigen_vec(40, t, 5, cv::Scalar(0));
for (int i = 0; i < t; i++) {
L_eigen_vec(cv::Rect(i, 0, 1, 40)) += eigen_vectors(cv::Rect(i, 0, 1, 40));
}
//计算特征脸 decenter_image x L_eigen_vec : (10304 x 40) x (40 x t)
cv::Mat eigen_faces = decenter_image * L_eigen_vec;
/*******人脸重建******/
std::string test_image_class_10_path = "../images/ORL_faces_test/s20/9.png";
cv::Mat test_image_class_10 = cv::imread(test_image_class_10_path, 0);
cv::imshow("origin_f", test_image_class_10);
// (10304 x 1)
test_image_class_10 = test_image_class_10.reshape(0, 10304);
test_image_class_10.convertTo(test_image_class_10, 5, 1.0 / 255);
//去均值
cv::Mat decenter_test_image_class_10 = test_image_class_10 - all_train_image_mean;
//将样本映射到特征子空间(t x 10304) x (10304 x 1) = (t x 1)
cv::Mat Y_test = cv::Mat(eigen_faces.t()) * decenter_test_image_class_10;
Y_test.convertTo(Y_test, 5, 1.0 / 255);
//样本重建(10304 x t) x (t x 1)
cv::Mat rebuild_image = eigen_faces * Y_test;
rebuild_image += all_train_image_mean;
rebuild_image.convertTo(rebuild_image, 0, 255);
cv::imshow("rebuild_image", rebuild_image.reshape(0, (92, 112)));
cv::waitKey(0);
/***人脸识别部分***/
//将训练集投影到特征子空间
//(t x 10304)*(10304 x 1)
std::vector<cv::Mat> train_sub;
for (auto iter = per_train_image.begin(); iter != per_train_image.end(); iter++) {
cv::Mat train_c;
cv::Mat((eigen_faces.t())* (*iter)).convertTo(train_c, 5, 1.0 / 255);
train_sub.push_back(train_c);
}
//读取测试集图像
std::string test_image_path = "../images/ORL_faces_test_3/s";
int right = 0;
for (int i = 1; i <= 40; i++) {
std::string test_image_sub_path = test_image_path + std::to_string(i);
cv::Mat per_image(112, 92, CV_32FC1, cv::Scalar(0));
for (int j = 0; j < 2; j++) {
std::string image_path = test_image_sub_path + "/" + std::to_string(j + 9) + ".png";
cv::Mat test_image = cv::imread(image_path, 0);
test_image = test_image.reshape(0, 10304);
test_image.convertTo(test_image, 5, 1.0 / 255);
test_image -= all_train_image_mean;
cv::Mat test_image_sub = cv::Mat(eigen_faces.t()) * test_image;
test_image_sub.convertTo(test_image_sub, 5, 1.0 / 255);
int result = compute_face_dist(train_sub, test_image_sub);
if (result == i)
right++;
}
}
std::cout << t << " recg right :" << right << "recg accracy:" << (float)(right / 80.0) << std::endl;
return;
}
6、参考资料
[1] M. A. Turk and A. P. Pentland, “Face recognition using eigenfaces,” Proceedings. 1991 IEEE Computer Society Conference on Computer Vision and Pattern Recognition, Maui, HI, USA, 1991, pp. 586-591, doi: 10.1109/CVPR.1991.139758.