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$ 保證優化問題有解，下面用拉格朗日乘子法對目標函式求解:
$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()