随着人类步入“大数据时代”，海量数据带来大量信息的同时，也制造了“维度灾难”。降维成为当前机器学习、模式识别、深度学习等领域的研究热点之一。长期以来，人们假设数据是全局线性结构或高斯分布的，这使得线性降维方法成为降维的主流。例如主成分分析(PCA)^[1]、线性判别分析(LDA)^[2]、多维尺度分析(MDS)^[3]、独立成分分析(ICA)^[4]等。

然而，实际中很多数据并非全局线性，流形学习(Manifold learning)方法应运而生，并被证明在处理非线性结构数据时更加行之有效。经典的流行学习方法包括等距特征映射(ISOMAP)^[5]、局部线性嵌入算法(LLE)^[6]、局部切空间排列算法(LTSA)^[8]等。

LLE作为一种经典的流形学习算法, 由于其易于理解、方便实现、全局最优无需迭代等优点，在人脸图像识别^[9]、肺结节检测^[10]、轴承故障诊断^[11]等方向得到了成功应用。但是，原始的LLE同样存在着诸多的问题：

1) 由于LLE假设数据均匀稠密采样，所以对那些受噪声污染、样本稀疏采样、曲率较大的数据，低维嵌入结果会被严重破坏；

2) 原始的流形学习算法均是无监督的，所以LLE没有很好地利用类别信息，降维结果不利于后期的识别分类；

3) 由于LLE是一种没有明确的映射表达式，所以对于新样本的泛化能力很弱。

针对上述问题，学者们从不同角度提出了改进算法。Donoho等人^[12]提出一种Hessian局部线性嵌入(HLLE)，用Hessian算子代替Laplacian算子，提升了处理不均匀采样数据的能力。Ridder^[13]提出了一种利用类别标签的有监督的局部线性嵌入算法(SLLE)，SLLE在计算邻近点间距离时，依据类别加入一个常数项，使得同类的样本在低维嵌入结果中离得更近，不同类的样本则离得更远。Wu等人^[10]在SLLE基础上提出了SC²SLLE算法，在邻域选择机制中加入了斯皮尔曼相关系数，进一步增强了分类性能。He等人^[14]提出了近邻保持嵌入算法(NPE)，NPE通过对LLE算法进行线性逼近，得到一个转换矩阵，使新样本可以方便地映射到低维空间。文献[15]提出了一种半监督样本外点插值算法(SOSI)，进一步提高了LLE的泛化性能。

现有的抗噪LLE算法^[16]虽在一定程度上提升了鲁棒性，但是处理受强噪声干扰、曲率过大的数据还是不够理想，且仅仅考虑了邻域内的欧式距离因素，忽略了数据本身的流形结构。为了进一步提升LLE算法的抗噪以及处理大曲率数据的能力，以文献[17]提出的权值构造方案为基础并加以改进，提出改进重构权值的局部线性嵌入算法(IRWLLE)，并通过人造数据和人脸数据实验来验证算法的有效性。

