从QR分解到PCA，再到人脸识别

PCA（Principal Component Analysis，主成分分析）是一种很常用的根据变量协方差对数据进行降维、压缩的方法。它的精髓在于尽量用最少数量的维度，尽可能精确地描述数据。

PCA对数据进行降维的过程可以用下面这个动图来解释（图片摘自http://stats.stackexchange.com/a/140579/93946）：

在上图中，一组位于直角坐标系的二维样本集，沿着斜线的方向有很强的相关性。所以如果我们将直角坐标系转换到斜向，也就是让横轴沿斜线方向，纵轴垂直于斜线方向。于是，在这个新的坐标系下，数据点在横轴上分布很分散，但是在纵轴方向比较集中。如果在误差允许范围内，我们完全可以将数据点在新纵轴上的坐标全部置为0，只用新横轴上的坐标来表示点的位置。这样，就完成了对数据降维的过程（即将原始直角坐标系的2个维度减小到新坐标系的1个维度）。对更高维的情况，处理过程与之类似。

PCA人脸识别

将PCA用于人脸识别的过程如下：

1.假设有400幅尺寸为100*100的图像，构成10000*400的矩阵 $X=[x_1,\dots,x_n]$ ；

2.计算均值 $\mu=\frac{1}{n}\sum_{j=1}^n x_j$ ，令 $H=\frac{1}{\sqrt{n-1}}[x_1-\mu,\dots,x_n-\mu]$ ；

3.根据定义，计算协方差矩阵 $\Sigma=HH^T$ ；

4.计算 $\Sigma$ 的特征值与特征向量，取前h个最大特征值所对应的特征向量，构成矩阵 $\Phi$ ；

5.矩阵 $\Phi$ 可对数据降维： $\Phi^T X=Y$ ，Y是h行400列的矩阵，也就是将数据从10000维降为h维。

这种做法一个明显的缺陷在于， $\Sigma$ 的维度为10000×10000，直接进行奇异值分解计算量非常大。利用QR分解，作间接的奇异值分解，可以减小计算量。

利用QR分解减小计算量

基于QR分解的PCA算法步骤如下：

1.已知 $\Sigma=HH^T$ ，其中 $\Sigma$ 为d*d，H为d*n，d代表原始数据的维数，n代表样本数，d远大于n；

2.对H作QR分解， $h=QR$ ，其中Q为d*t，R为t*n， $1\leq t \leq n$ ；

3.则 $\Sigma=QRR^T Q^T$ ，对 $R^T$ 作奇异值分解 $R^T=UDV^T$ ，其中U为n*t，V为t*t， $D=diag(\sigma_1,\dots,\sigma_t)$ ；

4.于是 $\Sigma=QVDU^T UDV^T Q^T=QVD^2 V^T Q^T=QV\Lambda V^T Q^T$ ，其中 $\Lambda=D^2$ ；

5.由于 $(QV)^T (QV)=V^T Q^T QV=V^T V=I$ ，所以QV可将 $\Sigma$ 对角化，QV为 $\Sigma$ 的特征向量矩阵， $\Lambda$ 为 $\Sigma$ 的特征值矩阵；

6.选取D前h个最大对角元所对应于V中的h个列，构成t*h的矩阵 $V_h$ ，则降维矩阵 $\Phi=QV_h$ ；

该过程涉及对一个很大的矩阵的QR分解，和对一个较小矩阵的奇异值分解。计算量与传统方法相比较的结果如下（图片摘自[1]）：

进一步，进行人脸识别的过程如下：

1.假设有c个类别，每类包含s个样本，则n=c∗s；

2.对X计算 $Y=\Phi^T X$ ，将Y（也称特征脸）按类别计算均值，得到c个长度为h的列向量 $v_1,\dots,v_c$ ；

3.对于未知类别的新样本x，计算 $y=\Phi^T x$ ，y的长度为h；

4.计算距离 $d(y,v_i),i=1,\dots,c$ ，取距离最小的i作为x的类标号。

距离度量d

1.欧式距离

2.曼哈顿距离

3.马氏距离

