本文比较通俗易懂，希望能把人脸识别这回事讲得清楚明了

本文第一部分先通俗地从人的角度讲清楚人脸识别是怎么一回事，以及一些专有名词的含义
本文第二部分讲一下算法及其数学原理
本文第三部分通过用python实践最简单不带任何优化的PCA算法以及PCA-LDA算法，实践非常简单，仅作入门了解机器学习算法

最后，本文不会出现太多公式，主要在代码之间体现数学运算过程

阅读本文第一部分只需要： 一颗好奇的心
阅读本文第二部分只需要：

1. 线性代数的基础
2. 了解概率论中“协方差”的概念

阅读本文第三部分只需要：

1. 有一定的面向对象编程的思想
2. 对python基础语法有一定的了解
3. 有一双会查python文档的手

1.人脸识别这回事

1.1我们是怎样认出一张脸的

我们先想想我们一般是怎样认出一张脸的，我们认出一张脸的第一步是，看他，看她。当然，作为一个脸盲症晚期患者，我是没办法做到“只是在茫茫人群中看了你一眼，就再也无法忘记你容颜 ”，所以需要看好多遍才认清一张脸的。

深层次地讲，我们为什么看一张脸很多次，就可以认出一个人的呢？很简单，我们要通过“学习”来认出一张脸，作为一个莫得感情的机器，就是来比对每一张脸的像素值是否接近。但是毕竟都大数据时代了，一张图片少说也有几百万像素，这计算量刷个脸还要几十秒可不太行，机器总要智能点吧。再认真想想，我们认清一张脸，真的会比对一张脸的每一个部分吗？显然不是，哆啦A梦的嘴巴特别大，汤姆猫的下巴特别尖，hello kitty眼睛特别小，我们一下子就把这三只猫给识别出来了，那么，我们凭借的是什么呢？特征！那，机器能不能也靠特征人脸识别呢？

1.2机器该如何认出一张脸

对于一张图片，想想它有什么特征呢？显然，每一个像素点就是一个"特征"，例如，对于一张200*200的图片，它有40000个特征，那么，对于机器而言，它要选择哪几个像素点来识别会更好，没有听错，就是“哪几个”。

所以，我们想象自己的大脑是一台机器，人脸识别的过程就是，第一步，先输入大量的图片，先认清楚这张脸长啥样，在机器学习中，这些数据被称为训练集；第二步，我们通过这么多的图片，来“学习”，就会发现，比较眼睛，鼻子，嘴巴……可以更容易分辨出这个人，在机器学习中，这就是模型的训练 ，举个例子，训练集输入了“哆啦A梦”，“汤姆猫”，“hello kitty”三只猫，然后机器发现了，通过比较嘴巴……可以将其分开，于是就选择了合适的特征了；第三步，机器已经通过“学习”知道了比较这些部位可以识别出一张脸，所以就需要我们输入一些图片来检测它的“学习”效果，上述方法被称为特征脸（eigenface）它的数学原理是主成分分析，这种方法单纯比较每张照片差别较大的特征，简称PCA。

假如我把情况变一下，只需要把“哆啦A梦”，“汤姆猫”，“hello kitty”三只猫分辨出来，那就肯定是，抓住他们不同类之间的特征，他们耳朵上的差别不大 ，所以一般地人脸识别的时候，需要会看五官，而在这种情况下，是不需要看耳朵的。我们仔细想想这个对比过程，假设我们分别有10张不同表情下的哆啦A梦，汤姆猫，hello kitty三只猫的图片，我们通过比较这这三只猫的图片，发现他们区别最大的特征，同时，我们还要确保这个特征在同一只猫的不同表情下都一样，比如，我发现哆啦A梦的嘴巴很大，汤姆猫正常情况下则不是很大，于是，机器就把比这个特征用来决策了……但是汤姆猫在见到杰克的时候，嘴巴张的比哆啦A梦还大，就有可能把汤姆归为哆啦A梦了，还是会造成误差的。为了解决这个问题，我们选取的“好特征”还需要类内差别比较小才行。综上，特征选择的另一个方法就是，使类内区别尽可能小，类间区别尽可能大。这种方法就是线性判别分析，简称LDA。

但是对比一下，想想这两种方法的优劣，前者明显适用性更广，也就是泛化能力较强（通过学习这30张哆啦A梦，汤姆猫，hello kitty的图片，为了学习到一般识别人脸需要的特征），而后者，明显是对于特定的数据集的拟合更好（通过学习这30张哆啦A梦，汤姆猫，hello kitty的图片，只为了把它们区分），所以有可能出现过拟合现象（“读死书”现象），因此后者更适用于有标签的有监督学习，前者适用于无监督或者半监督学习。

