LPP（局部保持投影）以及LE（拉普拉斯映射）的区别以及代码python解读

关于LPP与LE在降维上的区别在这篇文章上已经描述的十分清楚了：LE与LPP的简介（强烈建议看完这篇文章在看本文的代码解析）

假设已有数据集样本集合X={ ${\bold x_1}$ , ${\bold x_2}$ ,..., ${\bold x_n}$ }，且每个样本的维度为m。

（一）求解LPP映射矩阵的算法步骤如下：

1.对整个数据集构建邻接矩阵W，即查找每个样本 $x_i$ 在数据集中距离最近的k个样本，如样本 ${\bold x_j}$ 为样本 ${\bold x_i}$ 的k个最近邻接点之一，则在邻接矩阵W中的第i行j列元素 ${\bold W}_{ij}$ 以及第i列j行元素 ${\bold W}_{ji}$ 可通过以下公式计算得出：

${\bold W}_{ij}={\bold W}_{ji}=e^{-\frac{||{\bold x_i}-{\bold x_j}||^2}{t}}$ （1）

其中t为热核参数（设置多少自己看模型表现而定，一般可尝试1-10之间）若样本 ${\bold x_i}$ 不为样本 ${\bold x_j}$ 的k邻接点之一，或样本 ${\bold x_j}$ 不为样本 ${\bold x_i}$ 的k邻接点之一，则：

${\bold W}_{ij}={\bold W}_{ji}=0$ （2）

2.求解下面问题的特征值 $\lambda$ 以及特征向量 ${\bold v}$ ：

${\bold X}{\bold L}{\bold X}^T{\bold v} = \lambda {\bold X}{\bold D}{\bold X}^T {\bold v}$ （3）

其中的D是对角矩阵（diagonal matrix），其中的对角元素 ${\bold D_{ii}}$ 的值为第i行或者第i列的所有元素的值之和。其中的L为拉普拉斯矩阵，也是对称的半正定矩阵， ${\bold L}={\bold D}-{\bold W}$ 。

3.通过求解以上公式（3）中的特征问题，可以得到一系列的特征值以及该特征值对应的特征向量。然后我们根据这一系列的特征值进行升序排列，同时将对应的特征向量按列放入映射矩阵 ${\bold V}$ 中。若我们需要将原始数据从m维降至d维，则我们需要取映射矩阵的前d列。映射公式如下：(映射矩阵 ${\bold V}$ 为m*d的矩阵，样本 ${\bold x_i}$ 为m×1的向量， ${\bold s_i}$ 为从 ${\bold x_i}$ 降维得到的d*1的向量)

${\bold s_i}={\bold V}^T{\bold x_i}$ （4）

4.如果你没看文字开头链接的文章，你可能会有疑问，为什么求公式（3）的特征问题来得到映射矩阵。

公式（3）是由LPP的最小化问题转化而来，该最小化问题如下：

$\underset{{\bold v}}{argmin} {\bold v}^T{\bold X}{\bold L}{\bold X}^T{\bold v}$

$s.t. {\bold v}^T {\bold X}{\bold D}{\bold X}^T{\bold v}=1$

X和L已经确定的情况下，我们通过求解较小的特征值以及对应的特征向量来使该问题最小化。

