非线性角度2DPCA及其在人脸识别中的应用

乔 慧，周水生

西安电子科技大学 数学与统计学院，西安710126

主成分分析（PCA）[1]是常见的降维方法并且广泛应用于人脸识别领域。PCA 的目的是寻找方差最大的低维数据。但在人脸识别问题中，PCA首先需要将二维图像矩阵逐行或逐列排列成高维向量，这在一定程度上破坏了图像的二维结构以及图像的邻域结构信息。并且当样本数目相对较小且图像分辨率高时，重新排列导致训练样本维数过高，协方差矩阵过大，求解困难。为解决这一问题，Yang等[2]提出了二维主成分分析（2DPCA），此方法不需要将矩阵重新排列成向量，而是直接将图像矩阵并排，这样大大降低了图像协方差矩阵的规模，同时也较好地保留了图像的二维信息。2DPCA在图像处理领域具有广泛的应用[3-4]。

PCA和2DPCA都基于最小化L2范数平方之和，但在出现离群点时L2范数平方会使得实际投影方向偏离实际值，所以其鲁棒性较差。为提升鲁棒性，L1范数被用来代替L2范数，Ke和Kanade[5]提出L1-PCA并广泛地应用[6]，通过最小化L1范数来得到投影矩阵。在此基础上，Wang 等[7]又提出了2DPCA-L1，能够很好地利用图像的空间结构。但L1范数的一个缺点是它不具有旋转不变性[8]，而旋转不变性有助于避免算法的性能下降。基于这一问题，Wang 等[9]提出的F2DPCA，最小化重建误差的F范数，使得算法具有旋转不变性。但基于F范数的情形下，最小化重构误差不等价于方差的最大化。为了克服这一缺点，Gao 等[10]提出了角度2DPCA 方法，由于其本质是极小化角度的正切值，本文称之为Tan-2DPCA，该方法使用F 范数作为距离度量，最小化重构误差和方差之间的关系，不仅具有较强的鲁棒性，而且还具有旋转不变性。受到角度概念的启发，Zhou等[11]提出了Sin-2DPCA，将重构误差与样本之间的关系集成到目标函数中，从而更好地解释了角度的概念。实验结果表明，Sin-2DPCA 与Tan-2DPCA 具有相近的性能，都具有较强的抗噪性。

虽然上述的主成分分析方法都可以有效降低数据的维数，同时具有一定的抗噪性能，但它们都是提取图像的线性数据特征，而人脸图像中往往含有非线性特征。针对非线性问题，Schölkopf 等[12]和Zheng等[13]提出了KPCA，其主要思想是通过非线性映射将样本数据从输入空间映射到Hilbert空间，然后在该特征空间中应用PCA。不可避免地，KPCA还需要将二维图像矩阵转换成一维向量，因此Nhat等[14]提出了一种基于核的2DPCA方法（K2DPCA），该方法直接基于输入图像矩阵提取非线性主成分，其原理是将图像矩阵并排放置形成一个较大矩阵，将此矩阵的行或列作为样本并且非线性地映射到希尔伯特空间，然后应用PCA 方法。KPCA 和K2DPCA在非线性空间具有良好的性能，但仍存在以下两个缺点。一是鲁棒性不够强，KPCA基本保持了传统PCA 噪声的敏感性[15]。为缓解KPCA 中噪声敏感性问题，Xiao等[16]提出基于L1范数的KPCA优化问题，将L1范数引入KPCA。二是由于K2DPCA 中核矩阵的大小是图像样本的数量与行数或列数的乘积，因此其计算复杂度和存储复杂度显著增加。为降低复杂度，Sun 等[17]和Eftekhari等[18]对核矩阵进行了划分，将大的核矩阵分成若干小矩阵，通过计算多个小矩阵的特征向量得到投影矩阵。但求若干个小矩阵特征向量的整个过程仍然很复杂，无法完全解决大存储复杂度的问题。尤其是当训练样本较大时，小的核矩阵的数量增加，特征向量的计算耗时也很大。Wang等[19]研究了图像矩阵的矢量和矩阵混合表示方法，以降低核矩阵的大小。然而，这些方法可以减小核矩阵的大小，但在某种程度上破坏了图像的二维空间结构，当训练样本很大时，核矩阵的规模巨大，将核矩阵分成若干小矩阵的数量很大，这些方法并不能有效地降低核矩阵的规模。

本文针对KPCA 的鲁棒性以及算法的计算和存储复杂度作了研究。首先，结合角度2DPCA[10-11]鲁棒的优势，将角度的概念推广到非线性空间，提出Sin-K2DPCA。其次，针对Sin-K2DPCA计算复杂、存储复杂，提出基于Cholesky 分解的Chol+SinK2DPCA 以降低计算和存储复杂度。最后，给出实验结果及分析，并做出总结。