过拟合现象即是，由于一个模型训练的样本太少或者是算法问题，导致过于“死板”，举个例子，假如你告诉一个小孩，榕树叶是树叶，小孩子由于他没怎么见过树叶，就误认为树叶一定是无锯齿的，一定是绿色的，所以当他看到枫树叶是红色的时候，他就会认为枫树叶并不是树叶了。

过拟合：
只知道榕树叶是树叶–>树叶是是绿色的
遇到枫树叶—>枫树叶是红色的—>枫树叶不是树叶

当然，以上说法可能有失偏颇，毕竟通过LDA和PCA投影之后，得到的特征的意义是模糊的，也就是说，机器只能分辨出特征的重要性，至于利用的时候，就仅仅是为了保留信息的意义模糊的数据了。

1.3我们该如何“教”机器认出一张脸

好了，现在大概思路已经清楚下来了，我们先来把一大堆图片（训练集）输入给机器，然后机器分析这一大堆图片，他学会了怎么比较图片，然后他就开始用这种方法对付别的图片（也就是测试集），接下来的问题就是，如何通过数学方法“教会”机器。

2.理论:LDA与PCA这回事

2.1数学基础

对于本文中的数学方法，不同于考试，做题，侧重于介绍其含义而非其计算与证明，例如，做矩阵乘法的时候，我们要关心的是，这个矩阵的规模及每个元素的含义，至于值是多少、如何运算，还是留给计算机解决这个问题吧。

2.1.1数据预处理

这里只是简单介绍一下，若要深入了解，可以看这篇文章,我觉得作者讲得挺好的

标准化/归一化：把有量纲量变为无量纲量

1、min-max标准化（离差标准化）：

缺点：当有新数据加入时，可能会导致max和min的变化，需要重新定义

2、Z-score标准化（0-1标准化）：

2.1.2 利用投影来实现降维

先来直观理解，降维这回事
这是位于二维空间的一些数据点

把每个数据点平移到原点附近，即用每个数据点减去其中点，至于为什么要这么做，后面算方差的时候就清楚了，并进行“归一化”，所谓“归一化”即是把数据点的各个维度缩放到[-1,1]的范围中，归一化方法很多，在此不赘述

发现数据点之间大致分布，选取合适的方向，将数据点投影下去，即可

最终数据成功地从二维降到了一维

关于如何投影，假设读者有高中数学的基础，可以按高中数学来理解：
问题：有一数据点$(4,5)$在$e_1=(1,0),e_2=(0,1)$的平面中现想要将其投影到$e_1^{\prime }=(\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2}),e_2^{\prime }=(-\frac{\sqrt{2}}{2},\frac{\sqrt{2}}{2})$中
解：
设坐标为$$(x,y)$$
根据高中知识，有：$$x\ast \frac{\sqrt{2}}{2}+y\ast (-\frac{\sqrt{2}}{2})=4$$
$$x\ast \frac{\sqrt{2}}{2}+y\ast \frac{\sqrt{2}}{2}=5$$
则写成矩阵形式有：$$\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& -\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{c}x\\ y\end{array}\right]=\left[\begin{array}{c}4\\ 5\end{array}\right]$$
最后由于规范正交矩阵矩阵的逆与其转置相等，故有：
$$\left[\begin{array}{c}x\\ y\end{array}\right]={\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& -\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]}^T\left[\begin{array}{c}4\\ 5\end{array}\right]$$
故可将向量左乘某个矩阵的转置使其投影到其列空间

2.1.3 寻找最佳降维方向

从上面的例子可以看出，在投影的时候，关键在于投影到多少维，往哪个方向投影。我们可以再观察一下上面的图， 显然，投影方向是最分散的方向。试想，假如投影方向选择了“密集”的方向，那么带来的后果将是，严重的信息损失。
举个例子，假如有100个数据点，很不幸选择了一个“坏”的特征，即一个“坏”的投影方向，这些点被投影之后，绝大部分重合了，甚至有些靠的太近以至于无法区分，只剩下10个，那么这100个数据点所剩余的信息可能只剩于10个数据点所含的信息了。
那么，在数学上我们是利用什么来衡量数据的“分散程度”呢？
我们选用“方差”来衡量分散程度。

因此，我们只需要使所有投影后的数据的方差最大化即可
另外，对于投影而言，我们投影是为了用尽量少的特征描述全部的信息，假如这些特征有很强的相关性，那么我们有可能可以用其中一个取代另一个，所以我们还需要将数据维度之间的相关性减小，为了衡量不同维度之间的相关性，我们引用了“协方差”的概念