以下是来自哈工大大佬的python代码块，github地址：https://github.com/heucoder/dimensionality_reduction_alo_codes/tree/master/codes

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from sklearn.datasets import load_digits, load_irisdef rbf(dist, t = 1.0):'''rbf kernel function'''return np.exp(-(dist/t))def cal_pairwise_dist(x):'''计算pairwise 距离, x是matrix(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 distdef cal_rbf_dist(data, n_neighbors = 10, t = 1):dist = cal_pairwise_dist(data)n = dist.shape[0]rbf_dist = rbf(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 Wdef lpp(data,n_dims = 2,n_neighbors = 30, 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)  #建立邻接矩阵W，参数有最近k个邻接点，以及热核参数tD = np.zeros_like(W)for i in range(N):D[i,i] = np.sum(W[i]) #求和每一行的元素的值，作为对角矩阵D的对角元素L = D - W  #L矩阵XDXT = np.dot(np.dot(data.T, D), data)XLXT = np.dot(np.dot(data.T, L), data)'''np.linalg.eig用作计算方阵的特征值以及右特征向量np.linalg.pinv对矩阵进行求逆本人理解：在求解最小特征问题（公式3）时，左乘XDX^T的逆，而后求其特征值以及对应特征向量输出的eig_val,eig_vec为特征值集合以及特征向量集合（此时是无序的）'''eig_val, eig_vec = np.linalg.eig(np.dot(np.linalg.pinv(XDXT), XLXT))'''argsort返回的是特征值在数据eig_val排序后的序号。如数组a=[2,1,7,4],则np.argsort(a)返回的是[1,0,3,2]'''sort_index_ = np.argsort(np.abs(eig_val))eig_val = eig_val[sort_index_]#对特征值进行排序print("eig_val[:10]", eig_val[:10])#此时eig_val已经是升序排列了，需要排除前几个特征值接近0的数，至于为啥请看论文 "Locality Preserving Projections"j = 0while eig_val[j] < 1e-6: j+=1print("j: ", j)#返回需要提取的前n个特征值所对应的n个特征向量的序列号sort_index_ = sort_index_[j:j+n_dims]# print(sort_index_)eig_val_picked = eig_val[j:j+n_dims]print(eig_val_picked)eig_vec_picked = eig_vec[:, sort_index_] #获取该n个特征向量组成的映射矩阵data_ndim = np.dot(data, eig_vec_picked) #公式（4）return data_ndimif __name__ == '__main__':X = load_digits().datay = load_digits().targetdist = cal_pairwise_dist(X)max_dist = np.max(dist)print("max_dist", max_dist)data_2d = lpp(X, n_neighbors = 5, t = 0.01*max_dist)plt.figure(figsize=(12,6))plt.subplot(121)plt.title("LPP")plt.scatter(data_2d[:, 0], data_2d[:, 1], c = y)

（二）接下来是求解LE降维的步骤：

1.第一步构造邻接矩阵W，构造的思路与LPP一致，每个样本均与最近的k个邻接点建立权重关系：

${\bold W}_{ij}={\bold W}_{ji}=e^{-\frac{||{\bold x_i}-{\bold x_j}||^2}{t}}$ （5）

2.求解下面问题的特征值 $\lambda$ 以及特征向量 ${\bold y}$ ：

${\bold L}{\bold y}=\lambda {\bold D}{\bold y}$ (6)

其中，降维后的矩阵的列向量将会来自公式（6）求出的特征矩阵 ${\bold Y}=\{{{\bold y}_1, {\bold y}_2,..., {\bold y}_n}\}$ ，

3.通过求解公式（6），我们将得到的特征值按升序排列，并把特征向量按照对应特征值的大小依次作为列向量放入特征矩阵Y中。若我们需要将数据从m维降至d维，则只需要取特征矩阵的前d列。形成的n*d的矩阵即为降维后的矩阵，其中每一行数据 ${\bold s}_i$ 为样本 ${\bold x}_i$ 降维表示。

4.公式（6）由LE的最小化问题转化而来，该最小化问题表示如下：

$\underset{{\bold y}}{argmin} {\bold y}^T{\bold L}{\bold y}$ （7）

$s.t. {\bold y}^T {\bold D}{\bold y}=1$

在L和D是确定的情况下，我们通过求解公式（6）较小的特征值以及对应的特征向量来是问题（7）最小化。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_digits
from mpl_toolkits.mplot3d import Axes3Ddef rbf(dist, t = 1.0):'''rbf kernel function'''return np.exp(-(dist/t))def cal_pairwise_dist(x):'''计算pairwise 距离, x是matrix(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 distdef cal_rbf_dist(data, n_neighbors = 10, t = 1):dist = cal_pairwise_dist(data)n = dist.shape[0]rbf_dist = rbf(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 Wdef le(data,n_dims = 2,n_neighbors = 5, t = 1.0):''':param data: (n_samples, n_features):param n_dims: 目标维度:param n_neighbors: k个最近邻接点:param t: 热核参数t:return:'''N = data.shape[0]W = cal_rbf_dist(data, n_neighbors, t) #求邻接矩阵WD = np.zeros_like(W)for i in range(N):D[i,i] = np.sum(W[i])#建立对角矩阵D，对角元素均为邻接矩阵第i行或者第i列元素值之和'''本人理解：在求解最小特征问题（公式6）时，左乘矩阵D的逆，而后求其特征值以及对应特征向量输出的eig_val,eig_vec为特征值集合以及特征向量集合（此时是无序的）'''D_inv = np.linalg.inv(D) L = D - Weig_val, eig_vec = np.linalg.eig(np.dot(D_inv, L))sort_index_ = np.argsort(eig_val)eig_val = eig_val[sort_index_]print("eig_val[:10]: ", eig_val[:10])j = 0while eig_val[j] < 1e-6: #排除特征值接近于0，以及该特征值对应的特征向量j+=1print("j: ", j)sort_index_ = sort_index_[j:j+n_dims] #提取特征值大于1e-6开始的序号eig_val_picked = eig_val[j:j+n_dims]  print(eig_val_picked)eig_vec_picked = eig_vec[:, sort_index_]#根据序号从特征向量集合中得到映射后的降维矩阵# print("L: ")# print(np.dot(np.dot(eig_vec_picked.T, L), eig_vec_picked))# print("D: ")# D not equal I ???print(np.dot(np.dot(eig_vec_picked.T, D), eig_vec_picked))#LE不同于LPP，LE不存在Si=YXi的映射关系，特征矩阵中[j:j+n_dim]的列向量便是原始数据的降维表示。X_ndim = eig_vec_pickedreturn X_ndimif __name__ == '__main__':X = load_digits().datay = load_digits().targetdist = cal_pairwise_dist(X)max_dist = np.max(dist)print("max_dist", max_dist)X_ndim = le(X, n_neighbors = 20, t = max_dist*0.1)plt.scatter(X_ndim[:, 0], X_ndim[:, 1], c = y)plt.show()