邻域保持嵌入（NPE）

Posted 2022-06-30 woaiml

传统的线性降维方法，如主成分分析（PCA）、因子分析（FA）等，关注的是样本的方差，能学习线性流形的结构，却无法学习非线性流形。而经典的流形学习方法虽然能够学习非线性流形结构，但由于本身属于直推学习，无法进行新样本的泛化。另外一些基于核函数的降维方法，如KPCA，尽管可以处理非线性问题，但又忽略了流形的非线性结构。

NPE 作为局部线性嵌入（LLE）算法的线性逼近，它不仅能够捕捉流形的非线性结构，还保留了线性性质，能够进行新样本的泛化。因此，NPE 在效果令人满意的同时，还能够轻松应对新样本，在工业上得到了广泛的应用。

至于流形的问题，如流形是什么，为什么学习流形，可阅读我之前的博客“学习流形”的初认识。

NPE 算法思想及步骤

NPE 的思想和LLE相同，主要是在降维过程中保持流形的局部线性结构不变，从而来提取数据中的有用信息。这里的局部线性结构是通过重构权重矩阵来表征的，而重构权重矩阵就是邻域内邻居点对结点的线性重构的系数矩阵。

和其他经典流形学习算法类似，NPE 的算法步骤主要分为三步：

1. 构建邻域图

利用 K 最近邻（KNN）算法来构建邻域图。假设共有 m 个样本，则邻域图共有 m 个结点，其中 \\( \\mathbfx _ i \\) 表示第 i 个结点。如果 \\( \\mathbfx _ j \\) 是 \\( \\mathbfx _ i \\) 的 k 个最近邻居之一，那么将这两点连接，反之不连接。

2. 计算权重矩阵（提取流形特征）

设权重矩阵为 \\(\\mathbfW\\)，其中元素 \\(W _ij\\)代表结点 i 和结点 j 之间的边的权重，如果两点之间没有变，则对应的矩阵元素为 0。矩阵 \\(\\mathbfW\\) 的元素值主要通过最小化如下目标函数得到：

\\[ \\min \\sum _i \\left| \\left| \\mathbfx _i - \\sum _j W _ij \\mathbfx _j \\right| \\right| ^2, \\]

其中 \\(\\mathbfW\\) 应满足归一化约束：

\\[ \\sum _j W _ij = 1,j=1,2,...,m. \\]

3. 计算映射

通过求解广义特征向量问题来计算降维的线性映射：

\\[ \\mathbfXMX ^T \\mathbfa=\\lambda \\mathbfXX ^T \\mathbfa, \\]

其中数据集 \\( \\mathbfX = ( \\mathbfx _ 1,...,\\mathbfx _ m ) \\)，矩阵 \\( \\mathbfM = ( \\mathbfI - \\mathbfW ) ^T ( \\mathbfI - \\mathbfW ) \\)，矩阵 \\( \\mathbfI = diag(1,...,1) \\)。

按特征值从小到大的次序（ \\( \\lambda _ 0 \\leq ...\\leq \\lambda _ d-1 \\) ）将求解到的特征向量进行排列 \\( \\mathbfa _ 0,..., \\mathbfa _ d-1 \\)，这样降维后的嵌入坐标则为：

\\[ \\mathbfy _ i = \\mathbfA ^T \\mathbfx _ i, \\]

其中

\\[ \\mathbfA = ( \\mathbfa _0, \\mathbfa _1,...,\\mathbfa _d-1 ). \\]

NPE算法推导

NPE 算法主要是先提取流形的局部线性信息，然后通过保留这种信息来实现降维。具体来说，NPE 首先用局部的线性重构来表示流形局部的线性结构，其形式为均方差。设定样本数为 m ，维数为 n，降维后维数为 d，\\(\\mathbfQ (i)\\) 为样本 i 的 k 个近邻样本集合，则在高维空间中表征重构误差的目标函数为：

\\[ \\phi ( \\mathbfW ) = \\sum ^m _i=1 \\left| \\left| \\mathbfx _i - \\sum _j\\in \\mathbfQ(i) W _ij \\mathbfx _j \\right| \\right| ^2, \\]

由优化目标函数得来的矩阵 \\(\\mathbfW\\) 便包含了流形的局部信息。在降维的过程中，NPE 在参数空间（降维后的低维空间）中以保留和高维原始空间相同的局部线性重构为目标（即在目标函数中使用相同的权重矩阵 \\(\\mathbfW\\) ），来实现线性降维，其低维空间的目标函数为：

\\[ \\Phi ( \\mathbfy ) = \\sum ^m _i=1 \\left| \\left| y _i - \\sum ^m _j=1 W _ij y _j \\right| \\right| ^2. \\]

NPE 的算法推导主要分为两步，一是计算权重矩阵，二是计算映射和低维嵌入。具体推导如下：

首先将第一个目标函数表示为矩阵形式：

\\[ \\phi ( \\mathbfW ) = \\sum ^m _i=1 \\left| \\left| \\mathbfx _i - \\sum _j\\in \\mathbfQ(i) W _ij \\mathbfx _j \\right| \\right| ^2 \\]

\\[ = \\sum ^m _i=1 \\left| \\left| \\sum _j\\in \\mathbfQ(i) W _ij \\mathbfx _i - \\sum _j\\in \\mathbfQ(i) W _ij \\mathbfx _j \\right| \\right| ^2 \\]

\\[ = \\sum ^m _i=1 \\left| \\left| \\sum _j\\in \\mathbfQ(i) W _ij ( \\mathbfx _i - \\mathbfx _j ) \\right| \\right| ^2 \\]