2.1.4协方差与协方差矩阵

此外，将协方差公式中的b改成a则为方差的公式，至于为何分母为m-1而非m，涉及到数理统计中的无偏估计，不在本文的讨论范围内

了解了协方差与方差公式后，我们要处理数据需用矩阵，这时，我们便需要协方差矩阵。协方差矩阵的计算公式为:$$A^T\ast A$$举个例子

观察易得:

2.1.5特征值，特征向量以及对角化最大方差

观察上面的协方差矩阵，就会发现，协方差矩阵的对角线元素即为方差，故$trace(C)$即为方差之和

特征值与特征向量:
n阶矩阵A的n个特征值是其特征方程$$|E-A|=0$$的n个根$${\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n$$而属于特征值${\lambda }_i$的特征向量就是线性方程$$({\lambda }_{i}E-A)x=0$$的非零解。
$$Ax=\lambda x$$表明x在经过矩阵A的线性变换后方向不变，只是伸缩了λ倍。

将特征向量看做空间的基向量，则矩阵A就是这些基向量进行伸缩变换时使用的数学运算。给定一个矩阵A，就可以找出对应的基（特征向量），以及通过线性变换（矩阵A），对这些基进行的伸缩变换（特征值）。

我们的目的是为了最大化方差，最小化协方差，故可以将矩阵对角化即可最小化协方差以及最大化方差（当然这样不太严谨，可用下面的瑞利商问题求解或是用拉格朗日乘子法求解）

这时我们选取特征值排序前k大的特征向量，即可组成向量空间，将数据投影到其中即可

此外，协方差矩阵是实对称矩阵，他有以下性质：
实对称矩阵不同特征值对应的特征向量必然正交。
实对称矩阵一定可以对角化

2.1.6瑞利商问题以及广义瑞利商问题（这一部分不懂可跳过）

先来介绍一下埃尔米特(Hermitan)矩阵




瑞利商是指这样的函数
$$R(A,x)/R(A,x)=(x^{H}Ax)/(x^{H}x)$$
其中x为非零向量，而A为n×n的Hermitan矩阵。，即A^{H} =A。如果我们的矩阵A是实矩阵，则满足A^T=A的矩阵即为Hermitan矩阵。
瑞利商具有如下重要性质：
$${\lambda }_{min}\le (x^H{A}_x)/(x^{H}x)\le {\lambda }_{max}$$
它的最大值等于矩阵A最大的特征值，而最小值等于矩阵A的最小的特征值。

当向量x是标准正交基时，即满足$$x^{H}x=1$$时，瑞利商退化为： $$R(A,x)=x^{H}Ax$$

广义瑞利商是指这样的函数$$R(A,B,x)/R(A,B,x)=(x^{H}Ax)/(x^{H}Bx)$$
其中x为非零向量，而A,B为n×n的Hermitan矩阵，并且B为正定矩阵。
通过标准化可以将广义瑞利商转化为瑞利商的形式。令$$x=B^{(}-1/2)x^{\prime }$$则

而分子转化为：

此时$$R(A,B,x)$$转化为$$R(A,B,x^{\prime })$$:

此时转化为求$B^{-1/2}A{B}^{-1/2}$特征值和特征向量的问题

2.1.7奇异值分解

先来看看概念介绍（也可跳过直接看例子）：

对于矩阵A，求出特征值和特征向量后，有$$AW=\Lambda W$$其中$W$是$A$的特征向量组成的矩阵，$\Lambda$是$A$的特征值组成的对角矩阵。即$A$可以分解成$$W\Lambda W^{-1}$$再将$W$进行标准化和正交化处理后，$W$满足$$W^{T}W=I$$则A可以写为$$W\Lambda W^T$$
但以上情况要求A是n×n的方阵，而当A为m×n,m≠n时，如果想要对其进行矩阵分解，就要用到SVD。

SVD实际上是在A的列空间、零空间、行空间和左零空间中找到了两组正交基，使得列空间、零空间中的正交基V经过A变换后等于行空间、左零空间中的正交基U与特征值Λ的乘积，即$$AV=U\Lambda$$ $$A=U\Lambda V^T$$

m×n的矩阵A可以分解为$$A=U\Lambda V^T$$其中U为m×m矩阵，Λ为m×n对角阵，V为n×n矩阵。U和V均为酉矩阵，即满足$$U^{T}U=I，V^{T}V=I$$
由$$A^{T}A=V\Lambda U^{T}U\Lambda V^T=V{\Lambda }^2{V}^T$$解$A^{T}A$的特征值和特征向量可得$\Lambda$和$V^T$，再由$$A{A}^T=U\Lambda V^{T}V\Lambda U^T=U{\Lambda }^2{U}^T$$解$A{A}^T$的特征值和特征向量可得$U$。即可得到$$A=U\Lambda V^T$$

通俗一点理解奇异值分解：跟非方阵矩阵的对角化差不多，因为我们对方阵对角化的时候，是令det(A-λI)=0,而对于非方阵det必定为0，故无解。

还是来看个栗子吧
$$A=\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]$$
$$A^{T}A=\left[\begin{array}{ccc}0& 1& 1\\ 1& 1& 0\end{array}\right]\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]$$
对其对角化：
$${\lambda }_1=3,{\lambda }_2=1$$
$$u_1=\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\end{array}\right],u_2=\left[\begin{array}{c}-\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\end{array}\right]$$
$$A{A}^T=\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]\left[\begin{array}{ccc}0& 1& 1\\ 1& 1& 0\end{array}\right]=\left[\begin{array}{ccc}1& 1& 0\\ 1& 2& 1\\ 0& 1& 1\end{array}\right]$$
对其对角化：