在马氏距离中，x与y分布相同，且协方差矩阵为S。加入协方差矩阵的逆矩阵的作用是，将如下图（图片部分取自http://stats.stackexchange.com/a/62147/93946）中呈椭圆分布的数据归一化到圆形分布中，再来比较距离，可以抵消不同样本集在特征空间中的分布差异。

C++实现

环境要求：OpenCV（样本图像的读取），Armadillo（高性能线性代数C++函数库），Intel MKL（替代LAPACK为Armadillo提供矩阵分解计算）。

项目使用Visual Studio Ultimate 2012建立，不过核心代码只有一个cpp文件。

全部代码托管在github.com/johnhany/QR-PCA-FaceRec。

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

/************************************************************************/

/* QR-PCA-FaceRec by John Hany */

/* */

/* A face recognition algorithm using QR based PCA. */

/* */

/* Released under MIT license. */

/* */

/* Contact me at johnhany@163.com */

/* */

/* Welcome to my blog http://johnhany.net/, if you can read Chinese:) */

/* */

/************************************************************************/

#include <opencv2/opencv.hpp>

#include <armadillo>

#include <iostream>

using namespace std;

int component_num = 7;

string orl_path = "G:\\Datasets\\orl_faces";

enum distance_type {ECULIDEAN = 0, MANHATTAN, MAHALANOBIS};

//double distance_criterion[3] = { 10.0, 30.0, 3.0};

double distance_criterion[3] = { 1000.0, 1000.0, 1000.0};

bool compDistance(pair<int, double> a, pair<int, double> b);

double calcuDistance(const arma::vec vec1, const arma::vec vec2, distance_type dis_type);

double calcuDistance(const arma::vec vec1, const arma::vec vec2, const arma::mat cov2, distance_type dis_type);

