PCA 人脸识别

基于主成分分析（PCA）人脸识别的分析与实现。首先讨论了矩阵乘法的含义，接着尝试推导了主成分的计算过程，并最终将其应用于人脸图像，在多维空间中采用欧氏距离比较人脸数据。

1. 简介

PCA，即主成分分析。简单来说就是通过坐标系变换将位于原坐标系中的原数据变为在新坐标系下对应的新数据。通常这种变化伴随着数据维度的大幅降低，但是数据信息的损失较小。即减少不同维度之间的相关系，去除信息冗余，提取主成分的意思。

对于人脸图片识别来说，图片包含的信息量是十分巨大的。一张256*256分辨率的图像将每行像素收尾相接就会形成一个长度为65536大小的一维数组。换句话说一张256*256分辨率的图像拥有65536个维度。因此对于更高分辨率的图片更大的样本数量，如果不对数据降维则对如此巨大的数据进行处理在时间上就已经是不可取的。

2. 矩阵乘法的几何含义

先看下列式子：

$$\begin{align*} Y &= PX \\ \left[ \begin{matrix} c _1 \\ c_2 \\ \vdots \\ c_n \end{matrix} \right] &= \left[ \begin{matrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} \\ \end{matrix} \right] \times \left[ \begin{matrix} b _1 \\ b_2 \\ \vdots \\ b_n \end{matrix} \right] \\ \\ \left[ \begin{matrix} 2 \\ 0 \end{matrix} \right] &= \left[ \begin{matrix} 1 & 1 \\ -1 & 1 \\ \end{matrix} \right] \times \left[ \begin{matrix} 1 \\ 1 \\ \end{matrix} \right] \end{align*} \tag1$$

通常来讲一个矩阵乘以一个向量表示表示对这个向量伸缩旋转变换。例如上式所示向量 $X=[1,1]^T$ 经过矩阵相乘后就转变成向量 $Y=[2,0]^T$ 了。但是当矩阵$P$ 的每一行都是线性无关的向量时，那么矩阵$P$ 的所有横向量就可以作为一个新的坐标系的基向量。用它乘以一个其他向量就相当于将该向量的坐标转换为新坐标系下的新坐标，如图所示：

由于新基矩阵（黄色）维向量（横向量）的模长为$\sqrt2$ ,所以新坐标系中的向量的模长也要比原向量长$\sqrt2$ 倍。

3. PCA降维的数学运算过程

通常样本表示的列向量，如下所示表示一个含有m维数据的样本：

$$x=\lbrack x_1, x_2, \dots, x_m \rbrack^T \tag2$$

则原数据$X$可以表示为如下：（$m \times n$ 矩阵， $n$ 个样本，每个样本有$m$ 维数据）

$$X= \left[ \begin{matrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} \\ \end{matrix} \right] \tag3$$

原数据的协方差$C$ 可以表示为如下：

$$C = \frac{1}{m}XX^T \tag4$$

式$(2)、(3)、(4)$ 为已知初始数据，下面将给出我们的目的数据，设$Y$ 来表示降维后的数据，如下所示：

$$Y= \left[ \begin{matrix} b_{11} & b_{12} & \dots & b_{1n} \\ b_{21} & b_{22} & \dots & b_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ b_{k1} & b_{k2} & \dots & b_{kn} \\ \end{matrix} \right] (k<m) \tag5$$

降维数据$Y$ 的协方差可以表示为如下：

$$D = \frac{1}{m}YY^T \tag6$$

到此，我们希望降维数据$Y$ 的协方差矩阵$D$ 上主对角线上的元素尽可能大，而其他元素都为零。因为协防差矩阵表示的是不同维度的相关性，所以这样希望的话就意味着数据在同一维度上尽可能的分散并且维度与维度之间没有相关性。因此想要将原数据$X$ 变为新数据$Y$ 则需要一个新的基矩阵$P$ ，使得：

$$Y=PX \tag7$$

式$(7)$ 即为原数据与降维数据的联系，我们最终的目的就是找到满足$D$ 的基矩阵$P$ 。观察下述等式：

$$\begin{align*} D &= \frac{1}{m}YY^T \\ &= \frac{1}{m}(PX)(PX)^T \\ &= \frac{1}{m}PXX^TP^T \\ &= P(\frac{1}{m}XX^T)P^T \\ &= PCP^T \end{align*} \tag8$$

由于原数据$X$ 的协方差矩阵$C$ 为实对称矩阵，因此实对称矩阵$C$ 通过特征分解后如下所示：

$$\begin{align*} C &= E\Lambda E^T \\ &= [e_1, e_2, \dots, e_m] \left[\begin{matrix} \lambda_1 & 0 & \cdots & 0 \\ 0 & \lambda_2 & \cdots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & \cdots & \lambda_m \end{matrix}\right] \tag9 [e_1, e_2, \dots, e_m]^T \\ & (其中Ce_i=\lambda_ie_i, e_i是特征值\lambda_i所对应的特征向量，列向量) \end{align*}$$

在式$(9)$ 中我们可以看到矩阵$\Lambda$ 是非常符合降维数据$Y$ 的协方差矩阵$D$ 的要求，主对角线上的元素尽可能大，而其他元素都为零。因此我们联立式$(8)、(9)$ 得：

$$\begin{equation} \begin{cases} D = PCP^T\\ \Lambda = E^TCE \end{cases} \tag{10} \end{equation}$$

由式${10}$ 可得:

$$P=E^T \tag{11}$$