$${\lambda }_1=3,{\lambda }_2=1,{\lambda }_3=0$$
$$u_1=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right],u_2=\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]u_3=\left[\begin{array}{c}\frac{\sqrt{3}}{3}\\ -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{3}}{3}\end{array}\right]$$

观察一下$A{A}^T$，其实就是上面讲的协方差矩阵，再观察一下求解过程，非常巧，他们的非零特征值是一样的，而在PCA中，显然，特征值是0的是太小了，不要的。另外，有兴趣的话，可以根据秩与非零特征值的关系证明这并不是巧合

$$A=U\Lambda V^T=\left[\begin{array}{ccc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{3}& 0& -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\\ 0& 0\end{array}\right]\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\\ -\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]$$

然后我们试着删去为0的特征值，及其对应的特征向量，我们采用分块的思想进行运算：
(未删去)
$$A=U\Lambda =\left[\begin{array}{ccc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{3}& 0& -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\\ 0& 0\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}0& 1\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{3}}{3}\\ -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}0& 0\end{array}\right]$$
假如把它删去:
$$A=U\Lambda =\left[\begin{array}{cc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}\\ \frac{\sqrt{6}}{3}& 0\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}0& 1\end{array}\right]$$
观察不难发现，二者是完全等价的，删去一个为0的特征值，相当于删去一个零矩阵，对于原矩阵并没有损失，因此，SVD可用于对数据的压缩。具体有何作用，之后再讲。

2.2算法思路

令人头秃的数学基础终于结束了，现在是轻松愉快的算法实现了

2.2.1 PCA-Eigenface算法大体思路

1. 读入人脸，并进行预处理（归一化和中心化）。
2. 将人脸拉成一列或者一行，以一行为例，240张46*56图片中，假如这样处理，便可以组成一个240行2576列的图片，每一行的2576个元素依次为2576个特征，例如三个特征可以理解为三维空间的一个点，同理，这个例子就可理解为240个点处于一个2576的空间中。
3. 对这个矩阵求协方差矩阵。
4. 直接利用刚刚数学原理的成立条件，将特征向量按特征值大小来选取前k个特征向量，组成特征空间，例如，选取200个特征向量，即可组成240*200的特征空间。
5. 将测试集和训练集同时投影到特征空间中，然后使用KNN分类器进行比对分类，测试准确率
   

2.2.2 PCA-LDA算法的大体思路

1. PCA的第1~4步
2. 此时，已经得到了PCA投影的特征空间D（按上面的例子，是240*200的矩阵），再使用LDA，将240张的10类样本，分别求出类内散度矩阵，类间散度矩阵，为了最大化类间散度矩阵，最小化类内散度矩阵利用瑞利商问题进行求解，经过不太复杂的数学推理(见上2.1.5)，最终可转化为求$B^{-1/2}A{B}^{-1/2}$特征值和特征向量的问题
3. 假设PCA算法选取了200维的特征，形成一个$240\ast 200$的矩阵，则再利用LDA选择50个特征值，可形成一个$240\ast 50$的投影矩阵
4. 将测试集和训练集同时投影到特征空间中，然后使用KNN分类器进行比对分类，测试准确率
   

3. 实战:Python这回事

3.1python下的实现(这一部分采用无微不至的讲解)

