深入理解拉普拉斯特征映射

拉普拉斯特征映射(Laplacian Eigenmaps，LE)是一种降维方法，之前有讲过一种比较常见的降维算法： PCA。

LE在图嵌入中有一些应用，所以这里总结一下。

1. 一些定义

$G=(V,E)$，$|V|=n$，$W$为邻接矩阵，$D$为图的度矩阵（对角阵），有$D_{ii}={\sum }_{j=1}^n{w}_{ij}$。

图的拉普拉斯矩阵定义：
$$L=D-W$$
需要注意的是，对于有向图来说，$W$可以为出度矩阵，也可以为入度矩阵，视具体情况而定。我们下面讨论的都是无向图的拉普拉斯矩阵。

对于实对称矩阵$L$，有以下性质：对任意非零列向量$y=[y_1,...,y_n{]}^T$，我们有：
$$y^{T}Ly=\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_i-y_j{)}^2\ge 0$$
即$L$是半正定的。证明如下：
$$\begin{array}{rl}{y}^{T}Ly& =y^{T}Dy-y^{T}Wy\\ & =\sum _{i=1}^n{D}_{ii}{y}_i^2-\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}\\ & =\frac{1}{2}(\sum _{i=1}^n{D}_{ii}{y}_i^2-2\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}+\sum _{j=1}^n{D}_{jj}{y}_j^2)\\ & =\frac{1}{2}(\sum _{i=1}^n(\sum _{j=1}^n{w}_{ij})y_i^2-2\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}+\sum _{j=1}^n(\sum _{i=1}^n{w}_{ij})y_j^2)\\ & =\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n{w}_{ij}(y_i-y_j{)}^2\end{array}$$

2. 矩阵的迹求导

对LE的推导需要补充一些数学知识：
$$\begin{gathered} \frac{\partial tr(A^{T}BA)}{\partial A}=(B+B^T)A \\ \frac{\partial tr(AB{A}^{T}C)}{\partial A}=CAB+C^{T}A{B}^T \end{gathered}$$

3. LE

3.1 目标函数

在图嵌入中，我们的目的是得到每个顶点的低维向量表示，假设嵌入向量维度为$m$。

LE的思想：如果两个数据实例$i$和$j$很相似，那么它们在降维后的目标子空间中应该尽量接近。

而在图嵌入中，衡量两个节点是否相似的最直接度量为一阶邻近度，两个节点$v_i$和$v_j$间的一阶邻近度即两个节点间相连边的权重。

一阶邻近度通常意味着真实世界网络中两个节点的相似性。正所谓近朱者赤近墨者黑，比如在社交网络中成为朋友的人往往有相似的兴趣爱好，又比如万维网上相互链接的网页往往谈论类似的话题。由于这种重要性，现有的许多图嵌入算法在设计目标函数时都会保持一阶邻近度。

因此，LE的思想可以重新表述为：如果两个节点在原始图中彼此相邻，那么它们的低维嵌入向量应该彼此接近。

我们设节点$v_i$和$v_j$的向量表示分别为$y_i\in R^m$和$y_j\in R^m$，鉴于此，我们就可以得出LE的目标函数：
$$min\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n||y_i-y_j|{|}^2{w}_{ij}$$
从目标函数可以看出，如果两个节点$v_i$和$v_j$的一阶邻近度$w_{ij}$较大（很相似），那么我们在最小化目标函数时，就会更多地考虑二者间的差异。

假设$y_i$是一个列向量，那么$y_i^2=y_i^T{y}_i$。

下面开始对目标函数进行化简：
$$\begin{array}{rl}& \frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n||y_i-y_j|{|}^2{w}_{ij}\\ =& \frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n(y_i^T{y}_i{w}_{ij}-2y_i{y}_j{w}_{ij}+y_j^T{y}_j{w}_{ij})w_{ij}\\ =& \frac{1}{2}\sum _{i=1}^n(\sum _{j=1}^n{w}_{ij})y_i^T{y}_i-\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}+\frac{1}{2}\sum _{j=1}^n(\sum _{i=1}^n{w}_{ij})y_j^T{y}_j\\ =& \frac{1}{2}\sum _{i=1}^n{D}_{ii}{y}_i^T{y}_i-\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}+\frac{1}{2}\sum _{i=1}^n{D}_{ii}{y}_i^T{y}_i\\ =& \sum _{i=1}^n{D}_{ii}{y}_i^T{y}_i-\sum _{i=1}^n\sum _{j=1}^n{y}_i{y}_j{w}_{ij}\\ =& tr(Y^{T}DY)-tr(Y^{T}WY)\\ =& tr(Y^T(D-W)Y)\\ =& tr(Y^{T}LY)\end{array}$$
其中，$Y=(y_1,...,y_n{)}^T$，$tr(A)$表示对矩阵$A$求秩（对角线之和）。

3.2 约束条件

考虑到目标函数：
$$min\frac{1}{2}\sum _{i=1}^n\sum _{j=1}^n||y_i-y_j|{|}^2{w}_{ij}$$
我们不妨设想一种极端情况，假设所有节点都映射到了同一个位置，也就是所有节点的嵌入向量$y$相同，那么此时目标函数肯定有最小值0。又比如我们就假设所有节点的嵌入向量全部为0向量，此时目标函数也有最小值0。

以上两种情况是毫无意义的。此外，在上述情况下，$y$的维度也是任意的。因此，为了得到唯一解，我们需要对节点最终的嵌入向量$y$做一些限制。

