Laplacian Eigenmaps.

拉普拉斯特征映射(Laplacian Eigenmaps,LE)是一种基于图的降维算法，它从局部的角度去构建数据之间的关系，希望相互间有关系的点（即在图中相连的点）在降维后的空间中尽可能的靠近，从而在降维后仍能保持原有的数据结构。

1. 拉普拉斯特征映射的建模

拉普拉斯特征映射通过构建邻接矩阵为 $W$ 的图来重构数据流形的局部结构特征。$W$的构建形式是灵活的，比如使用kNN算法将每个点最近的K个点连上边。

根据假设，如果两个数据样本 $i$ 和 $j$ 很相似，那么 $i$ 和 $j$ 在降维后目标子空间中也应该尽量接近。设数据样本的总数为 $n$ ，目标子空间(即最终的降维目标)的维度为 $m$。

定义 $n×m$ 大小的降维目标矩阵$Y$，其中每一个行向量$y_i^\top$是数据样本 $i$ 在目标 $m$ 维子空间中的向量表示（即降维后的数据样本 $i$）。目标是让相似的数据样本 $i$ 和 $j$ 在降维后的目标子空间里仍尽量接近，故拉普拉斯特征映射优化的目标函数如下:

\[\mathop{\min}\quad \sum_{i,j} \mid\mid y_i-y_j \mid\mid^2 \cdot w_{ij}\]

2. 拉普拉斯特征映射的求解

求解上述目标函数：

\[\begin{aligned} & \sum_{i,j} \mid\mid y_i-y_j \mid\mid^2 \cdot w_{ij} \\ = & \sum_i \sum_j \left(y_i^\top y_i-2y_i^\top y_j+y_j^\top y_j\right)\cdot w_{ij} \\ = & \sum_i \left(y_i^\top y_i\right)\cdot \sum_j w_{ij}+\sum_j \left(y_j^\top y_j\right)\cdot \sum_i w_{ij} - 2\sum_i \sum_j y_i^\top y_j \cdot w_{ij} \\ = & 2\sum_i \left(y_i^\top y_i\right)\cdot d_{ii}- 2\sum_i \sum_j y_i^\top y_j \cdot w_{ij} \\ = & 2 \cdot \text{Trace}\left(Y^\top D Y \right)-2 \cdot \text{Trace}\left(Y^\top W Y \right) \\ = & 2 \cdot \text{Trace}\left[Y^\top (D-W) Y \right] \\ = & 2 \cdot \text{Trace}\left(Y^\top L Y \right) \end{aligned}\]

其中 $W$ 为图的邻接矩阵，对角矩阵$D$是图的度矩阵（$d_{ii}=\sum_j w_{ij}$），$L=D-W$为图的拉普拉斯矩阵。

为防止上述优化陷入平凡解$\forall y_i=\bold{0}$，对$Y$进行一些约束：如果降维的维度是$m$，则要求$\text{span}(y_1,y_2,…,y_n) \in \Bbb{R}^m$，即希望降维得到的所有样本的嵌入向量$y$能够尽可能地去填充$\Bbb{R}^m$空间。

由于度矩阵$D$可以衡量每个样本在图中的“重要性”，因此使用$D$对$Y$进行限制，即引入$Y^\top D Y = I$的约束，从而删除嵌入$y$中的任意比例因子，使每个$y_i$的尺度被固定住。

此时拉普拉斯特征映射优化的目标函数如下:

\[\begin{aligned} \mathop{\min}\quad &\text{Trace}\left(Y^\top L Y \right) \\ \text{s.t.}\quad &Y^\top D Y = I \end{aligned}\]

用拉格朗日乘子法对目标函数求解:

\[\begin{aligned} L(Y,\lambda) = & \text{Trace}\left(Y^\top L Y \right) - \sum_i \lambda_i \left(y_i^\top D y_i-I\right) \\ = & \text{Trace}\left(Y^\top L Y \right) - \text{Trace}\left(\Lambda \left(Y^\top D Y-I\right)\right) \\ \frac{\partial L}{\partial Y} = & (L + L^T)Y - \left(D Y\Lambda+ D^\top Y\Lambda^\top\right) \\ = & 2LY - 2D Y\Lambda= 0 \\ \Rightarrow LY = & D Y\Lambda \end{aligned}\]

对于任意向量$y_i$，上式可写为$Ly_i=\lambda_iDy_i$，这是一个广义特征值问题（$D^{-1}Ly_i=\lambda_iy_i$）。

将最优解$LY = D Y\Lambda$带回目标函数得（注意到$Y^\top D Y=I$）：

\[\begin{aligned} \text{Trace}\left(Y^\top L Y \right) = \text{Trace}\left(Y^\top D Y\Lambda \right)= \text{Trace}\left(\Lambda \right) \end{aligned}\]

最小化上述目标函数，即最小化特征值$\sum_i \lambda_i$之和。因此选择对$D^{-1}L$进行特征值分解后的$m$个最小非零特征值所对应的特征向量作为$Y$，即可达到降维的目的。

3. 拉普拉斯特征映射算法

拉普拉斯特征映射的完整流程如下：

\[W_{ij} = \exp\left(-\frac{\|x_i-x_j\|^2}{t}\right)\] 
def cal_pairwise_dist(x):
    '''计算任意两个点之间距离的平方
    (a-b)^2 = a^2 + b^2 - 2*a*b
    '''
    sum_x = np.sum(np.square(x), 1)
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    return dist

def cal_rbf_dist(data, n_neighbors = 10, t = 1):
    dist = cal_pairwise_dist(data)
    dist[dist < 0] = 0
    n = dist.shape[0]
    rbf_dist = np.exp(-(dist/t))

    W = np.zeros((n, n))
    for i in range(n):
        index_ = np.argsort(dist[i])[1:1+n_neighbors]
        W[i, index_] = rbf_dist[i, index_]
        W[index_, i] = rbf_dist[index_, i]

    return W

def le(data,
          n_dims = 2,
          n_neighbors = 5, t = 1.0):
    '''
    :param data: (n_samples, n_features)
    :param n_dims: target dim
    :param n_neighbors: k nearest neighbors
    :param t: a param for rbf
    :return:
    '''
    N = data.shape[0]
    W = cal_rbf_dist(data, n_neighbors, t)
    D = np.zeros_like(W)
    for i in range(N):
        D[i,i] = np.sum(W[i])

    D_inv = np.linalg.inv(D)
    L = D - W
    eig_val, eig_vec = np.linalg.eig(np.dot(D_inv, L))

    sort_index_ = np.argsort(eig_val)
    eig_val = eig_val[sort_index_]

    j = 0
    while eig_val[j] < 1e-6:
        j+=1

    sort_index_ = sort_index_[j:j+n_dims]
    eig_val_picked = eig_val[j:j+n_dims]
    eig_vec_picked = eig_vec[:, sort_index_]

    return eig_vec_picked