由于博主是先写好PCA代码，再添加LDA代码，故有些函数与变量命名可能有些奇怪

直接介绍PCA-LDA的实现，毕竟，把它删去一部分就可以变成LDA或者PCA了

3.1.1.介绍一下引入的模板

import math							##为了使用向下取整函数
import cv2							##为了读入并处理图片
import numpy						##以矩阵的形式处理数据
import operator						##用于使用operator.itemgetter
from sklearn import preprocessing	##emmm归一化的函数懒得写，所以，emmm直接调包

3.1.2.由于python定义常量不是很方便，所以直接定义变量也行

(博主当年刚开始写代码的时候，就是因为代码中出现了太多幻数，结果修改的时候出了bug，活生生熬到差不多凌晨四点

SIZE=(46,56)			##图片的规模
Tr_num_i=40				##数据集中类的数量
Tr_num_j=6				##数据集中每一类用于训练的数量
sum_num_j=10			##数据集中每一类的总数量
dimension_PCA=500		##第一步PCA特征选择的维数
dimension_LDA=10		##第二步LDA特征选择的维数

3.1.3读入人脸图片:

def create(number_1,number_2,start):image_matrix=[]for i in range(1,number_1+1):for j in range(start,number_2+1):image=cv2.imread('orl'+str(i)+'_'+str(j)+'.bmp',cv2.IMREAD_GRAYSCALE)       ##输入图片，i相同为同一类image=cv2.resize(image,SIZE)                                                ##调整大小image_array=image.tolist()                                                  ##图片转换成一维矩阵temp=numpy.array(image_array)image_array2=temp.flatten()image_matrix.append(image_array2)                                           ##把100张图片存起来在一个二维数组中return image_matrix

可以设置断点看看结果（规模为240*2576）

3.1.4对训练集中心化和归一化(标准化)

中心化函数

def centre(Matrix,num):mean_array = numpy.mean(Matrix, axis=0)  		##求均值mean_array=numpy.tile(mean_array,(num,1))       ##把向量复制成矩阵diff_matrix=Matrix-mean_array                   ##中心化return diff_matrix,mean_array					##返回中心化之后的矩阵和平均矩阵，返回平均矩阵只是在好奇平均脸长啥样

平均脸:

归一化(标准化)

Tra_Matrix=numpy.array(Tra_Matrix)				##将数据转为矩阵
Tra_Matrix = preprocessing.scale(Tra_Matrix)	##将数据标准化

看看归一化之后的结果

3.1.5 生成PCA投影空间

3.1.5.1.中心化并求协方差矩阵

def Train_eigenface(Matrix):diff_matrix, mean_array=centre(Matrix,Tr_num_i*Tr_num_j)	##中心化tran_matrix=numpy.transpose(diff_matrix)					##求矩阵的转置cov=numpy.matmul(tran_matrix,diff_matrix)                   ##求协方差矩阵

查看运行结果:
可以看出，$cov$是一个$2576\ast 2576$的矩阵

3.1.5.2对协方差矩阵对角化，并对特征向量和特征值排序

    eigenvalues,eigenvectors1=numpy.linalg.eig(cov)eigenvectors=list(eigenvectors1.T)index=numpy.argsort(eigenvalues)

运行查看结果:

值得注意的是index:

它是一个索引，尺寸为2576，如图所示 index[2321]=279 则说明，特征值排序为第2322(从0开始计数)的特征向量的下标为279

3.1.5.3根据排序结果进行特征选择，并形成PCA的特征空间

    for i in reversed(index):								##由于之前排序是由小到大，故将index倒序if(i<dimension_PCA):eigenvectors2.append(eigenvectors[i])			##一个一个地加入特征eigenfaces=numpy.matmul(eigenvectors2,diff_matrix.T)	##将训练集投影过去，并生成特征脸return eigenvectors2,eigenfaces.T,mean_array			##返回特征空间和特征脸

运行查看结果:

根据调试结果，不难看出，原数据从$240\ast 2576$被降维成了$240\ast 500$维（若是PCA，不需要这么高的维数，但是要使用LDA二次选择，因此，选择较高维数）

3.1.6生成LDA投影空间

3.1.6.1将$240\ast 500$的矩阵分割成40个$6\ast 500$的矩阵

def LDA(eigenfaces):mean_class_b1={}											##建立一个空的dicttemp=numpy.array(numpy.split(eigenfaces,Tr_num_i,axis=0))	##按8行为一个单位分割for i in range(0,Tr_num_i):if(i!=0):mean_class_b1=numpy.vstack((mean_class_b1,numpy.mean(temp[i],axis=0)))mean_class_b1=numpy.array(mean_class_b1)else:mean_class_b1=numpy.mean(temp[0],axis=0)mean_class_b1=numpy.array(mean_class_b1)
##把这个40*6*500的矩阵求平均值，并复制成240*500的矩阵

运行结果如下:

3.1.6.2求解$S_b$

其实就是直接把公式翻译成代码

mean_class_b1=numpy.array(mean_class_b1)						##转换类型diff_matrix1,mean_array1=centre(mean_class_b1,Tr_num_i)		##中心化sb=numpy.matmul(diff_matrix1.T,diff_matrix1)				##求Sb即类间散度矩阵

运行结果:

3.1.6.3求解$S_w$

    sw=0for i in range(0,Tr_num_i):diff_matrix2,mean_array2=centre(temp[i],Tr_num_j)	##中心化diff_matrix2=numpy.array(diff_matrix2)				sw_temp=numpy.matmul(diff_matrix2.T,diff_matrix2)	##计算每一个类的类内散度矩阵sw+=sw_temp											##求和sw=sw/40												##求均值

运行结果:

3.1.6.4求解$S_b^{-1}{S}_w$的特征值，并按3.1.5的方法排序，进行选择10维特征

    j_w=numpy.matmul(numpy.linalg.pinv(sw),sb)				##j_w即为Sb^(-1)Sweigenvalues,eigenvectors1=numpy.linalg.eig(j_w)			##对角化eigenvectors=list(eigenvectors1.T)index=numpy.argsort(eigenvalues)eigenvectors2=[]for i in reversed(index):if (i < dimension_LDA):eigenvectors2.append(eigenvectors[i])diff_matrix,mean_array=centre(eigenfaces.T,len(eigenfaces[0]))eigenfaces = numpy.matmul(eigenvectors2, diff_matrix)	##将PCA投影下的数据投影到LDA的投影空间中eigenvectors2=numpy.array(eigenvectors2)				return eigenfaces.T,eigenvectors2						##返回投影后的数据和特征空间

结果如下:
eigenvectors为$10\ast 500$

而返回的eigenface2为$240\ast 10$的矩阵，表示240张人脸的10个特征

3.1.7以同样的方式读入测试集，并将其投影到特征空间中

def test_martrix(Matrix,eigenvectors1,eigenvecotos2):diff_matrix1,mean_matrix1=centre(Matrix,Tr_num_i*(sum_num_j-Tr_num_j))diff_matrix1=numpy.array(diff_matrix1)test_eigenface1=numpy.matmul(diff_matrix1,eigenvectors1.T)			##投影到PCA空间中diff_matrix2, mean_matrix2 = centre(test_eigenface1, (Tr_num_i*(sum_num_j-Tr_num_j)))diff_matrix2 = numpy.array(diff_matrix2)test_eigenface=numpy.matmul(diff_matrix2,eigenvecotos2.T)			##投影到LDA-PCA空间中return test_eigenface

这个过程相当于$160\ast 2576$的矩阵右乘一个$2576\ast 500$的矩阵，得到一个$160\ast 500$的矩阵，表示$160$张人脸的$500$个特征。再由$160\ast 500$的矩阵右乘一个500*10的矩阵，得到一个$160\ast 10$的矩阵，表$160$张人脸的$500$个特征，表示$160$张人脸的$500$个特征
运行结果如下:

3.1.8 分类并求取准确率及召回率

我们通过以上步骤，已经将数据集（包括训练集和测试集）投影到特征空间中了，接下来，就需要输入一张测试集，假如图片类别属于A它能准确分类，即为准确

3.1.8.1 KNN分类器原理

原理：利用已经标签好的数据集，输入没有标签的新数据，将新数据的每个特征与样本集中数据所对应的特征进行比较，然后选取前K个最相似的数据进行投票，得票数最高的标签即为新数据的分类。

高维数据可以通过PCA算法降维到低维空间，然后利用两点间的距离公式计算距离，将得到的结果作为是否为相似的判别依据。

如图所示，红绿蓝三颜色是已经做好标记的数据，对新输入的黄色点数据可以利用K近邻算法对其进行预测。很明显当K取3时，距离黄色点最近的三点标签均为中等生，因此新输入数据标签也被标为中等生。
但若K取值不同时，对新数据的预测可能会发生变化。

3.1.8.2 KNN分类器的代码实现

输入并求出测试样本（1个）与训练集中所有的样本(240个)的距离之差

def knn(eigenfaces,test_eigenface,k):row=len(eigenfaces)									##获得训练集的行数，即人脸数量diff=numpy.tile(test_eigenface,(row,1))-eigenfaces	##将测试样本按行复制row次，使之变成同规格的矩阵，随后作差sqr_diff=diff**2									##求距离的平方distance=sqr_diff.sum(axis=1)**0.05					##开方

运行结果如下:

对距离进行排序，并投票

 sort_distance=distance.argsort()count={}for i in range(0,k):									##投票K次vote=math.floor((sort_distance[i])/Tr_num_j)+1	##由于训练集命名的情况，需向下取整count[vote]=count.get(vote,0)+1         			##原本应该是count[vote]++,这里相当于建立一个键值对，键是vote，值一直++sort_count=sorted(count.items(),key=operator.itemgetter(1),reverse=True)return sort_count[0][0]

运行结果：
sort_distance为
例如，sort_distance[0]=0就说明了，与测试样本最接近的训练样本即为第一个，k=5即说明了要取前k个，向下取整再加1是因为数据集是从1开始计数的而非0

得到了vote，再对vote计数得到count，如图所示：

count为1：5，说明了最近的5个类别均是1类，故把该样本分为1类

3.1.8.3 测试

测试即是用KNN分类器分出的类别与测试样本本身标签比较，倘若是相同，则准确，否则，为错误

def test(eigenfaces,test_eigenface,k):summ=0.0row=len(test_eigenface)for i in range(0,row):if(knn(eigenfaces,test_eigenface[i],k)==math.floor(i/(sum_num_j-Tr_num_j))+1):summ+=1.0return summ/((sum_num_j-Tr_num_j)*Tr_num_i)

结果如图所示

在选取10维特征的情况下，有142个分类准确的测试样本，故准确率为142/160

3.2实验结果

一个简单的结果统计

可见，对于LDA较PCA充分利用了标签的信息，故识别率高于传统的PCA

彩蛋：SVD优化

显然没有什么彩蛋，这只是个标题党
先来复习一下之前讲过的SVD原理
奇异值分解
先来看看概念介绍（也可跳过直接看例子）：

对于矩阵A，求出特征值和特征向量后，有$$AW=\Lambda W$$其中$W$是$A$的特征向量组成的矩阵，$\Lambda$是$A$的特征值组成的对角矩阵。即$A$可以分解成$$W\Lambda W^{-1}$$再将$W$进行标准化和正交化处理后，$W$满足$$W^{T}W=I$$则A可以写为$$W\Lambda W^T$$
但以上情况要求A是n×n的方阵，而当A为m×n,m≠n时，如果想要对其进行矩阵分解，就要用到SVD。

SVD实际上是在A的列空间、零空间、行空间和左零空间中找到了两组正交基，使得列空间、零空间中的正交基V经过A变换后等于行空间、左零空间中的正交基U与特征值Λ的乘积，即$$AV=U\Lambda$$ $$A=U\Lambda V^T$$

m×n的矩阵A可以分解为$$A=U\Lambda V^T$$其中U为m×m矩阵，Λ为m×n对角阵，V为n×n矩阵。U和V均为酉矩阵，即满足$$U^{T}U=I，V^{T}V=I$$
由$$A^{T}A=V\Lambda U^{T}U\Lambda V^T=V{\Lambda }^2{V}^T$$解$A^{T}A$的特征值和特征向量可得$\Lambda$和$V^T$，再由$$A{A}^T=U\Lambda V^{T}V\Lambda U^T=U{\Lambda }^2{U}^T$$解$A{A}^T$的特征值和特征向量可得$U$。即可得到$$A=U\Lambda V^T$$

通俗一点理解奇异值分解：跟非方阵矩阵的对角化差不多，因为我们对方阵对角化的时候，是令det(A-λI)=0,而对于非方阵det必定为0，故无解。

还是来看个栗子吧
$$A=\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]$$
$$A^{T}A=\left[\begin{array}{ccc}0& 1& 1\\ 1& 1& 0\end{array}\right]\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]=\left[\begin{array}{cc}2& 1\\ 1& 2\end{array}\right]$$
对其对角化：
$${\lambda }_1=3,{\lambda }_2=1$$
$$u_1=\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\end{array}\right],u_2=\left[\begin{array}{c}-\frac{\sqrt{2}}{2}\\ \frac{\sqrt{2}}{2}\end{array}\right]$$
$$A{A}^T=\left[\begin{array}{cc}0& 1\\ 1& 1\\ 1& 0\end{array}\right]\left[\begin{array}{ccc}0& 1& 1\\ 1& 1& 0\end{array}\right]=\left[\begin{array}{ccc}1& 1& 0\\ 1& 2& 1\\ 0& 1& 1\end{array}\right]$$
对其对角化：