一个最简单的限制就是：我们希望最终得到的所有节点的嵌入向量$y$能够尽可能地去填充$R^m$空间，而不是挤在一起。举个例子，假设$m=2$，那么我们不希望所有节点的嵌入向量$y$(一个点)挤在一起，或者就只是聚集在某一条直线附近，而是尽可能的去散开，因为如果聚在一起就表示图中所有节点具有类似的性质，这显然是不科学的，我们只是需要邻接的节点的嵌入向量相似即可。

为了理解限制条件，我查阅了拉普拉斯特征映射的论文，论文中给出的限制条件为：
$$Y^{T}DY=I$$
原文中关于约束条件$Y^{T}DY=1$的描述为：$Y^{T}DY=1$这一限制条件删除了嵌入中的任意比例因子。也就是说，每一个$y_i$的尺度被固定了。

至于为什么选则$D$而不是其他矩阵，原文中给出的解释为：度矩阵$D$可以很好地提供节点自身的一个度量。我们可以这样理解：$D$可以描述每一个节点在graph中的重要性，如果我们使用$D$对$y$的尺度进行限制，那么我们可以使得$y$在$R^m$上分布更加均匀，也就是前面我们提到的：我们希望最终得到的所有节点的嵌入向量$y$能够尽可能地去填充$R^m$空间，而不是挤在一起。

3.3 优化

有了$y_i^{T}D{y}_i=1$这一限制条件后，我们可以构造拉格朗日函数如下：
$$L=tr(Y^{T}LY)-\sum _{i=1}^n{\lambda }_i(1-y_i^{T}D{y}_i)$$
我们令$\Lambda$为一个对角阵，对角阵上的元素${\Lambda }_{ii}={\lambda }_i$，因此，上述拉格朗日函数可以变为矩阵形式：
$$L=tr(Y^{T}LY)-tr(\Lambda (Y^{T}DY-I))$$
对于第一项，由第二节中求导公式(1)可知：
$$\frac{\partial tr(Y^{T}LY)}{\partial Y}=(L+L^T)Y$$
对于第二项，由第二节中求导公式(2)可知：
$$\frac{\partial tr(\Lambda Y^{T}DY)}{\partial Y}=DY\Lambda +D^{T}Y{\Lambda }^T$$
因此我们有：
$$\begin{array}{rl}\frac{\partial L}{\partial Y}& =LY+L^{T}Y+DY\Lambda +D^{T}Y{\Lambda }^T\\ & =2LY+2DY\Lambda =0\end{array}$$
解得：
$$LY=-DY\Lambda$$
写成一般形式为：$L{y}_i={\lambda }_{i}D{y}_i$

3.4 广义特征值问题

上述问题转为广义特征值求解问题：
$$Ly=\lambda Dy$$

所谓广义特征值问题：对于实对称矩阵$A$和实对称正定矩阵$B$，如果有非零解$x$满足：
$$Ax=\lambda Bx$$
则称$\lambda$是$A$相对于$B$的特征值。

在本文中，$L$是实对称矩阵,$D$首先是实对称矩阵，其次$D$也是正定的，这点很好证明，便不再叙述。

那么到底如何求解广义特征值文题呢？我们知道，一般特征值的表达形式为：
$$Ax=\lambda x$$
对于一般特征值求解，我们可以转为：
$$(A-\lambda E)x=0$$

那么我们可以尝试将广义特征值问题转为此类问题。

首先，我们将矩阵$B$进行Cholesky分解：
$$B=G{G}^T$$
其中$G$为下三角矩阵。

然后有：
$$Ax=\lambda G{G}^{T}x$$
两边乘上$G^{-1}$：
$$G^{-1}Ax=\lambda G^{T}x$$
然后：
$$\begin{array}{rl}{G}^{-1}Ax& =\lambda G^{T}x\\ G^{-1}A((G^T{)}^{-1}{G}^T)x& =\lambda G^{T}x\end{array}$$
这时两边都出现了$G^{T}x$，我们将其合并：
$$(G^{-1}A(G^T{)}^{-1})(G^{T}x)=\lambda (G^{T}x)$$
我们令：
$$\begin{array}{rl}{G}^{-1}A(G^T{)}^{-1}& =T\\ G^{T}x& =s\end{array}$$
于是原问题变为一般特征值的求解问题：
$$Ts=\lambda s$$
然后，我们就能求出矩阵$T$的特征值$s$，进而解出$x$。

3.5 结果

经过3.4之后，得到了$L{y}_i=\lambda D{y}_i$中的解向量$y_i$，然后选取最小的$m$个非零特征值对应的特征向量作为节点的嵌入向量。

为什么要选取非零特征值的特征向量？根据$L$的性质，$L$的行和为0。因此，从$L{y}_i=\lambda D{y}_i$可以看出，$L$具有广义特征值0和对应的特征向量$1=(1,...,1{)}^T$，由如果$1=(1,...,1{)}^T$被选中，那么每一个节点的嵌入向量中的某一维度将全是1，那么嵌入向量将退化到更低一维的空间中。

为什么要选取$m$个最小的特征值对应的特征向量？参考网上说法：我们将$LY=-DY\Lambda$代回到原目标函数$tr(Y^{T}LY)$中，可以得到：
$$mintr(-Y^{T}DY\Lambda )$$
又由于有限制条件：
$$Y^{T}DY=I$$
所以最终目标函数为：
$$mintr(-\Lambda )$$
这里$tr(\Lambda )={\lambda }_1+...+{\lambda }_n$，目标函数即为特征值之和，为了使得目标函数最小化，当然应该选择最小的$m$个特征值。