拉普拉斯特征映射(Laplacian Eigenmaps)

1 介绍

  拉普拉斯特征映射(Laplacian Eigenmaps)是一种不太常见的降维算法,它看问题的角度和常见的降维算法不太相同,是从局部的角度去构建数据之间的关系。具体来讲,拉普拉斯特征映射是一种基于图的降维算法,它希望相互间有关系的点(在图中相连的点)在降维后的空间中尽可能的靠近,从而在降维后仍能保持原有的数据结构。

2 推导

  拉普拉斯特征映射通过构建邻接矩阵为 $W$ 的图来重构数据流形的局部结构特征。

  其主要思想是:如果两个数据实例 $i$ 和 $j$ 很相似,那么 $i$ 和 $j$ 在降维后目标子空间中应该尽量接近。设数据实例的数目为 $n$ ,目标子空间即最终的降维目标的维度为 $m$ 。 定义 $ n \times m$ 大小的矩阵 $Y$ ,其中每一个行向量 $y_{i}^{T}$ 是数据实例 $i$ 在目标 $m$ 维子空间中的向量表示(即降维后的数据实例 $i$ )。我们的目的是让相似的数据样例 $i$ 和 $j$ 在降维后的目标子空间里仍旧尽量接近,故拉普拉斯特征映射优化的目标函数如下:

    $\min \sum\limits _{i, j}\left\|y_{i}-y_{j}\right\|^{2} W_{i j}$

  下面开始推导:

    $ \begin{array}{l} \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n}&\left\|y_{i}-y_{j}\right\|^{2} W_{i j} \\ &=\sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n}\left(y_{i}^{T} y_{i}-2 y_{i}^{T} y_{j}+y_{j}^{T} y_{j}\right) W_{i j} \\ &=\sum\limits_{i=1}^{n}\left(\sum\limits_{j=1}^{n} W_{i j}\right) y_{i}^{T} y_{i}+\sum\limits_{j=1}^{n}\left(\sum\limits_{i=1}^{n} W_{i j}\right) y_{j}^{T} y_{j}-2 \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} y_{i}^{T} y_{j} W_{i j} \\ &=2 \sum\limits_{i=1}^{n} D_{i i} y_{i}^{T} y_{i}-2 \sum\limits_{i=1}^{n} \sum\limits_{j=1}^{n} y_{i}^{T} y_{j} W_{i j} \\ &=2 \sum\limits_{i=1}^{n}\left(\sqrt{D_{i i}} y_{i}\right)^{T}\left(\sqrt{D_{i i}} y_{i}\right)-2 \sum\limits_{i=1}^{n} y_{i}^{T}\left(\sum\limits_{j=1}^{n} y_{j} W i j\right) \\ &=2 \operatorname{trace}\left(Y^{T} D Y\right)-2 \sum\limits_{i=1}^{n} y_{i}^{T}(Y W)_{i} \\ &=2 \operatorname{trace}\left(Y^{T} D Y\right)-2 \operatorname{trace}\left(Y^{T} W Y\right) \\ &=2 \operatorname{trace}\left[Y^{T}(D-W) Y\right] \\ &=2 \operatorname{trace}\left(Y^{T} L Y\right) \end{array} $

  其中 $W $ 是图的邻接矩阵,对角矩阵 $D$ 是图的度矩阵 $\left(D_{i i}=\sum\limits_{j=1}^{n} W_{i j}\right)$ ,$ L=D-W$ 为图的拉普拉斯矩阵。

  变换后的拉普拉斯特征映射优化的目标函数如下:

    $\begin{array}{l}\min \operatorname{trace}\left(Y^{T} L Y\right)\\ \text { s.t. } Y^{T} D Y=I \end{array}$

  限制条件  $s . t . Y^{T} D Y=I$ 为了消除低维空间中的缩放因子,也为了保证 $D_{ii}$ 值较大的样本点在低维空间中更为重要,下面用拉格朗日乘子法对目标函数求解:

    $f(Y)=\operatorname{tr}\left(Y^{T} L Y\right)+\operatorname{tr}\left[\Lambda\left(Y^{T} D Y-I\right)\right]$

    $\begin{array}{l} \frac{\partial f(Y)}{\partial Y}&=L Y+L^{T} Y+D^{T} Y \Lambda^{T}+D Y \Lambda \\ &=2 L Y+2 D Y \Lambda=0 \end{array}$

    $\therefore L Y=-D Y \Lambda$

  其中,$\Lambda$ 为一个对角矩阵,另外 $L$ 、 $D$ 均为实对称矩阵,其转置与自身相等。对于单独的 $y$ 向量,上式可写为: $L y=\lambda D y$,这是一个广义特征值问题。通过求得 $m$ 个最小非零特征值所对应的特征向量,即可达到降维的目 的。

  关于这里为什么要选择 $m$ 个最小非零特征值所对应的特征向量。将 $L Y=-D Y \Lambda $ 带回到 $\min \operatorname{trace}\left(Y^{T} L Y\right)$ 中,由于有着约束条件 $Y^{T} D Y=I$ 的限制,可以得到 $ \min \quad \operatorname{trace}\left(Y^{T} L Y\right)=\min \quad t r a c e(-\Lambda)$ 。即为特 征值之和。我们为了目标函数最小化,要选择最小的 $m$ 个特征值所对应的特征向量。

