降维和聚类系列(二)：拉普拉斯特征映射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{egin{array}{l}sum_{j=1}^{n} w_{ij}, i=j \ 0, i eq jend{array} ight.$

$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{egin{array}\exp left(-frac{left|x_{i}-x_{j} ight|_{2}^{2}}{t} ight), if x_i and x_j are connected \ 0, otherwise end{array} ight.$

• simple-minded: ($t ightarrow inf$)

$mathrm{W}_{i j}=left{egin{array} {}1, if x_i and x_j are connected \ 0, otherwise end{array} ight.$ 

在上一小节中的邻接矩阵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$的距离也应该比较接近，因此有人提出如下的目标函数：

$minimizesum_{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] 从拉普拉斯矩阵说到谱聚类

[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）原理总结 (很详细)