1 相关算法

本章简要介绍相关的主成分分析算法，包括PCA[1]、KPCA[12]、2DPCA、K2DPCA[14]和角度2DPCA[10-11]。

1.1 PCA与KPCA

假设有m个样本X=[x1,x2,…,xm]，xi∈ℝn，将样本中心化后（仍记为xi），PCA 和KPCA 都是通过求解式（1）中的模型，从而得到投影矩阵W=[w1,w2,…,wp]，其中wi是具有最大方差的正交特征向量：

在KPCA中，St=K，其中K为核矩阵并满足Kij=k(xi,xj)=φ(xi)Tφ(xj)，φ为从ℝn→ℝf的非线性映射，计算K的前p个最大特征值对应的特征向量得到投影矩阵W，从而任一样本x降维得到y=WTKx，其中Kx=[k(x1,x),k(x2,x),…,k(xm,x)]T。

1.2 2DPCA与K2DPCA

设m个训练样本图像[A1,A2,…,Am]，Ai∈ℝs×t。2DPCA 通过求解式（2）中的模型得到投影矩阵R=[r1,r2,…,rp]，其中r1,r2,…,rp是正交的：

在2DPCA中，式（2）的最优解R通过计算协方差矩阵的前p个最大特征值对应的特征向量得到，从而对样本Ai降维可得Yi=AiR,Yi∈ℝs×p。令Ei=Ai-AiRRT为第i个训练样本的重构误差，则。因此，模型（2）等价于式（3）：

在K2DPCA 中，首先将m个训练样本按行排列为A=[A1,A2,…,Am]∈ℝs×mt，其次，将通过映射ψ：ℝs×t→ℝs×f投影到高维希尔伯特空间ℝs×f得到，进而得到核矩阵K，最后求K的前p个最大特征值对应的特征向量得到投影矩阵（其本质上相当于对A的每一行或列执行KPCA），从而对样本Ai降维得Yi=AiR,Yi∈ℝs×p。

1.3 角度2DPCA

2DPCA采用平方F范数作为距离度量，因此对异常值敏感，当存在异常值时导致解的偏差。为了提高2DPCA的鲁棒性和保持旋转不变性，通过最小化F范数度量下的重构误差与方差之间的比，提出了一个角度2DPCA。它的优化模型如下：

称这个模型为Tan-2DPCA，在式（4）中目标函数是样本与其方差夹角的正切值。

随后，基于相对重构误差的几何解释，通过最小化重构误差与样本的比值，提出了Sin-2DPCA[7]，模型如下：

对于问题（4）和（5），没有封闭解，文献[10]中采用迭代算法求解。当重构误差与样本之间的角度α较小时，有sinα≈tanα，所以Sin-2DPCA算法性能与Tan-2DPCA相似，较2DPCA、L1-2DPCA鲁棒性强，但算法收敛性仍是待解决的问题。Sin-2DPCA增强角度2DPCA的几何解释，降低了计算复杂度，节省了训练时间。与2DPCA相比，Tan-2DPCA 和Sin-2DPCA 具有两个优点。一是在F 范数下度量重构误差，当存在异常值时，实际特征偏离较小且具有旋转不变性，因此它提高了鲁棒性。二是角度2DPCA 的目标函数同时考虑重建误差和方差，从而捕获更有效的人脸图像信息进而提高识别率。但Sin-2DPCA优化了角度的几何解释，最小化重构误差与样本的比值，比Tan-2DPCA的计算更简便。

2 非线性角度主成分分析

基于Sin-2DPCA在线性空间中良好的鲁棒性，本文将此推广到非线性情形并对计算复杂度作出研究。由于非线性二维主成分分析的本质是将图像样本按行或列排列后形成大的矩阵，再以其行或列作为样本执行KPCA。因此，首先研究非线性一维角度主成分分析（Sin-KPCA），再推导得到非线性二维角度主成分分析（Sin-K2DPCA），最后研究Sin-K2DPCA 的计算复杂问题，提出基于Cholesky 分解[15]的Sin-K2DPCA（Chol+SinK2DPCA）。

2.1 Sin-KPCA

为了将“角度”概念引入非线性空间，首先研究线性角度PCA（Sin-PCA），其模型为：

根据F范数与矩阵的迹的关系，有如下等价关系：

据文献[10]中定理2，将式（8）中的G(Wt)TWt奇异值分解为QΣRT，其中Σ为对角阵，从而得到最优解是Wt+1=QRT。