int main(int argc, const char *argv[]) {

int class_num = 40;

int sample_num = 10;

int img_cols = 92;

int img_rows = 112;

cv::Size sample_size(img_cols, img_rows);

arma::mat mat_sample(img_rows*img_cols, sample_num*class_num);

//Load samples in one matrix `mat_sample`.

for(int class_idx = 0; class_idx < class_num; class_idx++) {

for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {

string filename = orl_path + "\\s" + to_string(class_idx+1) + "\\" + to_string(sample_idx+1) + ".pgm";

cv::Mat img_frame = cv::imread(filename, CV_LOAD_IMAGE_GRAYSCALE);

cv::Mat img_sample;

cv::resize(img_frame, img_sample, sample_size);

for(int i = 0; i < img_rows; i++) {

uchar* pframe = img_sample.ptr<uchar>(i);

for(int j = 0; j < img_cols; j++) {

mat_sample(i*img_cols+j, class_idx*sample_num+sample_idx) = (double)pframe[j]/255.0;

}

// cout << mat_sample.n_rows << endl << mat_sample.n_cols << endl << mat_sample(img_rows*img_cols/2, 0) << endl;

//Calculate PCA transform matrix `mat_pca`.

arma::mat H = mat_sample;

arma::mat mean_x = arma::mean(mat_sample, 1);

for(int j = 0; j < class_num * sample_num; j++) {

H.col(j) -= mean_x.col(0);

}

H *= 1.0/sqrt(sample_num-1);

arma::mat Q, R;

arma::qr_econ(Q, R, H);

arma::mat U, V;

arma::vec d;

arma::svd_econ(U, d, V, R.t());

// cout << "d" << endl << d << endl;

// arma::rowvec vec_eigen = d.head(component_num).t();

// cout << "vec_eigen" << endl << vec_eigen << endl;

arma::mat V_h(V.n_rows, component_num);

if(component_num == 1) {

V_h = V.col(0);

}else {

V_h = V.cols(0, component_num-1);

}

arma::mat mat_pca = Q * V_h;

//Calculate eigenfaces `mat_eigen_vec`.

arma::mat mat_eigen = mat_pca.t() * mat_sample;

// cout << "mat_eigen" << endl << mat_eigen << endl;

arma::mat mat_eigen_vec(component_num, class_num, arma::fill::zeros);

vector<arma::mat> mat_cov_list;

for(int class_idx = 0; class_idx < class_num; class_idx++) {

arma::vec eigen_sum(component_num, arma::fill::zeros);

for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {

eigen_sum += mat_eigen.col(class_idx*sample_num+sample_idx);

}

eigen_sum /= (double)sample_num;

mat_eigen_vec.col(class_idx) = eigen_sum;

mat_cov_list.push_back(arma::cov((mat_eigen.cols(class_idx*sample_num, class_idx*sample_num+sample_num-1)).t()));

// cout << mat_cov_list[class_idx] << endl;

}

// cout << "mat_eigen_vec" << endl << mat_eigen_vec << endl;

cout << "dis within class" << endl;

for(int class_idx = 0; class_idx < class_num; class_idx++) {

for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {

double dis = calcuDistance(mat_eigen.col(class_idx*sample_num+sample_idx), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS);

cout << dis << " ";

}

cout << endl;

}

cout << "dis between classes" << endl;

for(int class_idx = 0; class_idx < class_num; class_idx++) {

for(int sample_idx = 0; sample_idx < class_num; sample_idx++) {

double dis = calcuDistance(mat_eigen.col(sample_idx*sample_num), mat_eigen_vec.col(class_idx), mat_cov_list[class_idx], distance_type::MAHALANOBIS);

cout << dis << " ";

}

cout << endl;

}

//Classify new sample.

int correct_count = 0;

double max_dis = 0.0;

for(int class_idx = 0; class_idx < class_num; class_idx++){

for(int sample_idx = 0; sample_idx < sample_num; sample_idx++) {

arma::mat mat_new_sample(img_rows*img_cols, 1);

string filename = orl_path + "\\s" + to_string(class_idx+1) + "\\" + to_string(sample_idx+1) + ".pgm";

cv::Mat img_new_frame = cv::imread(filename, CV_LOAD_IMAGE_GRAYSCALE);

cv::Mat img_new_sample;

cv::resize(img_new_frame, img_new_sample, sample_size);

for(int i = 0; i < img_rows; i++) {

uchar* pframe = img_new_sample.ptr<uchar>(i);

for(int j = 0; j < img_cols; j++) {

mat_new_sample(i*img_cols+j, 0) = (double)pframe[j]/255.0;

}

arma::mat mat_new_eigen = mat_pca.t() * mat_new_sample;

vector<pair<int, double>> dis_list;

for(int new_class_idx = 0; new_class_idx < class_num; new_class_idx++) {

double dis = calcuDistance(mat_new_eigen.col(0), mat_eigen_vec.col(new_class_idx), mat_cov_list[new_class_idx], distance_type::MAHALANOBIS);

dis_list.push_back(make_pair(new_class_idx, dis));

}

sort(dis_list.begin(), dis_list.end(), compDistance);

if(dis_list[0].first == class_idx && dis_list[0].second <= distance_criterion[distance_type::MAHALANOBIS]) {

correct_count++;

}

if(dis_list.back().second > max_dis) {

max_dis = dis_list.back().second;

}

cout << "Maximum distance: " << max_dis << endl;

double correct_ratio = (double)correct_count / (class_num * sample_num);

cout << "Correctness ratio: " << correct_ratio * 100.0 << "%" << endl;

cin.get();

return 0;

}

bool compDistance(pair<int, double> a, pair<int, double> b) {

return (a.second < b.second);

}

double calcuDistance(const arma::vec vec1, const arma::vec vec2, distance_type dis_type) {

if(dis_type == ECULIDEAN) {

return arma::norm(vec1-vec2, 2);

}else if(dis_type == MANHATTAN) {

return arma::norm(vec1-vec2, 1);

}else if(dis_type == MAHALANOBIS) {

arma::mat tmp = (vec1-vec2).t() * (vec1 - vec2);

return sqrt(tmp(0,0));

}

return -1.0;

}

double calcuDistance(const arma::vec vec1, const arma::vec vec2, const arma::mat cov2, distance_type dis_type) {

if(dis_type == ECULIDEAN) {

return arma::norm(vec1-vec2, 2);

}else if(dis_type == MANHATTAN) {

return arma::norm(vec1-vec2, 1);

}else if(dis_type == MAHALANOBIS) {

arma::mat tmp = (vec1-vec2).t() * cov2.i() * (vec1 - vec2);

return sqrt(tmp(0,0));

}

return -1.0;

}

分类测试

测试样本采用The AT&T Database of Faces，原称The ORL Database of Faces，包含取自40个人的样本，每人10幅，共400幅图像，图像尺寸92*112。

样本库的链接为http://www.cl.cam.ac.uk/research/dtg/attarchive/facedatabase.html。

1.欧式距离降维及分类效果：

即使将h设为最大的400（样本数），其分类正确率也只能达到99%。

2.曼哈顿距离降维及分类效果：

在h为128时，分类正确率可以达到100%，降维能力略好于欧式距离。

3.马氏距离降维及分类效果：

在h仅仅为7的时候，其分类正确率就已经达到100%了。采用马氏距离的PCA方法可以将人脸数据的维度从10000左右降到7，降维效果非常优秀。在h超过9时，分类过程中计算的最大马氏距离超出了双精度浮点数double的上限。

参考文献

[1] Sharma A, Paliwal K K, Imoto S, et al. Principal component analysis using QR decomposition[J]. International Journal of Machine Learning and Cybernetics, 2013, 4(6): 679-683.

[2] Raj D. A Realtime Face Recognition system using PCA and various Distance Classifiers[J]. 2011.

[3] Turk M, Pentland A. Eigenfaces for recognition[J]. Journal of Cognitive Neuroscience, 1991, 3(1): 71-86.

最新最旧

孙泽金

您好，如果我直接对QR分解完后的R矩阵进行PCA处理是不是理论上也是可行的呢

zwg

您好！71行代码H /= 1.0 / sqrt(sample_num – 1);是不是应该为class_num*sample_num?

肖雅

你好，我想请问一下，你这篇博客的思想是直接对整张图片进行匹配是吧。这样的话当人脸角度发生变化，或者拍照的环境变化时，结果应该会受到比较大的影响，是不是去分别比对眼睛，鼻子，嘴巴，结果会更加可信呢？

John Hany

你好，文中的方法其实完全是统计意义上的分类，只有样本在外观上很接近时才能保证有很高的精度。你所说的方法属于基于特征的分类，方法也有很多，效果一般也更好。比较经典的从五官位置辅助分类的方法比如香港中文大学的DeepID（2014年）：http://mmlab.ie.cuhk.edu.hk/pdf/YiSun_CVPR14.pdf，他们所采用的五官定位算法：http://www.cv-foundation.org/openaccess/content_cvpr_2013/papers/Sun_Deep_Convolutional_Network_2013_CVPR_paper.pdf。

当然，人脸姿态变化比较大时，五官的变化也很大甚至发生遮挡。最近很流行的思路是用GAN（及其变种）来训练生成网络，顺便也能获得一个高精度的分类器，比如上礼拜中科院的TP-GAN：https://arxiv.org/pdf/1704.04086.pdf

眉刘海

楼主，有做过Android opencv的人脸识别么？

做过，用的是opencv官方文档里常见的Haar分类器。只不过目前OpenCV Android应用开发的代码更新工作还没有进行到人脸识别的章节。

bruno

第一个动图好棒

milliar wang

楼主，可以转载你的这篇文章吗？

当然可以：）只要标明来源，被转载也是我的荣幸呢！

七月恋堇

你好，我用的vs2010+opencv2.4.9，编译时出现

fatal error C1083: 无法打开包括文件:“armadillo”: No such file or directory

想问一下什么原因

你的机器上是否正确配置了armadillo呢？另外还需要安装lapack或者intel mkl

renyong

我没有用Intel MKL，而是用LAPACK为Armadillo提供矩阵分解计算。编译通过了，运行时提示warning：svd():decomposition failed，error:Mat::cols(): indices out of bounds or incorrectly used .测试时arma::svd_econ(U, d, V, R.t());返回0。请问这是什么原因。

应该是Mat下标越界导致计算失败吧，能（用断点，或vs的call stack）确定是哪一行出错的吗？

我设置断点单步调试在78行时arma::svd_econ(U, d, V, R.t());svd_econ()返回值已经为0了，

运行到89行：V_h = V.cols(0, component_num-1);程序出错

后跳转到了debug.hpp中的这一段代码template<typename T1>
arma_cold
arma_noinline
static
void
arma_stop_logic_error(const T1& x)
{
#if defined(ARMA_PRINT_ERRORS)
{
get_stream_err1() << "\nerror: " << x << std::endl;
}
#endif

throw std::logic_error( std::string(x) );
}

component_num设置为多少？它应该是大于1的

我设置7，不行

错误信息：

Unhandled exception at 0x000007FEFD82A49D in test1.exe: Microsoft C++ exception: std::logic_error at memory location 0x00000000002EDC18.

这有可能是环境配置的问题吗

数据图像的尺寸是多大的？

92×112，跟你用的是同一个人脸库啊

而且也正确读入了？是不是前面哪个步骤出了问题，用断点看一下每一行输出矩阵的尺寸是否和公式符合。实在不行，把你的代码发到我的邮箱johnhany@163.com，我来看看

我单步调试发现应该是arma::svd_econ(U, d, V, R.t());失败，导致U，V被重置，然后出错。但我不知道为什么会分解失败，代码已经发到你的邮箱，请把我看一下

我在我电脑上装上intel的MKL，然后添加到vs2013中，就可以成功运行了。我想之前的错误应该是Armadillo自带的库分解能力问题造成的。

非常感谢之前的热情帮助

不客气，能运行起来就好！

刘同学

你好，我在我的vs上运行出现跟你一样的错误，我已经下载mkl并添加到vs中，还没有运行成功，求详细的解决办法？

您好，您这个程序是基于arm的吗？如果是，有X64下的吗

vs项目本身就是针对x64建立的，而且代码本身是可以跨平台的，只要编译器配置得当，两种平台都是可以的

艾尔伯特

能给个qq或者微信交流吗，如果可以发我邮箱，谢谢了

抱歉回复的很晚，最近我在抓紧整理一批Android Studio 2上用OpenCV的教程，很快就会放出来！

您好，我正在进行android平台上的人脸识别开发，我用JNI的方法，需要把识别后的图像（Mat类型）存入jintArray型数组，再把jintArray传回java层显示，你有什么好办法啊

超级赛亚人

这数学知识用的深啊~ 奇异值分解什么的矩阵里面的内容吗？

奇异值分解在一般的矩阵理论书籍中都会有介绍，可以看成是对非方阵进行的特征值分解。通俗地讲，矩阵乘以自身的转置，这个乘积矩阵的特征值开根号就是原矩阵的奇异值。

周

请问有木有android opencv 人脸识别的教程

可以先参考一下OpenCV Android SDK里面的样例（OpenCV-android-sdk\samples\face-detection）

wpDiscuz

Peanut 在 Windows7+VS2012下64位OpenCV3.0+CUDA7.5的编译和部署在确认一切配置都是正确的之后，我将OpenCV版本从3.4. ...
Dominic 在 N条线段求交的扫描线算法这是属于退化情况的一种吧，扫描线算法不适合求这样的退化情况。 ...
John Hany 在 Android Studio 2.3利用CMAKE进行OpenCV 3.2的NDK开发可能你在Java代码里显式地调用了opencv，需要按照这篇 ...
yangrise 在 Android Studio 2.3利用CMAKE进行OpenCV 3.2的NDK开发博主您好，请问出现opencv manager packag ...
Ling 在 OpenCV实现基于傅里叶变换的旋转文本校正谢谢博主*v*
John Hany 在 Windows7+VS2012下64位OpenCV3.0+CUDA7.5的编译和部署从这个错误信息还是不好判断根本原因出在哪里，这个错误其实也是 ...
王小我在 Windows7+VS2012下64位OpenCV3.0+CUDA7.5的编译和部署请教大神再cmake成功后，尝试编译opencv-core时 ...
John Hany 在 Tracking-Learning-Detection原理分析在那一段话的紧接着下面两段文字解释的就是如何定义蕨。在那副数 ...

JohnHany的博客

编程手绘文字

从QR分解到PCA，再到人脸识别

33