引言

近年来，随着科学技术的发展，数据的数量和维度呈现出快速增长的趋势。降维不仅可以简化模型、降低模型的训练时间，而且可以提高模型的泛化能力。因为高维数据难以表达且需要较高的计算和存储代价，所以降维方法成为机器学习和模式识别领域中的热门研究方向，并已在人脸识别、手写体识别等系统中取得了良好的效果。

降维方法可分为线性和非线性两种。线性降维方法，如主成分分析法(Principal component analysis, PCA)^[1]、线性判别分析(Linear discriminant analysis, LDA) ^[2]等，几何表达简单、计算成本较低，是分析高维数据的基础。PCA由Karl Pearson于1901年提出, 该方法使用正交线性变换将一组观测的可能相关变量转换为一组线性不相关变量，这些不相关变量叫做主成分。拥有最大方差，即尽可能不相关的数据称为第一主成分，并以此类推。PCA算法的核心是对数据协方差矩阵进行特征分解。LDA算法与方差分析和回归分析密切相关，它试图找到一组特征的线性组合，能将多个对象尽量分开，即找到最佳投影方向。

如果所需的特征在一个低维的流形子空间中，那么传统的线性降维算法并不能很好地提取出该非线性结构，因此很多非线性的流形学习算法被提出，包括等距映射(Isomap)^[3]、拉普拉斯特征映射(Laplacian eigenmaps, LE)^[4]及局部线性嵌入(Locally-linear embedding, LLE)^[5]等。

但是以上流形学习算法只考虑了训练样本，处理过程中未能得到一个显式映射，因此不能处理新样本，即所谓的“out of sample”问题。针对这个问题，He等人提出了保局算法(Locality preserving projection, LPP)^[6]和局部保持嵌入算法(Neighborhood preserving embedding, NPE) ^[7]:(1) LPP算法首先整合邻近点信息建立一个图，该图可看做一个连续流形结构的线性离散近似^[6]，然后使用拉普拉斯算子计算出转换矩阵将数据点映射到低维子空间; LPP算法使原始高维空间中距离相近点在降维后依旧保持近距离。(2) NPE算法的目的在于保持局部流形结构。在一定范围内给定一些数据点，首先建立权重矩阵来描述数据点之间的关系，并将每个点用邻近点线性表示; 然后找到一个最佳嵌入使得该邻近结构在低维空间中得以保持。以上两个利用局部信息进行降维的算法虽然没有“out of sample”问题，但非局部信息在此过程中被忽略，即算法没有考虑相隔较远的点之间的关系。2007年，Yang等人提出了非监督判别投影(Unsupervised discriminant projection，UDP)算法^[8]。该方法考虑了样本点在低维表示后的局部和非局部性，最小化局部散度的同时最大化非局部散度，即使非局部散度与局部散度的比值尽可能大。但UDP是非监督的学习方法，也就是不能利用样本的类别信息，不利于分类；而且在处理过程中也出现小样本问题^[9]：当训练样本个数远小于样本维数时，就会出现矩阵奇异问题。2012年，楼宋江等人提出的小样本有监督拉普拉斯判别分析算法(Supervised Laplacian discriminant analysis, SLDA)^[9], 通过将类内拉普拉斯散度矩阵投影在去除零空间的总体拉普拉斯散度矩阵上，避免了小样本问题。但是由于该方法中没有明确最大化类间散度，因而导致在投影空间的可分性不够好。2015年，丁春涛等人提出的基于双邻接图的判别近邻嵌入算法(Double adjacency graphs-based discriminant neighborhood embedding, DAG-DNE)^[10], 融入了构建双邻接图的想法，解决了判别近邻嵌入(Discriminant neighborhood embedding, DNE)图嵌入降维方法在构造权重矩阵时数据不平衡的问题。