非线性情形下，将任一向量xi通过非线性映射φ：ℝn→ℝf，映射到高维特征空间ℝf，引入核矩阵K为Kij=k(xi,xj)=φ(xi)Tφ(xi) ，则非线性角度PCA，即Sin-KPCA模型为：

算法1给出Sin-KPCA的求解步骤。

算法1 Sin-KPCA

输入：{x1,x2,…,xm},xi∈ℝn，误差上界ε；

初始化投影矩阵Ut0，满足(Ut0)TUt0=Ip,t0=0。

1. 计算核矩阵Kij=k(xi,xj)=φ(xi)Tφ(xj),i,j=1,2,…,m并中心化K；

3. SVD分解G(Ut)Ut=QΛRT得到Ut+1=QRT；

4. 若‖Ut-Ut-1‖F ＞ε,则令t←t+1，返回2；

5. 否则停止迭代，输出投影矩阵U；

6. 将U带入计算Wφ,则任一样本x将通过(Wφ)Tφ(x)=UKx,Kx=[k(x1,x),k(x2,x),…,k(xm,x)]投影到低维特征空间。

2.2 Sin-K2DPCA

Sin-KPCA 算法进行人脸识别时，需要将二维图像矩阵重新排列成向量，破坏了图像的二维结构，根据文献[14]提出的K2DPCA算法，K2DPCA的本质是在所有图像样本按行或列排列成大的矩阵后，以大矩阵的列或行作为训练样本来执行KPCA。因此，Sin-K2DPCA 同样可经过Sin-KPCA 得到，具体模型为式（11）。与Sin-KPCA不同的是，式（11）中的K是将图像按行或列排列后得到的核矩阵，U是通过求解新的关于核矩阵K的极小化问题对应的解。

给定训练样本图像A=[A1,A2,…,Am] ，将Ai∈ℝs×t通过非线性映射ψ:ℝs×t→ℝs×f映射到高维空间ℝs×f中，令为矩阵A的第j行，j=1,2,…,s,则：

其中，φ是与2.1 节中相同的映射。令i=1,2,…,m,j=1,2,…,s,l=s(i-1)+j是第i个样本图像矩阵的j行的转置。从而得到核矩阵K∈ℝms×ms由Kij=k(zi,zj)=φ(zi)Tφ(zi),i,j=1,2,…,ms构成。进一步，以K的每一行或列作为样本执行Sin-KPCA。

与Sin-KPCA 不同的是，Sin-K2DPCA 中核矩阵维数是训练样本的数量与样本矩阵行或列的乘积，即ms×ms维或mt×mt维的。因此，当训练样本规模较大时，Sin-K2DPCA 在实际计算时核矩阵的规模会成倍增加，计算复杂且存储困难，甚至会超过计算机负荷。为解决此问题，下面提出Cholesky分解的方法避免计算整个核矩阵，减少计算量与存储空间。

2.3 基于Cholesky分解的Sin-K2DPCA

为缓解Sin-K2DPCA 算法的计算复杂度和存储复杂度较大这一问题，本节基于Cholesky 分解的算法，将核矩阵低秩分解以便减少计算复杂度和存储复杂度。通过使用Cholesky 分解方法，得到大规模核矩阵K的低秩近似。在整个低秩近似过程中，只需计算核矩阵的对角线元素和部分精选列，并仅需存储低秩矩阵L，这降低了计算复杂度并节省了存储空间。

当有m个s×t维的训练样本时，核矩阵是ms×ms或mt×mt维的，直接计算整个核矩阵计算复杂度过大。为避免这个问题，采用选主元的Cholesky方法迭代确定基本集B⊂M:={1,2,…,m} ，从而得到近似矩阵T。迭代过程中，记第i次的迭代误差为Ei=K-，将Ei对角线上最大元素的列号加入到基本集Bi中，从而得到Bi+1，进一步构造K的低秩近似矩阵,详细见文献[20]。这样最大程度地极小化trace(Ei)，减小近似误差，得到的近似矩阵是在迹范数意义下最优的低秩近似。

通过选主元的Cholesky 分解方法，可以将Sin-K2DPCA中的核矩阵近似分解为K≈LLT,其中rank(K)=r(r≪m),L∈ℝm×r。对近似核矩阵中心化为Kc=K-，其中1 ∈ℝM×M是分量全部为1 的矩阵，是以P的行均值向量为列排成的矩阵。在实际计算中，无需直接计算大规模的核矩阵，也不需要对K进行SVD分解。记,可通过PTP来得到投影矩阵。事实上，将核矩阵低秩分解后等价于将Sin-K2DPCA转化到线性空间中，用线性的办法求解，这样将减小了计算和存储复杂度。

