1. 鄰接矩陣,度矩陣,拉普拉斯矩陣
給定一個無向圖:
我們可以用鄰接矩陣(Adjacent Matrix)表示它:
把這個鄰接矩陣記為W,W中的1表示有連接,0表示沒有連接,例如第一行第二列的1表示圖的節點1和節點2有連接,第一行第三列的0表示圖的節點1和節點3沒有連接在一起。因為是無向圖,所以W是對稱矩陣。
我們把W中的每一行相加,放到對角線上,然后把對角線以外的元素填0,得到度矩陣(Degree Matrix):
把度矩陣記為D,則
$\mathrm{D}_{i j}=\left\{\begin{array}{l}\sum_{j=1}^{n} w_{ij}, i=j \\ 0, i \neq j\end{array}\right.$
$D_{ii}$表示了第i個節點的度,反映了這個節點的重要性。
拉普拉斯矩陣(Laplacian Matrix)定義為L=D-W:
2. 把數據樣本構建成圖
給定樣本$x_1,x_2,...,x_n$,把這些樣本看作圖的一個個節點,那我們怎么把這些樣本節點用邊連接起來呢,哪些節點之間要有邊,哪些節點之間不需要有邊?因為拉普拉斯特征映射的核心思想是只保持局部的幾何結構不變,所以我們只需要把節點與其周圍若干個節點連接就可以,有下面2種較常用的選節點方法:
- $\varepsilon$-neighborhoods: 如果樣本$x_i$和$x_j$距離小於$\varepsilon$, 就把樣本$x_i$和$x_j$用邊連接起來。
- k-nearest neighbors: 樣本$x_i$只連接與它距離最近的k個樣本。
連接起來之后,我們需要給邊賦權重,有以下兩種常用方法:
- Heat Kernel:
$\mathrm{W}_{i j}=\left\{\begin{array}\\exp \left(-\frac{\left\|x_{i}-x_{j}\right\|_{2}^{2}}{t}\right), if\ x_i\ and\ x_j\ are\ connected \\ 0, otherwise \end{array}\right.$
- simple-minded: ($t \rightarrow \inf$)
$\mathrm{W}_{i j}=\left\{\begin{array} {}1, if\ x_i\ and\ x_j\ are\ connected \\ 0, otherwise \end{array}\right.$
在上一小節中的鄰接矩陣W用的是simple-minded方法。
3. The unnormalized graph Laplacian
關於性質1,李宏毅機器學習教學視頻semi-supervise視頻的50:16~58:36這段里有講到一個實際例子:
如果n=3,也就是只有3個節點,也可以直接把公式展開,各項消消掉,合並同類項等,就能比較清楚地知道上述的性質1是成立的:
4. The normalized graph Laplacian
5. The Algorithm
輸入樣本$x_1,x_2,...,x_n$,聚類數目k
- 根據第2節的方法構建鄰接矩陣W,根據第1小節的方法構建度矩陣D;
- 根據第3或4小節的方法構建拉普拉斯矩陣L;
- 求解廣義特征值問題:$L f=\lambda D f$;
- 把特征值從小到大排列為$\lambda_0,\lambda_1, \lambda_2...,\lambda_n$, 對應的特征向量為$f_0, f_1, f_2, ..., f_n$,其中最小的特征值$\lambda_0=0$,它對應的特征向量$f_0$是一個元素全為1的向量;
- 把$f_1, f_2, ..., f_k$這k個特征向量按列排列,形成大小為n*k的矩陣U;
- 把U的行記為$y_i (i=1,2,..,n)$,$y_i$就是樣本$x_i$從高維空間映射大k維空間后的坐標,即$y_i$是$x_i$的embedding。
- 對$y_i (i=1,2,..,n)$進行聚類,比如用k-means進行聚類。$y_i$的聚類結果就是$x_i$的聚類結果。
(注:1. 很多文獻中提到的相似度矩陣(similarity matrix)應該就是鄰接矩陣。2. 譜聚類算法有很多種變體,這里敘述的只是其中一種)
6. Why This Algorithm works
拉普拉斯特征映射只考慮保持局部特征映射,因此一個自然的想法是,原高維空間中距離較近的樣本$x_i$和$x_j$,映射為$y_i$和$y_j$,那么$y_i$和$y_j$的距離也應該比較接近,因此有人提出如下的目標函數:
$minimize\sum_{i=1}^{n} \sum_{j=1}^{n} (f_{i}-f_{j})^{2} W_{i j}$
如果$x_i$和$x_j$距離較近,則$W_{i j}$較大,這個時候如果$y_i$和$y_j$距離較遠,則目標函數就會很大,因此最小化目標函數可以讓$y_i$和$y_j$距離接近。
根據第三節的公式,有$\sum_{i=1}^{n} \sum_{j=1}^{n} (f_{i}-f_{j})^{2} W_{i j}=f^T W f$, 所以我們的目標函數可轉化為
$minimize f^T L f, s.t. f^T D f=1$
約束是為了防止向量f的任意縮放。
用拉格朗日乘子法把有約束優化問題轉化為無約束優化問題:
$minimize f^T L f + \lambda f^T D f$
函數對f求偏導,令該偏導為0,得到$(L-\lambda D)f=0$,即$Lf=\lambda D f$
7. 實例
下面的程序來自:https://github.com/heucoder/dimensionality_reduction_alo_codes/tree/master/codes/LE
該程序中有兩個laplacian eigenmaps的例子,第一個是把三維的swiss roll拉平成2維,第二個例子是把64維的手寫數字樣本降成2維。
import numpy as np import matplotlib.pyplot as plt from sklearn.datasets import load_digits from mpl_toolkits.mplot3d import Axes3D ''' author: heucoder email: 812860165@qq.com date: 2019.6.13 ''' 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) # np.square(x)把矩陣x中的每個元素都各自平方 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): ''' :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] # 舍棄小於1e-6的特征值對應的特征向量 eig_val_picked = eig_val[j:j+n_dims] print(eig_val_picked) eig_vec_picked = eig_vec[:, sort_index_] X_ndim = eig_vec_picked return X_ndim if __name__ == '__main__': # example 1: swiss roll 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() # example 2: hand-written digits # 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()
運行結果:
swiss roll:
digits:
參考資料:
[1] 從拉普拉斯矩陣說到譜聚類
[3] A Tutorial on Spectral Clustering
[4] Laplacian Eigenmaps for Dimensionality Reduction and Data Representation
[5] CMU機器學習課程
[6] Laplacian Eigenmap中的核心推導過程
[7] 漫談 Clustering (4): Spectral Clustering
[8] 譜聚類(spectral clustering)原理總結 (很詳細)