$${\lambda }_1=3,{\lambda }_2=1,{\lambda }_3=0$$
$$u_1=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right],u_2=\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]u_3=\left[\begin{array}{c}\frac{\sqrt{3}}{3}\\ -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{3}}{3}\end{array}\right]$$

观察一下$A{A}^T$，其实就是上面讲的协方差矩阵，再观察一下求解过程，非常巧，他们的非零特征值是一样的，而在PCA中，显然，特征值是0的是太小了，不要的。另外，有兴趣的话，可以根据秩与非零特征值的关系证明这并不是巧合

$$A=U\Lambda V^T=\left[\begin{array}{ccc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{3}& 0& -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\\ 0& 0\end{array}\right]\left[\begin{array}{cc}\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\\ -\frac{\sqrt{2}}{2}& \frac{\sqrt{2}}{2}\end{array}\right]$$

然后我们试着删去为0的特征值，及其对应的特征向量，我们采用分块的思想进行运算：
(未删去)
$$A=U\Lambda =\left[\begin{array}{ccc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{3}& 0& -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}& \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\\ 0& 0\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}0& 1\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{3}}{3}\\ -\frac{\sqrt{3}}{3}\\ \frac{\sqrt{3}}{3}\end{array}\right]\left[\begin{array}{cc}0& 0\end{array}\right]$$
假如把它删去:
$$A=U\Lambda =\left[\begin{array}{cc}\frac{\sqrt{6}}{6}& \frac{\sqrt{2}}{2}\\ \frac{\sqrt{6}}{3}& 0\\ \frac{\sqrt{6}}{6}& -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\\ 0& 1\end{array}\right]=\left[\begin{array}{c}\frac{\sqrt{6}}{6}\\ \frac{\sqrt{6}}{3}\\ \frac{\sqrt{6}}{6}\end{array}\right]\left[\begin{array}{cc}\sqrt{3}& 0\end{array}\right]+\left[\begin{array}{c}\frac{\sqrt{2}}{2}\\ 0\\ -\frac{\sqrt{2}}{2}\end{array}\right]\left[\begin{array}{cc}0& 1\end{array}\right]$$
观察不难发现，二者是完全等价的，删去一个为0的特征值，相当于删去一个零矩阵，对于原矩阵并没有损失，因此，SVD可用于对数据的压缩。