首先，LLE算法基本思想是任意一个样本点与其局部邻域内的近邻点构成一个局部线性平面，每一个近邻点对应一个权值，在重构误差最小化的原则下，该样本点可由这些近邻点线性表示。通过重构权值矩阵把高维观测空间与低维嵌入空间联系起来，在低维空间中，固定每个样本和它的近邻点的重构权值不变，通过最小化重构误差来得到高维向量在低维空间的嵌入结果。给定$D$维原始数据集$\mathit{\boldsymbol{X}} = \left( {{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots, {\mathit{\boldsymbol{x}}_N}} \right) \in {{\bf{R}}^{D \times N}}$采样于$d $维流形$M $, $d<D $, $N$为样本个数，设$N$的全局低维坐标为$\mathit{\boldsymbol{Y}} = ({\mathit{\boldsymbol{y}}_1}, {\mathit{\boldsymbol{y}}_2}, \cdots, {\mathit{\boldsymbol{y}}_N}) \in {{\bf{R}}^{d \times N}}$。

对于每一样本${\mathit{\boldsymbol{x}}_i}\left( {i = 1, 2, \cdots, N} \right)$，首先选择与其距离最近的$k$个样本作为局部邻域${\mathit{\boldsymbol{X}}_i} = \left( {{\mathit{\boldsymbol{x}}_{i1}}, {\mathit{\boldsymbol{x}}_{i2}}, \cdots, {\mathit{\boldsymbol{x}}_{ik}}} \right) \in {{\bf{R}}^{D \times k}}$，${\mathit{\boldsymbol{x}}_{ij}}\left( {j = 1, 2, \cdots, k} \right)$为${\mathit{\boldsymbol{x}}_i}$的第$j$个近邻样本。LLE认为${\mathit{\boldsymbol{x}}_i}$可由其近邻点线性重构，即${\mathit{\boldsymbol{x}}_i} = {w_{i1}}{\mathit{\boldsymbol{x}}_{i1}} + {w_{i2}}{\mathit{\boldsymbol{x}}_{i2}} + \cdots + {w_{ik}}{\mathit{\boldsymbol{x}}_{{\rm{ik}}}}$，其中${w_{ij}}\left( {j = 1, 2, \cdots, k} \right)$为近邻点${\mathit{\boldsymbol{x}}_{ij}}$所对应的重构权值，利用最小化重构误差来计算${\mathit{\boldsymbol{x}}_i}$对应的重构权值向量${\mathit{\boldsymbol{w}}_i} = {\left[{{w_{i1}}, {w_{i2}}, \ldots, {w_{ik}}} \right]^{\bf{T}}}$，即

$ \begin{array}{l} \mathop {{\rm{arg\ min}}}\limits_{{w_i}} {\left\| {{\mathit{\boldsymbol{x}}_i} - \sum\limits_{j = 1}^k {{w_{ij}}{\mathit{\boldsymbol{x}}_j}} } \right\|^2}\\ \;\;\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\sum\limits_{j = 1}^k {{w_{ij}} = 1} \end{array} $

(1)

式中，${{w_{ij}}}$是${{\mathit{\boldsymbol{x}}_{ij}}}$对应的重构权值。

若样本点${\mathit{\boldsymbol{x}}_j}$不是${\mathit{\boldsymbol{x}}_i}$的近邻，则${w_{ij}} = 0$。令$\mathit{\boldsymbol{B}} = \left[{\left( {{\mathit{\boldsymbol{x}}_1}-{\mathit{\boldsymbol{x}}_i}} \right), \left( {{\mathit{\boldsymbol{x}}_2}-{\mathit{\boldsymbol{x}}_i}} \right), \cdots, \left( {{\mathit{\boldsymbol{x}}_k}-{\mathit{\boldsymbol{x}}_i}} \right)} \right] $，再另$\mathit{\boldsymbol{G}} = {\mathit{\boldsymbol{B}}^{\rm{T}}}\mathit{\boldsymbol{B}}$为半正定$k$阶Gram方阵，${\mathit{\boldsymbol{e}}_k}$为$k$维全1列向量，由约束条件$\sum\limits_{j = 1}^k {{w_{ij}} = 1} $，则式(1)可化为

$\begin{array}{l} {\rm{arg}}\;\mathop {{\rm{min}}}\limits_{{w_i}} {\left\| {\mathit{\boldsymbol{B}}{\mathit{\boldsymbol{w}}_i}} \right\|^2} = \mathop {{\rm{arg}}\;{\rm{min}}}\limits_{{w_i}} (\mathit{\boldsymbol{w}}_i^{\rm{T}}\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{w}}_i})\\ \;\;\;\;\;\;\;\;{\rm{s}}{\rm{.t}}{\rm{.}}\;\;\;\mathit{\boldsymbol{e}}_k^{\rm{T}}{\mathit{\boldsymbol{w}}_i} = 1 \end{array} $

(2)

利用Lagrange乘数法计算此条件极值，得到邻域${\mathit{\boldsymbol{X}}_i}$对应的权重${\mathit{\boldsymbol{w}}_i} = {\mathit{\boldsymbol{e}}_k}\mathit{\boldsymbol{Ge}}_k^{\rm{T}}{\mathit{\boldsymbol{G}}^{ - 1}}{\mathit{\boldsymbol{e}}_{k}}$，实际操作中，只需先求解$\mathit{\boldsymbol{G}}{\mathit{\boldsymbol{w}}_i} = {\mathit{\boldsymbol{e}}_k}$，然后再将解向量${\mathit{\boldsymbol{w}}_i}$归一化，得到的上述权值${\mathit{\boldsymbol{w}}_i}$。对每一个样本重复以上步骤，从而得到总体的权重矩阵$\mathit{\boldsymbol{W}} = \left[{{\mathit{\boldsymbol{w}}_1}, {\mathit{\boldsymbol{w}}_2}, \cdots, {\mathit{\boldsymbol{w}}_N}} \right]$。

LLE算法以局部线性结构为出发点，以重构权值矩阵作为高维观测空间与低维嵌入空间之间联系的纽带，并且低维嵌入结果具有平移、旋转、缩放不变性。LLE算法最后归结为稀疏矩阵的特征分解问题，具有解析的全局最优解，无需迭代，在保持局部特性的流形学习算法中具有里程碑式的意义。但LLE算法也存在抗噪性差、无法处理大曲率数据、对$k$值选取敏感、对数据要求高(匀称稠密采样)等问题。