受到DAG-DNE算法的启发，本文提出了一种新的基于拉普拉斯判别分析的有监督降维方法，称为适用于小样本的双邻接图判别分析方法(Double adjacent graph-based discriminant analysis for small sample size, DAG-DA)。DAG-DA在构建邻接图时考虑了样本的不平衡问题，每个样本点都至少建立两个关联，并将局部拉普拉斯散度矩阵与非局部拉普拉斯散度矩阵的差投影到总体拉普拉斯散度矩阵的秩空间中，然后在该空间中进行特征问题的求解，求得投影矩阵。通过在简单人工数据集、Yale人脸库和ORE人脸库的实验证明了DAG-DA算法的优越性。

1 适用于小样本的SLDA算法

针对小样本问题，已有相应的一些解决方案^{[9, 11-12]}，其中监督化拉普拉斯判别分析算法和本文的工作紧密相关。该算法是一种有监督的流形学习算法，其目标是最小化局部散度，同时最大化非局部散度，即非局部散度与局部散度的比值尽可能大^[13]，表达式为

$ \max J\left( \mathit{\boldsymbol{W}} \right) = \frac{{{J_{\rm{N}}}\left( \mathit{\boldsymbol{W}} \right)}}{{{J_{\rm{L}}}\left( \mathit{\boldsymbol{W}} \right)}} \Leftrightarrow \min J\left( \mathit{\boldsymbol{W}} \right) = \frac{{{J_{\rm{L}}}\left( \mathit{\boldsymbol{W}} \right)}}{{{J_{\rm{T}}}\left( \mathit{\boldsymbol{W}} \right)}} $

(1)

式中：W为投影矩阵，J_L为局部散度，J_N为非局部散度，J_T为总体散度，J_T=J_L+J_N。

以上优化问题(1)可以转化为如下目标函数

$ {\mathit{\boldsymbol{S}}_{\rm{L}}}\mathit{\boldsymbol{w = }}\lambda {\mathit{\boldsymbol{S}}_{\rm{T}}}\mathit{\boldsymbol{w}} $

(2)

式中：λ为特征分解得到的特征值，w为特征向量，S_L为局部拉普拉斯散度矩阵，S_T为丢弃零空间的总体拉普拉斯散度矩阵。

2 适用于小样本的DAG-DA算法 2.1 双邻接图判别分析

由于适用于小样本的监督化拉普拉斯判别分析算法在分子上仅考虑将局部拉普拉斯散度矩阵进行投影, 忽略了非局部拉普拉斯散度矩阵。为了强调非局部信息，本文考虑在最优化函数里引入非局部散度, 表达式为

$ \min J\left( \mathit{\boldsymbol{W}} \right) = \frac{{{J_{\rm{L}}}\left( \mathit{\boldsymbol{W}} \right) - {J_{\rm{N}}}\left( \mathit{\boldsymbol{W}} \right)}}{{{J_{\rm{T}}}\left( \mathit{\boldsymbol{W}} \right)}} = \frac{{{J_D}\left( \mathit{\boldsymbol{W}} \right)}}{{{J_{\rm{T}}}\left( \mathit{\boldsymbol{W}} \right)}} $

(3)

局部散度为

$ {J_{\rm{L}}} = \frac{1}{{2nn}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{F_{ij}}{{\left( {{y_i} - {y_j}} \right)}^2}} } $

(4)

非局部散度为

$ {J_{\rm{N}}} = \frac{1}{{2nn}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{B_{ij}}{{\left( {{y_i} - {y_j}} \right)}^2}} } $

(5)

式(4, 5)中：F为相似度矩阵；B为差异度矩阵。与SLDA方法从数据点的K近邻中找到同类和异类样本点来构建邻接图不同，本文在构造邻接图时借用双邻接图的思想，从同类和异类样本点间分别找到K个近邻点，再构建F和B。

设样本集{(x_i, y_i)}_i=1ⁿ，其中x_i∈ R^d，y_i∈{1, 2, …, c}为x_i的类别标签，n为样本的个数，d为样本的维数，c为类别总数。对于这个样本集，利用样本的类别信息，对每一个数据点都构造相似度矩阵F和差异度矩阵B，根据样本分布赋予不同的权值，表达式为