在PCA的降维中，我们使用了协方差矩阵，但是协方差矩阵是$Cov=A^{T}A$
但是在 SVD中我们发现一个性质就是，$A{A}^T$与$A^{T}A$的非零特征值是相同的。
首先 $$A^{T}Ax=0$$与$$Ax=0$$同解，故他们的秩是一样的
同理有
$$rank(A)=rank(A^T)=rank(A{A}^T)=rank(A^{T}A)$$
设$$A^{T}Ax=\lambda x(\lambda 为对角矩阵)$$
则有$$A{A}^{T}Ax=A\lambda x$$
又因为$\lambda$是对角矩阵
则有$$A{A}^T(Ax)=\lambda (Ax)$$
故$A{A}^T$与$A^{T}A$的非零特征值是相同的。

那么，这有什么用呢？
回到上面的例子，有2576个特征，也就是说，对角化$ATA$的时候，我们计算出了2576个特征值，但是令计算机失望的就是，这2576个特征值里，只有240个是有用的。

抽象地讲，就是向量张成的空间小于其所在的空间，而PCA-LDA的本质就是最小化向量张成的空间，并使用它来取代向量所在的空间。举个例子，有2个点在一个三维空间中，在不改变原点的情况下，这两个点就是在一个平面内，这就说明了，两个三维向量张成的空间为二维空间，从矩阵的形式上看，假如这两个向量是线性无关的，则秩为2。
回到人脸识别的例子假如说有240个2576维的数据点，那么，这些数据点所在的空间即为2576维的空间。而这些数据点张成的空间为240维，也就是说，假如以一个240维的空间描述这些数据点，也是完全无碍的，是可以描述全部信息的。

本质地说，PCA使用的是右奇异矩阵对数据的列进行降维，而左奇异矩阵则是对数据的行进行降维

因此，使用SVD进行降维可大大优化计算量

部分图源网络，侵删
若有写的不清楚的地方，可在评论区讨论
（由于本人水平有限，疏忽错误在所难免，还请各位指正）