针对原始LLE算法存在的不足，为了克服原始的权值只考虑距离因素而忽视结构因素的缺点，采用测地线(Geodesic)距离来描述结构，引出一种新的权值构造方法，加入结构权值$w_{ij}^S$和距离权值$w_{ij}^D$, 定义重构权值${w_{ij}}$为$w_{ij}^D$与$w_{ij}^S$乘积，即

$ {w_{ij}} = w_{ij}^S \cdot w_{ij}^D $

(3)

$ w_{ij}^S = {\mathit{D}_E}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)/{D_G}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right) $

(4)

$ w_{ij}^D = {\rm{exp}}\left[{-{D_G}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)/{d_m}} \right] $

(5)

式中，${D_E}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)$、${D_G}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)$分别为样本${\mathit{\boldsymbol{x}}_i}$与其邻域内样本点${\mathit{\boldsymbol{x}}_{ij}}$间的欧氏距离和测地距离，${d_m}$为邻域内所有测地距离的中值。本文测地距离的计算选用经典的Dijkstra算法，其常用于ISOMAP算法之中。

对于距离权值$w_{ij}^D$来说，在某局部邻域中${d_m}$是固定的，某个近邻点$ {\mathit{\boldsymbol{x}}_{ij}}$距样本点${\mathit{\boldsymbol{x}}_{i}}$的测地距离越远，则$ - {D_G}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)/{d_m}$就越小，这个近邻点${\mathit{\boldsymbol{x}}_{ij}}$所对应的$w_{ij}^D$也就越小，这符合“近邻点离得越远，对${\mathit{\boldsymbol{x}}_{i}}$重构的贡献值越小”的思想，除以中值${d_m}$在一定程度上减少了噪声对权值的影响。分析结构权值$w_{ij}^S$, 选择欧氏距离${D_E}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)$与测地距离${D_G}\left( {{\mathit{\boldsymbol{x}}_i}, {\mathit{\boldsymbol{x}}_{ij}}} \right)$的比值来衡量${\mathit{\boldsymbol{x}}_i}$的局部邻域的线性程度，近邻点${\mathit{\boldsymbol{x}}_{ij}}$越偏离线性平面其对重构${\mathit{\boldsymbol{x}}_i}$的贡献度越小，则该近邻点的结构权值$w_{ij}^S$就越小，进一步强调了结构对权值的重要性，提升了抗噪性。

IRWLLE算法分为以下3个环节:

1) 选择局部邻域。对于每一样本${\mathit{\boldsymbol{x}}_i}\left( {i = 1, 2, \cdots, N} \right)$，根据$k$近邻方法以欧氏距离为测度，选择与其距离最近的$k$个样本作为局部邻域${\mathit{\boldsymbol{X}}_i} = \left( {{\mathit{\boldsymbol{x}}_{i1}}, {\mathit{\boldsymbol{x}}_{i2}}, \cdots, {\mathit{\boldsymbol{x}}_{ik}}} \right) \in {{\bf{R}}^{D \times k}}$，其中${\mathit{\boldsymbol{x}}_{ij}}\left( {j = 1, 2, \cdots, N} \right)$为${\mathit{\boldsymbol{x}}_i}$的第$j$个近邻。

2) 计算重构权值矩阵。对于每一个局部邻域${\mathit{\boldsymbol{X}}_i}$，计算出邻域内的欧氏距离和测地距离，并按式(3)—式(5)计算相应的新重构权值矩阵${\mathit{\boldsymbol{w}}_i}$。重复以上步骤得到总体的重构权值矩阵$\mathit{\boldsymbol{W}}$。

3) 计算低维嵌入。固定之前到的重构权值矩阵$\mathit{\boldsymbol{W}}$，计算输入样本集${\mathit{\boldsymbol{X}}_i}$的低维嵌入坐标$\mathit{\boldsymbol{Y}}$。通过最小化重构误差来计算$\mathit{\boldsymbol{Y}}$，即

$ \begin{array}{l} \mathop {\;\;\;{\rm{arg\ min}}}\limits_Y {\sum\limits_{i = 1}^N {\left\| {{y_i} - \sum\limits_{j = 1}^k {{w_{ij}}{\mathit{\boldsymbol{y}}_j}} } \right\|} ^2}\\ {\rm{s}}{\rm{.t}}{\rm{.}}\;\;\sum\limits_{i = 1}^N {{\mathit{\boldsymbol{y}}_i} = 0}, \frac{1}{N}\sum\limits_{i = 1}^N {{\mathit{\boldsymbol{y}}_i}\mathit{\boldsymbol{y}}_i^{\rm{T}} = {\mathit{\boldsymbol{I}}_d}} \end{array} $