$ {F_{ij}} = \left\{ \begin{array}{l} \exp \left\{ { - \frac{{{\mathit{\boldsymbol{x}}_i} - \mathit{\boldsymbol{x}}_j^2}}{t}} \right\}\;\;\;{\mathit{\boldsymbol{x}}_i} \in N_k^ + \left( j \right)\;或\;{\mathit{\boldsymbol{x}}_j} \in N_k^ + \left( i \right)\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(6)

$ {B_{ij}} = \left\{ \begin{array}{l} - \exp \left\{ { - \frac{{{{\left\| {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right\|}^2}}}{\sigma }} \right\}\;\;\;{\mathit{\boldsymbol{x}}_i} \in N_k^ - \left( j \right)\;或\;{\mathit{\boldsymbol{x}}_j} \in N_k^ - \left( i \right)\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;其他 \end{array} \right. $

(7)

式中：N_k⁺(i)和N_k⁺(j)分别表示样本x_i和x_j的同类近邻集合，且x_i和x_j类别相同；N_k^-(i)和N_k^-(j)分别表示样本x_i和x_j的异类近邻集合，且x_i与x_j类别不相同；t和σ为可变参数。由式(6, 7)可以看出，当同类样本距离越近，它们的相似度越高，反之亦然。当异类样本距离越近，它们的差异度越小，反之亦然。

由式(4, 6)可以得到

$ {J_{\rm{L}}} = \frac{1}{{2nn}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{F_{ij}}{{\left( {{y_i} - {y_j}} \right)}^2}} } = \frac{1}{{nn}}{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{L}}}\mathit{\boldsymbol{W}} $

(8)

式中S_L为局部拉普拉斯散度矩阵

$ {\mathit{\boldsymbol{S}}_{\rm{L}}} = \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{F_{ij}}\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right){{\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right)}^{\rm{T}}}} } = \mathit{\boldsymbol{XF'}}{\mathit{\boldsymbol{X}}^{\rm{T}}} - \mathit{\boldsymbol{XF}}{\mathit{\boldsymbol{X}}^{\rm{T}}} = \mathit{\boldsymbol{XL}}{\mathit{\boldsymbol{X}}^{\rm{T}}} $

(9)

式中F′为一个对角矩阵，对角线上的元素为：$ F\prime {_{ii}} = \sum\limits_j {{F_{ij}}} $，L=F′－ F，L即为拉普拉斯类内矩阵。

同样，由式(5，7)可以得到

$ {J_{\rm{N}}} = \frac{1}{{2nn}}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{B_{ij}}{{\left( {{y_i} - {y_j}} \right)}^2}} } = \frac{1}{{nn}}{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{N}}}\mathit{\boldsymbol{W}} $

(10)

$ {\mathit{\boldsymbol{S}}_{\rm{N}}} = \frac{1}{2}\sum\limits_{i = 1}^n {\sum\limits_{j = 1}^n {{B_{ij}}\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right){{\left( {{\mathit{\boldsymbol{x}}_i} - {\mathit{\boldsymbol{x}}_j}} \right)}^{\rm{T}}}} } = \mathit{\boldsymbol{XB'}}{\mathit{\boldsymbol{X}}^{\rm{T}}} - \mathit{\boldsymbol{XB}}{\mathit{\boldsymbol{X}}^{\rm{T}}} = \mathit{\boldsymbol{XN}}{\mathit{\boldsymbol{X}}^{\rm{T}}} $

(11)

式中：B′为一个对角矩阵，对角线上的元素为$B\prime {_{ii}} = \sum\limits_j {{B_{ij}}} $；N=B′－B为拉普拉斯类间矩阵。

由式(3, 8, 10)可推导出目标函数等价于

$ \min \frac{{{J_{\rm{D}}}\left( \mathit{\boldsymbol{W}} \right)}}{{{J_{\rm{T}}}\left( \mathit{\boldsymbol{W}} \right)}} = \frac{{{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{D}}}\mathit{\boldsymbol{W}}}}{{{\mathit{\boldsymbol{W}}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{T}}}\mathit{\boldsymbol{W}}}} = \frac{{{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{D}}}\alpha }}{{{\mathit{\boldsymbol{\alpha }}^{\rm{T}}}{\mathit{\boldsymbol{S}}_{\rm{T}}}\mathit{\boldsymbol{\alpha }}}} \Leftrightarrow {\mathit{\boldsymbol{S}}_{\rm{D}}}\mathit{\boldsymbol{w = }}\lambda {\mathit{\boldsymbol{S}}_{\rm{T}}}\mathit{\boldsymbol{w}} $