\\[ = \\mathbfW ^T _i ( \\mathbfx _i - \\mathbfx _j )( \\mathbfx _i - \\mathbfx _j ) ^T \\mathbfW _i \\]

\\[ = \\sum ^m _i=1 \\mathbfW ^T _i \\mathbfZ _i \\mathbfW _i \\]

其中 \\(\\mathbfW _ i\\) 是有样本 i 的 k 个近邻点对应的权重系数组成的列向量，矩阵 \\(\\mathbfZ _ i = ( \\mathbfx _ i - \\mathbfx _ j )( \\mathbfx _ i - \\mathbfx _ j ) ^T ,j\\in \\mathbfQ(i)\\)。

考虑到归一化约束，我们再将其约束等式化为矩阵形式：

\\[ \\sum _j\\in \\mathbfQ(i) W _ij = \\mathbfW ^T _i \\mathbf1 _k =1, \\]

其中 \\(\\mathbf1 _k\\) 为 k 维全 1 向量。

利用拉格朗日乘子法将矩阵形式得目标函数和约束等式合并为一个目标：

\\[ \\mathbfL (\\mathbfW) = \\sum ^m _i=1 \\mathbfW ^T _i \\mathbfZ _i \\mathbfW _i + \\lambda'(\\mathbfW ^T _i \\mathbf1 _k - 1). \\]

对 \\(\\mathbfW\\) 求导并令其值为 0，则：

\\[ \\mathbfW _i = -\\dfrac12 \\lambda \\mathbfZ ^-1 _i \\mathbf1 _k. \\]

利用归一化约束，对结果进行归一化，则最终权重矩阵 \\(\\mathbfW _i\\) 为：

\\[ \\mathbfW _i = \\dfrac-\\dfrac12 \\lambda \\mathbfZ ^-1 _i \\mathbf1 _k\\mathbfW ^T _i \\mathbf1 _k \\]

\\[ = \\dfrac\\mathbfZ ^-1 _i \\mathbf1 _k\\mathbf1 ^T _k (\\mathbfZ ^-1 _i) ^T \\mathbf1 _k. \\]

得到权重矩阵后，我们就可以来计算降维的映射了。首先定义：

\\[ z _i = y _i - \\sum ^m _j=1 W _ij y _j, \\]

再将上式转换为矩阵形式：

\\[ \\mathbfz = \\mathbfy - \\mathbfWy = (\\mathbfI - \\mathbfW) \\mathbfy. \\]

则第二个目标函数可简化为：

\\[ \\Phi(\\mathbfy) = \\sum ^m _i=1 \\left( y _i - \\sum ^m _j=1 W _ij y _j \\right) ^2 \\]

\\[ = \\sum ^m _i=1 (z _i) ^2 \\]

\\[ = \\mathbfz ^T \\mathbfz \\]

\\[ = \\mathbfy ^T (\\mathbfI - \\mathbfW) ^T (\\mathbfI - \\mathbfW) \\mathbfy. \\]

由于 NPE 算法是线性降维方法，所以降维后的嵌入可表示为：

\\[ \\mathbfy ^T = \\mathbfa ^T \\mathbfX, \\]

这里 \\(\\mathbfy\\) 只是 m 个样本降到一维的嵌入坐标，\\(\\mathbfa\\)是其对应的映射。

那么第二个目标函数可进一步表示为：

\\[ \\Phi(\\mathbfy) = \\mathbfa ^T \\mathbfX (\\mathbfI - \\mathbfW) ^T (\\mathbfI - \\mathbfW) \\mathbfX ^T \\mathbfa \\]

\\[ = \\mathbfa ^T \\mathbfX \\mathbfM \\mathbfX ^T \\mathbfa. \\]

并且参数空间的嵌入坐标也有约束：

\\[ \\mathbfy ^T \\mathbfy = 1 \\implies \\mathbfa ^T \\mathbfX \\mathbfX ^T \\mathbfa. \\]

同第一个目标函数的求解，这里利用拉格朗日乘子法将第二个目标函数和约束条件合并，并对其求导，可得下列的广义特征向量问题：

\\[ \\mathbfXMX ^T \\mathbfa=\\lambda \\mathbfXX ^T \\mathbfa. \\]

其中 \\(\\mathbfM\\) 可由求得的权重矩阵 \\(\\mathbfW\\) 得出，那么上述的广义特征向量便可以求解。最后解得的特征向量，按次序排列便是映射矩阵了。

代码实现

这里分别给出 python 和 matlab 的代码实现。

python

首先创建一个典型的“瑞士卷”数据集：

import numpy

def swissroll(N=1000):
    tt = numpy.array((5*numpy.pi/4)*(1+2*numpy.random.rand(N)))
    height = numpy.array((numpy.random.rand(N)-0.5))
    noise = 0.0
    X = numpy.array([(tt+noise*numpy.random.randn(N))*numpy.cos(tt), 10*height, (tt+noise*numpy.random.randn(N))*numpy.sin(tt)])
    return X

接着，再利用 shogun 库来进行 NPE 算法的实现：

import modshogun as sg

# load data
feature_matrix = swissroll()
# create features instance
features = sg.RealFeatures(feature_matrix)

# create Neighborhood Preserving Embedding converter instance
converter = sg.NeighborhoodPreservingEmbedding()

# set target dimensionality
converter.set_target_dim(2)
# set number of neighbors
converter.set_k(10)
# set number of threads
converter.parallel.set_num_threads(2)
# set nullspace shift (optional)
converter.set_nullspace_shift(-1e-6)

# compute embedding with Neighborhood Preserving Projections method
embedding = converter.embed(features)

matlab

至于 matlab 代码在何晓飞教授的主页上有，由于篇幅过长，这里就不贴出来了。