3 步骤

  使用时算法具体步骤为:

  步骤1:构建图

    使用某一种方法来将所有的点构建成一个图,例如使用KNN算法,将每个点最近的K个点连上边。K是一个预先设定的值。

  步骤2:确定权重

    确定点与点之间的权重大小,例如选用热核函数来确定,如果点 i 和点 j 相连,那么它们关系的权重设定为:

    $W_{i j}=e^{-\frac{\left\|x_{i}-x_{j}\right\|^{2}}{t}}$

    另外一种可选的简化设定是 $W_{i j}=1$ 如果点 $i$ ,$ j$ 相连,否则 $W_{i j}=0 $ 。

  步骤3:特征映射

    计算拉普拉斯矩阵 $L$ 的特征向量与特征值: $L y=\lambda D y $

    使用最小的 $m$ 个非零特征值对应的特征向量作为降维后的结果输出。

 

Code:

# coding:utf-8

import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import load_digits from mpl_toolkits.mplot3d import Axes3D def make_swiss_roll(n_samples=100, noise=0.0, random_state=None): #Generate a swiss roll dataset.
    t = 1.5 * np.pi * (1 + 2 * np.random.rand(1, n_samples)) x = t * np.cos(t) y = 83 * np.random.rand(1, n_samples) z = t * np.sin(t) X = np.concatenate((x, y, z)) X += noise * np.random.randn(3, n_samples) X = X.T t = np.squeeze(t) return X, t def 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 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 = 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 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_] print("eig_val[:10]: ", eig_val[:10]) j = 0 while eig_val[j] < 1e-6: j+=1

    print("j: ", j) sort_index_ = sort_index_[j:j+n_dims] 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)) X_ndim = eig_vec_picked return X_ndim if __name__ == '__main__': X, Y = make_swiss_roll(n_samples = 2000)  #生成瑞士卷数据集
    X_ndim = le(X, n_neighbors = 5, t = 20) fig = plt.figure(figsize=(12,6)) ax1 = fig.add_subplot(121, projection='3d') ax1.scatter(X[:, 0], X[:, 1], X[:, 2], c = Y) ax2 = fig.add_subplot(122) ax2.scatter(X_ndim[:, 0], X_ndim[:, 1], c = Y) plt.show() X = load_digits().data y = load_digits().target dist = 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.savefig("LE2.png") plt.show()

 

 2 total views,  1 views today

页面下部广告