(12)

式中：S_D =S_L－S_N，α∈range (S_T)。

2.2 秩空间

首先对S_T进行特征分解。由于总体拉普拉斯散度矩阵S_T非负定，所以对它进行特征值分解后的特征值可能为零，可能为正。取非零的特征值$ \mathit{\boldsymbol{ \boldsymbol{\tilde \varLambda} }}$所对应的特征向量U₁组成S_T的主空间α。即S_T特征分解后，特征值为

$ \mathit{\boldsymbol{ \boldsymbol{\varLambda} }} = {\rm{diag}}\left( {{\lambda _1},{\lambda _2}, \cdots ,{\lambda _r},{\lambda _{r + 1}}, \cdots ,{\lambda _n}} \right) $

其中特征值λ₁≥λ₂…≥λ_r，λ_r+1=λ_r+2=…=λ_n=0。特征值对应的特征向量矩阵为

$ \mathit{\boldsymbol{U}} = \left\{ {{\mathit{\boldsymbol{u}}_1},{\mathit{\boldsymbol{u}}_2}, \cdots ,{\mathit{\boldsymbol{u}}_r},{\mathit{\boldsymbol{u}}_{r + 1}}, \cdots ,{\mathit{\boldsymbol{u}}_n}} \right\} $

其中特征值λ_i对应的特征向量为u_i。

令非负特征值矩阵$\mathit{\boldsymbol{ \boldsymbol{\tilde \varLambda} }} $ =diag(λ₁, λ₂, …, λ_r), 部分特征向量矩阵U₁= {u₁, u₂, …, u_r}, U_⊥= {u_r+1, u_r+2, …, u_n }。令$ {\mathit{\boldsymbol{W}}_{\rm{1}}}\mathit{\boldsymbol{ = }}{\mathit{\boldsymbol{U}}_{\rm{1}}}\mathit{\boldsymbol{ \boldsymbol{\tilde \varLambda} }}$，将差异拉普拉斯散度矩阵S_D投影到W₁上为

$ {{\mathit{\boldsymbol{\tilde S}}}_D} = \mathit{\boldsymbol{W}}_1^{\rm{T}}{\mathit{\boldsymbol{S}}_{\rm{D}}}{\mathit{\boldsymbol{W}}_1} $

(13)

因为差异拉普拉斯散度矩阵S_D也是非负定的，所以对它进行特征值分解后的特征值可能为零，也可能为正。设分解后d个最小的特征值为μ₁≤μ₂≤…≤μ_d，其对应向量分别为φ₁, φ₂, …, φ_d，令W₂= {φ₁, φ₂, …, φ_d}。最终的投影矩阵W =W₁W₂。DAG-DA算法的步骤如下:

输入：训练样本集合{(x_i, y_i)}_i=1ⁿ，令样本矩阵为X=[x₁, x₂, …, x_n]∈ R ^d×n，降维子空间维数r。

输出：投影矩阵W∈ R^d×r。

(1) 构建相似度矩阵S和差异性矩阵B；