本节提出的基于选主元的Cholesky 分解的Sin-K2DPCA 算法，记为Chol+SinK2DPCA，具体如算法2所示。

该算法进行Cholesky更新时计算复杂度为Ο(msr)，进而奇异值分解时的计算复杂度为Ο(msr2)，故算法的总的计算复杂度为Ο(msr)+Ο(msr2)。与Sin-K2DPCA相比，此方法避免了计算整个核矩阵，只需计算低秩矩阵P，从而减小了计算量和存储空间。因此Chol+SinK2DPCA 不仅提升了算法在核空间的鲁棒性，而且降低了计算复杂度，能够有效实现大规模数据的降维。

算法2 Chol+SinK2DPCA

输入：{A1,A2,…,Am},Ai∈ℝs×t,误差上界ε,低秩上界r。

初始化投影矩阵Ut0，满足(Ut0)TUt0=Ip,t0=0；

1. Cholesky分解核矩阵得到低秩近似并中心化K≈PPT；

3. SVD分解G(Ut)Ut=QΛRT得到Ut+1=QRT；

5. 否则停止迭代，输出投影矩阵U。

3 实验结果及分析

本章通过实验对KPCA、K2DPCA 和提出的方法Sin-K2DPCA、Chol+SinK2DPCA进行比较，所有算法都用MATLAB2017a 实现，并在CPUi7-4790、3.60 GHz 和8 GB RAM 的计算机上运行。为更好地发挥线性核与高斯核的优点，本实验采用线性与高斯相结合的核函数。对于所有人脸识别实验，使用基于欧氏距离的1-NN算法来分类。

3.1 ORL数据库

本实验比较了KPCA、K2DPCA、Sin-K2DPCA 和Chol+SinK2DPCA 在ORL 标准人脸数据库中的性能。参数选择不是本文研究的重点，因此核参数的确定采用简单的网格方法，即为不同的方法选择集合{2-12,2-11,…,21}中性能最好的参数。对该数据库，KPCA 采用γ=2-6，所有非线性二维主成分分析采用γ=2-1。

ORL数据库包含40个人，每人有10幅不同的图像，部分图像是在不同时间拍摄的，包括不同的灯光、面部表情（睁/闭眼、微笑/不笑）和面部细节（眼镜/不戴眼镜），在本实验中所有灰度图像分辨率为48×48。随机选取10%的图像，将其用矩形遮挡，形成方块噪声。噪声是灰度噪声图像，像素在0到255之间随机生成，噪声大小为图像大小的5%～15%，加噪后的部分图像如图1所示。随机选取50%的图像进行训练，其余的进行测试，即每次实验训练图像200张，测试图像200张。

图1 ORL数据库的原始图像与加噪图像

考虑到Cholesky分解算法中参数r的选取，一般不需要直到满足trace(K-LLT)≤ε终止，通常采用核矩阵低秩逼近的秩r终止分解方法。通常情形下，参数r=20。对KPCA、K2DPCA、Sin-K2DPCA和Chol+SinK2DPCA，实验分析了主成分个数p的变化对测试识别精度的影响以及最好识别精度的训练时间，实验结果如图2 和表1所示。

图2 加噪的ORL数据集上测试识别率

表1 ORL数据集上测试识别率比较

图2表示在加噪声的ORL数据集中，测试识别率随主成分个数的变化情况。由图2 可以看出，K2DPCA、Sin-K2DPCA与Chol+SinK2DPCA算法明显优于KPCA算法。随着主成分个数p的增加，K2DPCA、Sin-K2DPCA与Chol+SinK2DPCA的识别率均逐渐增大，当主成分个数p=25 时趋于稳定。而KPCA 的识别率还有上升的趋势，下文的实验中将KPCA的主成分个数取更大的值p=100，从而得到KPCA 的最高识别精度。通过对比，可以看出角度的主成分分析方法有较好的抗噪效果，并且Cholesky 分解的Sin-K2DPCA 方法的识别精度是最高的。可能是因为在选主元的过程中，增强了抗噪性。

在表1中，列出由KPCA、K2DPCA、Sin-K2DPCA和Chol+SinK2DPCA 算法在原始（无噪声）数据集和噪声数据集上的平均识别精度（实验10次的平均值）以及均方差，给出了原始数据集和噪声数据集上识别精度的差异，最后一列给出了算法在取得最好识别精度时的训练时间。

