从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.

共有31条评论

zwg说道：

2017年6月4日 15:52

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

回复
zwg说道：

2017年6月4日 15:51

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

回复
肖雅说道：

2017年4月13日 09:44

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

回复
1. John Hany说道：
  
  2017年4月17日 20:43
  
  你好，文中的方法其实完全是统计意义上的分类，只有样本在外观上很接近时才能保证有很高的精度。你所说的方法属于基于特征的分类，方法也有很多，效果一般也更好。比较经典的从五官位置辅助分类的方法比如香港中文大学的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
  
  回复
眉刘海说道：

2016年11月21日 16:30

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

回复
1. John Hany说道：
  
  2016年11月21日 16:37
  
  做过，用的是opencv官方文档里常见的Haar分类器。只不过目前OpenCV Android应用开发的代码更新工作还没有进行到人脸识别的章节。
  
  回复
bruno说道：

2016年8月27日 18:56

第一个动图好棒

回复
milliar wang说道：

2016年8月3日 19:05

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

回复
1. John Hany说道：
  
  2016年8月3日 23:32
  
  当然可以：）只要标明来源，被转载也是我的荣幸呢！
  
  回复
七月恋堇说道：

2016年7月12日 15:57

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

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

想问一下什么原因

回复
1. John Hany说道：
  
  2016年7月13日 13:39
  
  你的机器上是否正确配置了armadillo呢？另外还需要安装lapack或者intel mkl
  
  回复
renyong说道：

2016年7月4日 22:03

我没有用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。请问这是什么原因。

回复
1. John Hany说道：
  
  2016年7月5日 11:33
  
  应该是Mat下标越界导致计算失败吧，能（用断点，或vs的call stack）确定是哪一行出错的吗？
  
  回复
  1. renyong说道：
    
    2016年7月5日 15:17
    
    我设置断点单步调试在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) );
    }
    
    回复
    1. John Hany说道：
      
      2016年7月7日 13:55
      
      component_num设置为多少？它应该是大于1的
      
      回复
      1. renyong说道：
        
        2016年7月7日 14:44
        
        我设置7，不行
        
        错误信息：
        
        Unhandled exception at 0x000007FEFD82A49D in test1.exe: Microsoft C++ exception: std::logic_error at memory location 0x00000000002EDC18.
        
        这有可能是环境配置的问题吗
      2. John Hany说道：
        
        2016年7月7日 15:34
        
        数据图像的尺寸是多大的？
      3. renyong说道：
        
        2016年7月7日 17:27
        
        92×112，跟你用的是同一个人脸库啊
      4. John Hany说道：
        
        2016年7月7日 18:05
        
        而且也正确读入了？是不是前面哪个步骤出了问题，用断点看一下每一行输出矩阵的尺寸是否和公式符合。实在不行，把你的代码发到我的邮箱johnhany@163.com，我来看看
      5. renyong说道：
        
        2016年7月8日 11:23
        
        我单步调试发现应该是arma::svd_econ(U, d, V, R.t());失败，导致U，V被重置，然后出错。但我不知道为什么会分解失败，代码已经发到你的邮箱，请把我看一下
2. renyong说道：
  
  2016年7月12日 11:51
  
  我在我电脑上装上intel的MKL，然后添加到vs2013中，就可以成功运行了。我想之前的错误应该是Armadillo自带的库分解能力问题造成的。
  
  非常感谢之前的热情帮助
  
  回复
  1. John Hany说道：
    
    2016年7月13日 13:41
    
    不客气，能运行起来就好！
    
    回复
renyong说道：

2016年7月2日 16:27

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

回复
1. John Hany说道：
  
  2016年7月2日 22:42
  
  vs项目本身就是针对x64建立的，而且代码本身是可以跨平台的，只要编译器配置得当，两种平台都是可以的
  
  回复
艾尔伯特说道：

2016年6月4日 21:00

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

回复
1. John Hany说道：
  
  2016年6月14日 21:29
  
  抱歉回复的很晚，最近我在抓紧整理一批Android Studio 2上用OpenCV的教程，很快就会放出来！
  
  回复
艾尔伯特说道：

2016年6月4日 20:55

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

回复
超级赛亚人说道：

2016年5月21日 12:01

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

回复
1. John Hany说道：
  
  2016年5月21日 13:44
  
  奇异值分解在一般的矩阵理论书籍中都会有介绍，可以看成是对非方阵进行的特征值分解。通俗地讲，矩阵乘以自身的转置，这个乘积矩阵的特征值开根号就是原矩阵的奇异值。
  
  回复
周说道：

2016年5月17日 10:42

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

回复
1. John Hany说道：
  
  2016年5月17日 14:40
  
  可以先参考一下OpenCV Android SDK里面的样例（OpenCV-android-sdk\samples\face-detection）
  
  回复

JohnHany的博客

编程手绘文字

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

共有31条评论

发表评论取消回复

共有31条评论

发表评论 取消回复

发表评论取消回复