(2) 计算局部拉普拉斯散度矩阵S_L，非局部拉普拉斯散度矩阵S_N，以及它们的和S_T与差S_D；S_L=XLX^T, 其中L=S′－S，S′为对角矩阵，$ \mathit{\boldsymbol{S}}{\prime _{ii}} = {\sum _j}{\mathit{\boldsymbol{S}}_{ij}}$；S_N=XNX^T, 其中N=B′－B，B′为对角矩阵，$\mathit{\boldsymbol{B}}{\prime _{ii}} = {\rm{ }}{\sum _j}{\mathit{\boldsymbol{B}}_{ij}};{\mathit{\boldsymbol{S}}_{\rm{T}}} = {\mathit{\boldsymbol{S}}_{\rm{L}}} + {\mathit{\boldsymbol{S}}_{\rm{N}}} $，且S_D=S_L－S_N；

(3) 对矩阵S_T进行特征分解，取其非零特征值λ_i和非零特征值对应的向量u_i，分别组成矩阵Λ和U，其中Λ为对角矩阵且Λ_ii=λ_i，U的第i列为u_i，令W₁=UΛ；

(4) 将S_D在W₁上投影，得到$ {{\mathit{\boldsymbol{\tilde S}}}_{\rm{D}}} = \mathit{\boldsymbol{W}}_1^{\rm{T}}{\mathit{\boldsymbol{S}}_{\rm{D}}}{\mathit{\boldsymbol{W}}_1}$，对${{\mathit{\boldsymbol{\tilde S}}}_{\rm{D}}} $进行特征分解，取最小的r个特征值对应的特征向量，组成一个矩阵W₂；

(5) 得到最终投影矩阵：W=W₁W₂。

2.3 DAG-DA算法和SLDA算法比较

(1) 在用邻接图来保持样本的局部结构信息时，SLDA算法找到近邻后，利用类别信息再构建类内和类间邻接图，这样可能会出现只有类内近邻或只有类间近邻的数据不平衡问题。而在计算权值时，同等距离下类间近邻比类内近邻的权值更大，这显然不利于保持数据的结构。DAG-DA算法分别在类内和类间找近邻点，从热核函数可以看出:若是类内近邻点，则赋予较大的权值；若是类间近邻点，则权值较小。这样可以扩大类内距离缩小类间距离的影响，同时也解决了数据不平衡的问题，更好地保持了数据的流形结构。

(2) 在优化问题选取目标函数时，SLDA算法的目标函数(式(1))分子上仅考虑了将局部拉普拉斯散度矩阵进行投影，而忽略了非局部拉普拉斯散度矩阵; 而DAG-DA算法在最优化函数(式(3))里引入了非局部散度，强调了类间散度，使得降维后的数据可分性更好。

(3) 在对总体拉普拉斯散度矩阵进行特征分解后，SLDA算法选取的投影空间为特征向量乘以一个小于1的数; 而DAG-DA算法直接将非零向量作为投影空间，相当于少乘一个缩小因子，使数据更加准确。

3 实验结果与分析

为验证DAG-DA算法的有效性，在各类数据上对DAG-DA算法进行实验，并与SLDA算法、LPP算法和DNE算法进行比较。这4种算法均可以对近邻个数进行调节，其中SLDA算法和DAG-DA算法中的t和σ为可变参数。本文实验均采用最近邻分类器对降维后的数据进行分类。

3.1 人工数据

随机生成两类符合均匀分布的人工数据，一类是0~1之间的随机数，另一类是0.7~1.7之间的随机数，分别用SLDA算法和DAG-DA算法学习得到投影矩阵。图 1显示了SLDA算法和DAG-DA算法学习投影矩阵后投影后的效果。从图 1中看到，图 1(b)中的类间点距离更大，类内点距离更小，也就是DAG-DA算法相对比SLDA而言，降维效果有了明显的提升，也更适合于分类。在对优化函数进行等价替代后，SLDA算法的目标函数里并没有考虑类间距离，仅将局部拉普拉斯散度矩阵投影到总空间。DAG-DA算法对其进行改进，将局部拉普拉斯散度矩阵与非局部拉普拉斯散度矩阵的差进行投影，是一个更理想的目标函数。

