在计算协方差和特征向量的方法上，书上使用的是一种被作者称为compacttrick的技巧，以及奇异值分解（SVD），这些都是什么东西呢？
如何把PCA运用在多张图片上？
所以，我们需要进一步的了解，同时，为示例对多张图片进行PCA，我选了一个跟书相似但更有趣的例子来做——人脸识别。
一个特征脸（Eigenface，也叫标准脸）其实就是从一组人脸图像应用PCA获得的主成分特征向量之一，下面我们能验证，每个这样的特征向量变换为二维图像后看起来有点像人脸，所以才被称为特征脸，计算特征脸是进行人脸识别的首要步骤，其计算过程其实就是PCA。
特征脸计算步骤
准备一组（假设10张）具有相同分辨率（假设:100×100）的人脸图像，把每张图像打平（numpy.flatten）成一个向量，即所有像素按行串联起来,每个图像被当作是10000维空间的一个点。再把所有打平的图像存储在10000×10矩阵X中，矩阵的每一列就是一张图片，每个维为一行，共10000维
对X的每一维（行）求均值，得到一个10000×1的向量，称为平均脸（因为把它变换为二维图像看起来像人脸），然后将X减去平均脸，即零均值化
计算X的协方差矩阵C=XX^(X^表示X的转置)
计算C的特征值和特征向量，这组特征向量就是一组特征脸。在实际应用中，我们只需要保留最主要的一部分特征向量做为特征脸即可。
有了上个笔记的基础知识后，上面计算过程不难理解。但在实现代码之前，我们先来看看上面提到的计算X的协方差矩阵C=XX^引发的一个问题。
对于上面举例的矩阵X，它有10000行（维），它的协方差矩阵将达到10000×10000，有10000个特征向量，这个计算量是很大的，消耗内存大，我尝试过不可行。
这个数学问题终究还是被数学解决了，解决的方法可以见维基百科特征脸内容描述。简言之，就是把本来计算XX^的协方差矩阵（设为C）变换为计算X^X的协方差矩阵（C'），后者的结果是10×10（10为样本数量），很快就可以算出来。当然，通过这个变换分别算出来的特征向量不是等价的，也需要变换一下：设E为从C算出来的特征向量矩阵，E'为从C'算出来的特征向量矩阵，则E=XE'，最后再把E归一化。
这个技巧就是书上PCA示例使用的，被称为compacttrick的方法。
但要看明白书上的示例代码，还要搞清一点：
对原始图像数据集矩阵的组织方式，我们用行表示维度，列表示样本，而书上和《GuidetofacerecognitionwithPython》（见底部参考链接）使用的是行表示样本，列表示维度。就是因为这两种组织方式的不同，导致了PCA算法的代码看起来有些不同。这一点很容易让人困惑，所以写到这里，我应该特别的强调一下。
我之所以在上个笔记，包括上面对特征脸的计算步骤描述，都认定以行表示维度，列表示样本的方式，是为了与数学原理的详解保持一致（注：下面的代码示例还是使用这种方式），当我们明白了整个原理之后，我们便知道使用这两种矩阵表达方式都可以，两者实现的PCA代码差别也很小，下面会讲到。
网上有不少用于研究的人脸数据库可以下载，我在参考链接给出了常被使用的一个。下载解压后在目录orl_faces下包含命名为s1,..,s40共40个文件夹，每个文件夹对应一个人，其中存储10张脸照，所有脸照都是92×112的灰度图，我把部分照片贴出来：
接下来，我们按照特征脸计算步骤中的第1点所述，把这400张图像组成矩阵（图像组织不分先后），代码：
defgetimpaths(datapath):paths=[]fordirinos.listdir(datapath):try:forfilenameinos.listdir(os.path.join(datapath,dir)):paths.append(os.path.join(datapath,dir,filename))except:passreturnpathsimpaths=getimpaths('./orl_faces')m,n=np.array(Image.open(impaths[0])).shape[0:2]#图片的分辨率，下面会用到X=np.mat([np.array(Image.open(impath)).flatten()forimpathinimpaths]).Tprint'X.shape=',X.shape#X.shape=(10304,400)我们把每个图像都打平成行向量，把所有图像从上到下逐行排列，最后转置一下，便得到一个10304×400的矩阵，其中10304=92×112
PCA函数我尽量使用与书上相同的变量命名，方便对比：
defpca(X):dim,num_data=X.shape#dim:维数，num_data:样本数mean_X=X.mean(axis=1)#求出平均脸，axis=1表示计算每行（维）均值，结果为列向量X=X-mean_X#零均值化M=np.dot(X.T,X)#使用compacttrick技巧，计算协方差e,EV=np.linalg.eigh(M)#求出特征值和特征向量print'e=',e.shape,eprint'EV=',EV.shape,EVtmp=np.dot(X,EV).T#因上面使用了compacttrick，所以需要变换print'tmp=',tmp.shape,tmpV=tmp[::-1]#将tmp倒序，特征值大的对应的特征向量排前面，方便我们挑选前N个作为主成分print'V=',V.shape,Vforiinrange(EV.shape[1]):V[:,i]/=np.linalg.norm(EV[:,i])#因上面使用了compacttrick，所以需要将V归一化returnV,EV,mean_X执行PCA并画图对上面得到的X矩阵调用pca函数，并画出平均脸和部分特征脸：
V,EV,immean=pca(X)#画图plt.gray()plt.subplot(2,4,1)#2行4列表格，第一格显示`平均脸`plt.imshow(immean.reshape(m,n))#以下选前面7个特征脸，按顺序分别显示到其余7个格子foriinrange(7):plt.subplot(2,4,i+2)plt.imshow(V[i].reshape(m,n))plt.show()显示效果图如下：
希望不会被这些特征脸吓到：）这些所谓的脸事实上是特征向量，只不过维数与原始图像一致，因此可以被变换成图像显示出来，不同的特征脸代表了与均值图像差别的不同方向。当然，我们求特征脸，并不是为了显示他们，而是保留部分特征脸来获得大多数脸的近似组合。因此，人脸便可通过一系列向量而不是原始数字图像进行保存，节省了很多存储空间，也便于后续的识别计算。
与书上pca的实现对比上面我给出的pca函数代码，是按照我们一路学习PCA的思路写出来的，虽然跟书上pca实现很接近，但还是有几个点值得分析：
如上提到，我们对X矩阵的组织是以行为维、列为样本的方式，即一个列对应一张打平的图像，而书上的例子是以行为样本、列为维的方式，每一行对应一张打平的图像，而且参考链接里的例子也都是以后者进行组织的，但没关系，我们只需要对上面的代码作一点修改即可：
defpca_book(X):num_data,dim=X.shape#注意：这里行为样本数，列为维mean_X=X.mean(axis=0)#注意：axis=0表示计算每列（维）均值，结果为行向量X=X-mean_X#M=np.dot(X.T,X)#把X转置后代入，得到M=np.dot(X,X.T)#跟书上一样e,EV=np.linalg.eigh(M)#求出特征值和特征向量#tmp=np.dot(X,EV).T#把X转置后代入，得到tmp=np.dot(X.T,EV).T#跟书上一样V=tmp[::-1]#以下是对V归一化处理，先省略，下面讲所以我们看到，其实算法的本质是一样的，只不过要注意的地方是维数和样本数反过来了，另外，对X的运算换成X的转置即可。类推的，如果X使用我们的上面的组织方式，调用pca_book函数的代码为V,EV,immean=pca_book(X.T)
归一化算法不同。因为使用书上的方法，在对特征值求平方根(np.sqrt(e))的时候会产生一个错误（负数不能开平方根），所以我这里使用的归一化方法是从《GuidetofacerecognitionwithPython》抄来的。
书上的pca算法多了一个判断分支，当dim<=num_data即维数少于样本数的时候直接使用SVD（奇异值分解）算法，显然在一般的人脸识别的例子里，不会被用到，因为单张92×112图像打平后维数为10304，而样本数为400，远远低于维数。
归一化原先以为归一化是一种比较简单的运算，一了解才发现原来是一种不简单的思想，在机器学习中常被使用，看回上面的代码：
S=sqrt(e)[::-1]#计算e的平方根再对结果倒序排列foriinrange(V.shape[1]):V[:,i]/=S我在网上找到了关于compattrick后如何对求得的向量归一化的数学推荐，截图如下：这跟左奇异值有关，属于SVD中的内容，有兴趣的话自行研究。当我使用这种方法实现时，程序运行出现错误：FloatingPointError:invalidvalueencounteredinsqrt，发现是对负数开平方根产生了错误，也就是说对协方差矩阵求得的特征值中包含了负数。那么，如果要使用这种归一化方法的话，是否只要排除掉负的特征值及其对应的特征向量就可以了？
SVD（奇异值分解）我们的代码示例的PCA方法使用的是特征分解，线性代数中，特征分解（Eigendecomposition），又称谱分解（Spectraldecomposition）,是将矩阵分解为由其特征值和特征向量表示的矩阵之积的方法，但需要注意只有对可对角化矩阵才可以施以特征分解。而SVD（singularvaluedecomposition）能夠用于任意m乘n矩阵的分解，故SVD适用范围更广。但是，如果矩阵维数很大，如我们之前所举例10000维的时候，计算SVD也是很慢的，所以我们看到书上的例子增加了一个分支判断，当维度<样本数的时候，才使用SVD，否则使用compattrick方法的PCA。回顾一下上篇笔记举的PCA应用例子：把一张二维图像变换成一维，维数为2，对于这个例子，直接使用SVD是比较合适的，这样PCA函数将变得很简洁。
这里不会详细地讨论如何实现一个好的人脸识别算法，而是为了示例PCA的运用，所以只是简单的介绍一下。上面我们已经得出了400张人脸的特征脸（特征向量），首先第一个问题，我们得选择多少个特征向量作为主成分？特征值贡献率假设我们选择k个特征向量，其对应的特征值之和与所有特征值之和的比值，就是这k个特征值的贡献率。所以主成分的选择问题就转化为选择k个特征向量，使得特征值的贡献率大于等于某个值（如90%）即可。我们把选定的k个特征向量组成的矩阵设为W。
识别步骤第一步：将每个人脸样本图像减去平均图像后，投影到主成分上
W=EV[:,k]#假设k已经根据特征值贡献率算出来了projections=[]#存放每个人脸样本的投影forxiinX.T:#X为我们之前组织的所有人脸样本的10000×400矩阵projections.append(np.dot(W.T,xi-mean_X))#mean_X为之前我们求得的平均图像第二步：设要识别的图像为D，将D也投影到主成分上得到Q，然后计算Q与各个样本人脸投影的欧几里得距离，得出最小的欧几里得距离
defeuclidean_distance(p,q):#求欧几里得距离p=np.array(p).flatten()q=np.array(q).flatten()returnnp.sqrt(np.sum(np.power((p-q),2)))minDist=np.finfo('float').maxQ=np.dot(W.T,D-mean_X)foriinxrange(len(projections)):dist=euclidean_distance(projections[i],Q)ifdist 人脸识别的方法有很多种，基于特征脸的识别只是其中一种。但要实现一个可用的人脸识别，还有很多问题要考虑。另外PCA本身对某些特定情况的原始数据集也存在一些缺点。至此，关于PCA，将不再深入探讨。在PCA的学习过程中，深感一种技术应用的背后，必有惊艳的数学原理支撑，体会了一把数学之美:)，但因本人数学水平有限（后悔大学时没好好学），对以上理解必存在错漏和不详之处，所以也是很欢迎你的批评指正。 Linux爱好者，技术积累主要在Linux、Qt、Android，后以Android开发为主，从上层（kotlin,java）到底层(jni,linux)有一定的工作经验和理解，擅长快速学习和知识关联梳理，整合不同技术资源为客户提供合适的解决方案。