(6)

式中，${{\mathit{\boldsymbol{I}}_d}}$为$d$阶单位矩阵。

令${\mathit{\boldsymbol{y}}_i} = \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{I}}_i}, \mathit{\boldsymbol{M}} = {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)$，则式(6)可化为

$ \mathop {{\rm{arg\ min}}}\limits_Y {\sum\limits_{i = 1}^N {\left\| {\mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{I}}_i} - \mathit{\boldsymbol{Y}}{\mathit{\boldsymbol{w}}_i}} \right\|} ^2} = \mathop {{\rm{arg\ min}}}\limits_Y \;{\rm{tr}}\left( {\mathit{\boldsymbol{YM}}{\mathit{\boldsymbol{Y}}^{\rm{T}}}} \right) $

(7)

则上述问题最终转化为求解矩阵$\mathit{\boldsymbol{M}}$的特征分解问题，即

$ \mathit{\boldsymbol{M}}{\mathit{\boldsymbol{Y}}^{\rm{T}}} = \lambda {\mathit{\boldsymbol{Y}}^{\rm{T}}} $

(8)

求得$\mathit{\boldsymbol{M}}$的最小的$d+1$个特征值对应的特征向量，注意到最小特征值为0，它所对应的特征向量为全1向量，所以必须剔除，剩下的$d$个特征向量即为低维嵌入结果$\mathit{\boldsymbol{Y}}$。

经过以上分析，IRWLLE算法的流程为

输入：参数$k、d$，高维数据集

$ \mathit{\boldsymbol{X}} = \left( {{\mathit{\boldsymbol{x}}_1}, {\mathit{\boldsymbol{x}}_2}, \cdots, {\mathit{\boldsymbol{x}}_N}} \right) \in {{\bf{R}}^{D \times N}} $

输出：低维嵌入结果

$ \mathit{\boldsymbol{Y}} = \left( {{\mathit{\boldsymbol{y}}_1}, {\mathit{\boldsymbol{y}}_2}, \cdots, {\mathit{\boldsymbol{y}}_N}} \right) \in {{\bf{R}}^{d \times N}} $

1)选择局部邻域。对于某个样本${{\mathit{\boldsymbol{x}}_i}}$，根据欧氏距离选取$k$个最近临点构成局部邻域${{\mathit{\boldsymbol{X}}_i}}$。

2) 计算重构权值矩阵。固定${{\mathit{\boldsymbol{X}}_i}}$，分别计算结构权值$w_{ij}^S$和距离权值$w_{ij}^D$得到${{\mathit{\boldsymbol{w}}_i}}$，重复以上过程得到总体的重构权值矩阵${\mathit{\boldsymbol{W}}}$。

3) 计算低维嵌入。根据${\mathit{\boldsymbol{W}}}$，求矩阵$\mathit{\boldsymbol{M}} = {\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)^{\rm{T}}}\left( {\mathit{\boldsymbol{I}} - \mathit{\boldsymbol{W}}} \right)$的特征值分解，取第2到$d+1$小的特征值对应的特征向量组成低维嵌入${\mathit{\boldsymbol{Y}}}$。

为了验证IRWLLE算法的有效性，分别对人造数据和人脸数据进行实验，并与原始LLE算法和Hessian局部线性嵌入(HLLE)算法进行比较。

选取了流行学习经典的Swiss roll、S-curve和Helix这3种人工数据集来作为实验对象，通过与LLE算法和HLLE算法对比，分析改进算法的性能。为了验证各个算法的抗噪性和处理稀疏采样数据的能力，将3种人工数据分为400和1 000数据点、原始和加噪4种数据。

首先对于Swiss roll数据，它是一个3维的卷筒结构，好的流形算法可以恢复出Swiss roll的几何结构，例如卷曲曲面或平面，同时还能保持各个数据点的近邻关系，即颜色的渐变性。

图 1是Swiss roll数据通过各算法的处理结果。通过对比可以看出，IRWLLE算法面对稀疏采样的数据虽在一定程度上存在数据点重叠现象，但总体上仍然很好地保持了数据点的近邻关系，颜色渐变平稳，没有出现颜色交杂现象，而且对Swiss roll的结构展开也比较理想。而原始的LLE算法和HLLE算法无论怎样调节近邻参数$k$的取值，虽然也将数据降维2维，但是数据点间的近邻关系被严重破坏，颜色混杂，不利于分类学习。因此，可以看出IRWLLE算法在处理稀疏采样的数据时，对比原LLE算法具有明显的优势。