由表1 可得，在原始数据集上，这三种算法的差异可以忽略不计，而在噪声数据集上，Sin-K2DPCA 与基于Cholesky 分解的Sin-K2DPCA 方法明显优于KPCA和K2DPCA。这是因为角度K2DPCA 模型最小化了相对重构误差，并且采用F 范数作为距离度量，这不仅对异常值具有较强的鲁棒性，还保留了旋转不变性。此外，第三列数值表明，所有算法的性能都受到噪声的影响，而Chol+SinK2DPCA方法是所有比较算法中最好的一种，具有最高的精度和最小的减小量。由最后一列知，在识别精度提高的情况下，Chol+SinK2DPCA 的训练时间与Sin-K2DPCA 相比缩短了近一半。由于将核矩阵进行选主元的Cholesky 分解成K≈PPT，在进行SVD求解时，避免了计算整个核矩阵，只需对分解的矩阵P做运算即可，此时的计算复杂度为Ο(9 600×202)，这样大大减小了计算复杂度，从而缩短了训练时间。

3.2 Yale数据库

为了进一步研究所提出算法的性能，在Yale数据集进行了实验，该数据集包含了15 个个体的GIF 格式的165张灰度图像。每人有11个图像，且有不同的面部表情或光照。实验中，每张灰度图像的分辨率为48×48。本实验中，KPCA 的核参数为γ=2-6，所有K2DPCA 核参数为γ=2-3。

将样本随机添加如实验3.1 中加入随机矩阵噪声，示例如图3。随机选取每人的6 张图像作为训练样本，剩余的5张作为测试，即每次实验有90个训练图像和75个测试图像。实验结果如图4和表2。

图3 Yale数据库原始图像和加噪图像

图4 加噪声的Yale数据集的测试精度

表2 加噪声的Yale数据集的测试精度比较

表2 为四种方法在Yale 数据库上的识别率和训练时间的比较（所有结果都是10次随机测试的平均值）。

图4 再次支持了前面的实验结果，在此数据集中，Chol+SinK2DPCA的性能要比其他算法的好。由表2可知，在抗噪性较好的同时，Chol+SinK2DPCA 的计算速度为Sin-K2DPCA的一半。

3.3 Extended YaleB数据库

最后选取Extended YaleB 大型人脸数据库，说明Chol+SinK2DPCA 方法对大规模数据的有效性。本实验中，KPCA 的核参数为γ=2-6，所有K2DPCA 核参数为γ=2-3。

Extended YaleB 数据库由2 414 张正面图像组成，这些图像是从38 个不同光照的个体中采样的，每人有64 幅图像。在该数据库中，每张图像被调整为48×42像素。在实验中随机选择每人5张图像，添加方形遮挡的，噪声的位置是随机的，噪声大小与图像大小的比值为0.05～0.15，如图5。在实验中，随机选择每个人的25 幅图像进行训练，其余图像进行测试。此时，训练样本总数为950（38×25），KPCA需要计算核矩阵的维数是950×950，K2DPCA 需要计算的核矩阵维数45 600×45 600，这超出了实验所用电脑的计算和存储能力，但采用Cholesky 分解的Sin-K2DPCA 方法可以有效实现。本实验仍然采用线性与高斯结合的核函数。

图5 Extended YaleB原始图像和加噪图像

由图6 可以看出，在加噪声和干净数据上，Chol+SinK2DPCA的识别精度相差1.22%，这表明该方法具有良好的抗噪性能。同时，由于不必计算整个核矩阵，而是将其分解减小了计算机的存储空间以及计算复杂度，克服了原始K2DPCA算法不能识别大规模样本的缺点。

图6 测试识别率随主成分个数的变化

以上实验结果充分体现了Sin-K2DPCA 与Chol+SinK2DPCA 的优越性。主要原因是，在角度K2DPCA中，基于F 范数最小化重构误差，使得在异常值对结果的影响变小，在ORL 与Yale 数据库中，角度2DPCA 方法的鲁棒性明显提升。此外，对大规模的核矩阵选主元的Cholesky 分解，不需计算整个核矩阵，只需计算与存储低秩近似的矩阵，大大减小了计算复杂度和存储复杂度，在Extend YaleB数据集中解决了原始方法由于核矩阵过大而不能运行的问题。

4 结束语

本文首先提出了一种新的K2DPCA算法，将核空间中的相对重构误差降到最小，提高了K2DPCA的鲁棒性。在此基础上，提出了Chol+SinK2DPCA算法。该算法降低了计算复杂度，对孤立点的鲁棒性优于Sin-K2DPCA和K2DPCA。对于大型数据集，Chol+SinK2DPCA 克服了由于计算量大K2DPCA 不能实现的问题。同时对部分人脸数据库的实验结果表明，该算法对小规模训练集具有较强的鲁棒性。