除了图 1能直观地看出两种算法的分类效果外，表 1还给出了两种算法在人工数据测试后的具体结果。从表 1中可以看出，用DAG-DA算法降维比SLDA算法得到的类内距离和更小，类间距离和更大，说明降维后分类效果更佳。

3.2 Yale人脸数据库

Yale人脸数据库^[14]共有165张人脸灰度图像，每张图像数字为32×32的像素矩阵，来自15个不同的人，每人11张。每个人的各张图像有不同的光照条件(中央光照、左侧光照、右侧光照)和不同的表情(眼镜/无眼镜、高兴、常态、悲伤、睡眠、惊讶和眨眼)。Yale人脸数中不同条件下的一些人脸图像样本如图 2所示。该实验中，随机选择60%样本作为训练样本，其余的40%作为测试样本。

在实验时，为了减少大量的计算，先采用PCA方法对测试数据进行初始降维，初始降维的投影矩阵为A。再采用DAG-DA，SLDA，LPP，DNE四种算法进行再次降维，二次降维矩阵为W，最终投影矩阵为P = WA。利用P把训练样本和训练样本投影到低维空间，再用近邻分类器进行分类，得到各个方法降维后的识别率。在构造邻接图时，选取k=3，通过循环找到最佳的t和σ。图 3显示了DAG-DA算法，SLDA算法，LPP算法和DNE算法在不同维数下的识别率。表 2显示了不同算法在Yale人脸库上的最高识别率和特征维度。

3.3 ORL人脸数据库

ORL人脸数据库^[15]共有400张人脸灰度图像，每张图像为32×32的像素矩阵，来自40个不同的人，每人10张。每个人的各图像有不同的人脸大小、人脸旋转程度和不同的表情(睁眼/闭眼、高兴/悲伤、戴眼镜/不戴眼镜等)。ORL人脸数据库中不同条件下的一些人脸图像样本如图 4所示。该实验中，随机选择60%样本作为训练样本，其余的40%作为测试样本。

同样，为了减少大量的计算，先采用PCA方法对测试数据进行初始降维，再采用DAG-DA，SLDA，LPP和DNE四种算法进行再次降维，最后用近邻分类器进行分类，得到各个方法降维后的识别率。本文选取k=3，通过循环找到最佳的t和σ。图 5显示了DAG-DA算法，SLDA算法，LPP算法和DNE算法在不同维数下的识别率。表 3显示了不同算法在ORL人脸库上的最高识别率和特征维度。

从实验结果可以看出，DAG-DA算法达到的最佳识别率最高。SLDA算法的识别率高于LPP算法，因为SLDA算法考虑了非局部信息和样本的类别信息。DNE算法也是监督化学习方法，它构建了类内与类间的邻接图，让同类数据在低维空间紧凑，不同类数据尽可能分开，识别率高于LPP算法。DAG-DA算法不仅考虑了非局部信息和样本类别信息，而且考虑了数据不平衡问题，并在目标函数内加入了非局部拉普拉斯散度矩阵，所以最佳识别率不仅优于LPP算法，也比SLDA和DNE算法更胜一筹。

4 结束语

本文为了处理监督化拉普拉斯判别分析方法中的邻接图数据不平衡问题，以及目标函数中类间信息缺失问题，提出了适用于小样本的双邻接图判别分析算法DAG-DA，并在Yale和ORL人脸数据集上验证了该方法的有效性。

但是，本文算法在实现过程中需要寻找热核函数的最佳参数值，同时去掉了总体拉普拉斯散度矩阵的零空间，这两步使得该算法的时间复杂度比LPP和DNE算法高。下一步工作将尝试寻找时间复杂度较低且降维效果好的算法。