图 2是400个混杂有高斯白噪声的数据通过各算法的处理结果。对比各结果发现, LLE和HLLE算法得到的结果仍然很不理想，噪声的加入使得数据点间的近邻关系更加混乱；反观IRWLLE算法，对比图 1(b)和图 2(b)发现噪声对IRWLLE算法的结果几乎未造成影响，仍然具有很强的鲁棒性，这主要是由于噪声数据过于偏离流形，结构权值减少了其权重。

对于S-curve数据和Helix曲线，限于篇幅只展示IRWLLE算法的抗噪能力，故固定取1 000个数据点，进行实验分析。

图 3(a)是原始的S-curve数据，图 3(b)是加入强噪声之后的S-curve数据，图 3(c)-(e)是各个算法对图 3(b)的处理结果。对比发现，加入强噪声之后LLE算法结果的近邻关系被破坏较为严重，边缘的深蓝色和中部的红色交杂在一起，降维效果最差；HLLE由于引入了Hessian算子在一定程度上减小了噪声数据的干扰，但是被错误排列的数据点仍然较多；观察图 3(c)可以发现，IRWLLE算法面对噪声数据最为鲁棒，3维“S”形流形被充分展开，而且近邻关系保持完好，只有极个别深蓝色点被错误排列，效果最佳。而且在实验中LLE和HLLE算法均对$k$值的选取敏感，不同取值的效果差距明显，参数选择困难，但是由于IRWLLE方法采用比值的形式，距离很远的近邻点权值很小，所以降维结果对$k$值的选取依赖性较小。

图 4(a)是原始的Helix数据，图 4(b)是加入强噪声之后的Helix数据，肉眼以难以区分颜色的渐变关系。图 4(c)-(e)是各算法对图 4(b)的处理结果。对比图 4(b)(c)可以看出，虽然加噪数据已经十分“嘈杂”，但是IRWLLE算法依然能够发现数据的内在流形结构，对噪声干扰的鲁棒性非常强并且最大程度上保持了数据点间的近邻关系；反观图 4(d)(e)可以清楚地发现，LLE、HLLE算法对这种曲率极大且加入噪声的数据没有处理能力，不仅流形结构被严重破坏，而且近邻关系也显得杂乱无章。

采用ORL和Yale人脸数据库进行人脸识别实验。ORL人脸数据库由剑桥大学AT&T实验室采集和整理，共400幅灰度人脸图像，每幅图像分辨率为92×112像素，包含40个人，每人10幅图像，呈现不同的表情、角度以及面部饰品的变化。选取每人的前5幅图像作为训练集，后5幅图像作为测试集。Yale人脸数据库包含了15人，每人11幅，共165幅人脸图像，图像分辨率为100×100像素。每人10幅图像，前5幅图像作为训练集，后5幅图像作为测试集。图 5、图 6分别列举一些了ORL和Yale数据样本，其中每一行为一名志愿者，每行的前5幅为训练图像，后5幅为测试图像。

分别采用LLE和IRWLLE算法对整理好的ORL和Yale图像数据进行降维，并采用最近邻分类器进行人脸识别，改变参数$k$和$d$来寻找每种降维算法所能达到的最大识别率，其中$k$是某样本${\mathit{\boldsymbol{x}}_i}$的近邻值个数，$d$是高维数据经过降维之后的维数，识别率为识别正确的测试样本数与总测试样本数之比。

表 1展示了LLE和IRWLLE算法在ORL和Yale人脸数据上的降维和识别效果的比较。对比发现，IRWLLE算法在识别率上较原始LLE有了显著提升。

LLE算法是一种广泛研究和应用的流形降维算法，针对LLE算法存在的抗噪性差、难以处理大曲率数据、难以处理稀疏采样数据和对参数选择敏感问题，提出一种改进的局部线性嵌入算法IRWLLE，引入一种有别于原始LLE算法的重构权值矩阵的计算方法，将测地距离与欧氏距离之比定义为结构权值；将测地距离与中值测地距离之比定义为距离权值，通过结构权值与距离权值的乘积综合考虑了流形的结构和距离因素，思想明确易于实现。

通过对国际通用的经典人工数据和人脸数据的实验，表明IRWLLE算法具有良好的降维效果和识别率，有效地克服了原有算法的不足。但需指出的是，该改进方法并未引进样本的类别信息，也未考虑样本外问题，这些都是以后工作的重点。