由式$(11)$ 可知，我们所求的基矩阵即为原数据$X$ 的协方差矩阵$C$ 经过特征值分解后的由特征向量组成的矩阵$E^T$ 。在矩阵$E^T$ 中，每一行均是一个特征向量，也是新坐标系中的基向量。在$E^T$ 中有$m$ 行（原数据有$m$ 维）。根据$\Lambda$ 中特征值得大小从$E^T$ 中删去特征值较小的维度保留特征值较大的维度从而构成新的基矩阵$P$ ,因此在$P$ 应当少于$m$ 行，减少维度从而达到降维的效果。

$$\begin{align*} Y &= PX \\ &= \left[\begin{matrix} c_{11} & c_{12} & \dots & c_{1m} \\ c_{21} & c_{22} & \dots & c_{2m} \\ \vdots & \vdots & \ddots & \vdots \\ c_{k1} & c_{k2} & \dots & c_{km} \\ \end{matrix} \right]_{k \times m} \times \left[ \begin{matrix} a_{11} & a_{12} & \dots & a_{1n} \\ a_{21} & a_{22} & \dots & a_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \dots & a_{mn} \\ \end{matrix} \right]_{m \times k} \\ & = \left[ \begin{matrix} d_{11} & d_{12} & \dots & d_{1n} \\ d_{21} & d_{22} & \dots & d_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ d_{k1} & d_{k2} & \dots & d_{kn} \\ \end{matrix} \right]_{k \times n} \\ \end{align*} \tag{12}$$

由式$(12)$ 可以看到，降维数据$Y$ 由原数据$X$ 的$m$ 维降低至$k$ 维$(k<m)$ 。

4. 计算特征值与特征向量

通常来讲图像的维数是非常高的，比如一张 256×256 分辨率的图像拥有 65536 个维度。而图像样本的数量通常则远远维度数。所以对于原数据$X$ 的协方差矩阵$C$ 是非常大的。例如 20 张 256×256 的样本，那么$X$ 的大小为 65536×20，其协方差矩阵的大小为 65536×65536。

对于矩阵相乘我们已知：

$$r(AB)\leqslant min \{ r(A),r(B) \}$$

因此协方差矩阵$C$ 的秩 $r(C) \leqslant n(n为样本数)$ 一遍来讲$r(C) = n$ 。由于$C$ 是实对称矩阵且相似于对角矩阵，故$C$ 的秩就等于其非零特征值的个数，也就意味值当样本数量远远小于维数时，其协方差矩阵的特征值大多数是为 0 的。

因此考虑矩阵$H$ ：

$$H=\frac{1}{m}X^TX \tag{13}$$

$H$ 的大小为$n×n(n为样本数)$ ，其大小远远小于$m×m(m为维数)$ 。现在对$H$ 进行特征分解，其特征值和特征向量可以表示为： $H\alpha_i=\lambda_i\alpha_i$ 。考虑下列等式：

$$\begin{align*} H\alpha_i &= \lambda_i\alpha_i \\ X^TX\alpha_i &= \lambda_i\alpha_i \\ XX^T(X\alpha_i) &= \lambda_i(X\alpha_i) \\ C(X\alpha_i) &= \lambda_i(X\alpha_i) \end{align*} \tag{14}$$

由式$(14)$ 可知，若 $\lambda_i$ 和 $\alpha_i$ 是矩阵$H$ 的特征值和特征向量，那么若 $\lambda_i$ 和 $X\alpha_i$ 也是协方差矩阵$C$ 的非零特征值和其对应的特征向量。这样就可以求出$C$ 所有的非零特征值和其对应的向量。在 PCA 中就可以根据这些非零特征值进行筛选主成分了。通过矩阵$H$ 可以极大的减少所需的计算量。

5. 人脸识别过程

到目前为止，我们就可以将人脸图像的样本数据进行降维处理得到由特征脸$b_i$ 组成的特征集$Y$ 。

$$\begin{align*} Y &= [b_1, b_2, \dots, b_n]\\ \\ &= \left[ \begin{matrix} b_{11} & b_{12} & \dots & b_{1n} \\ b_{21} & b_{22} & \dots & b_{2n} \\ \vdots & \vdots & \ddots & \vdots \\ b_{k1} & b_{k2} & \dots & b_{kn} \\ \end{matrix} \right] \tag{15} \end{align*}$$

对与一张新的图像，可以表示成一个$1×n$ 大小的向量$v$ 。首先将其投影到降维空间也就是用矩阵$P$ 去乘以向量$v$ ,从而得到向量$u$ 。

$$u=P \times v \tag{16}$$

再通过分别计算向量$u$ 当每一个特征脸$bi$ 的欧氏距离，取欧氏距离最短的特征脸$b\{min}$ 所对应的样本即为本次所识别的结果。

6. 代码

源代码参考自互联网$^{[1]}$,增加了些许调整和注释并使用automake进行构建。当前程序识别精确度有限，若使用分类器代替简单的欧氏距离判断则可提高识别精确度。

代码编译：

$> cd build && ../src/configure && make

代码运行：

$> ./pcafacerecog -a <存放训练图像的目录> -e <存放测试图像的目录>

7. 参考资料

1. PCA人脸识别学习及C语言实现，http://blog.csdn.net/jinshengtao/article/details/18599165
2. 基于PCA的人脸识别步骤，http://blog.csdn.net/yutianzuijin/article/details/10823985
3. PCA的数学原理，http://blog.codinglabs.org/articles/pca-tutorial.html
4. 浅谈协方差矩阵，https://www.cnblogs.com/chaosimple/p/3182157.html