降維和聚類系列(二)：拉普拉斯特征映射Laplacian Eigenmaps，譜聚類，實例代碼

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

1. 根據第2節的方法構建鄰接矩陣W，根據第1小節的方法構建度矩陣D；
2. 根據第3或4小節的方法構建拉普拉斯矩陣L；
3. 求解廣義特征值問題：$L f=\lambda D f$；
4. 把特征值從小到大排列為$\lambda_0,\lambda_1, \lambda_2...,\lambda_n$, 對應的特征向量為​$f_0, f_1, f_2, ..., f_n$，其中最小的特征值$\lambda_0=0$，它對應的特征向量$f_0$是一個元素全為1的向量；
5. 把​$f_1, f_2, ..., f_k$這k個特征向量按列排列，形成大小為n*k的矩陣U；
6. 把U的行記為$y_i (i=1,2,..,n)$，$y_i$就是樣本$x_i$從高維空間映射大k維空間后的坐標，即$y_i$是$x_i$的embedding。
7. 對$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: [e-mail]
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] 從拉普拉斯矩陣說到譜聚類

[2] Laplacian Eigenmaps

[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）原理總結 (很詳細)