(五)等度量映射（Isomap）

一、问题

低维流形嵌入到高维空间之后，直接再高维空间中计算直线距离具有误导性，因为高维空间中的直线距离在低维嵌入流形上是不可达的。低维嵌入流形上两点间的距离是“测地线”距离，是两点之间的本真距离。直接在高维空间中计算直线距离是不恰当的。

二、解决

如何计算测地线距离呢？流形在局部上与欧氏空间同胚，对于每个点基于欧氏距离找出其近邻点，建立一个近邻连接图(1指定近邻点的个数，k近邻图；2指定距离的阈值，$ϵ$近邻图)，即解决近邻连接图上最短路径的问题（使用Dijkstra或者Floyd算法求最短路径）。由此得到距离矩阵，接下来使用MDS方法来获得样本点在低维空间中的坐标。MDS原理及python实现

三、算法描述

输入：样本集$D=(x_{(1)},x_{(2)},...x_{(n)}{)}^T;$
    \qquad \,\,\, 近邻参数$k$;
    \qquad \,\,\, 低维空间维数$d^{{}^{\prime }}$
输出：样本集$D$在低维空间中的投影$Z=(z_{(1)},z_{(2)},...z_{(n)}{)}^T;$
过程：
for $i=1,2,..n$ do
\qquad 确定$x_{(i)}^T$的$k$近邻，$x_{(i)}^T$与$k$近邻之间的距离设置为欧氏距离，与其他点之间的距离设置为无穷大；
end for
调用最短路径算法计算任意两样本点之间的距离，得到距离矩阵，作为MDS算法的输入；
return      \,\,\,\, MDS算法的输出

四、python实现

# coding:utf-8
import numpy as np
from sklearn.datasets import make_s_curve
import matplotlib.pyplot as plt
from sklearn.manifold import Isomap
from mpl_toolkits.mplot3d import Axes3D

def floyd(D,n_neighbors=15):    #k近邻
    Max = np.max(D)*1000
    n1,n2 = D.shape
    k = n_neighbors
    D1 = np.ones((n1,n1))*Max
    D_arg = np.argsort(D,axis=1) #返回从小到大的索引值，每一列进行排序
    print(D_arg)
    for i in range(n1):
        D1[i,D_arg[i,0:k+1]] = D[i,D_arg[i,0:k+1]]  #找出与i最近的k个数
    for k in range(n1):
        for i in range(n1):
            for j in range(n1):
                if D1[i,k]+D1[k,j]<D1[i,j]:
                    D1[i,j] = D1[i,k]+D1[k,j]  #需再次理解 
    return D1

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)   #square计算各元素平方      1 按行相加  0按列相加
    dist = np.add(np.add(-2 * np.dot(x, x.T), sum_x).T, sum_x)
    #返回任意两个点之间距离的平方
    return dist

def my_mds(dist, n_dims):    #距离矩阵分解，得到降维之后的数据
    # dist (n_samples, n_samples)
    dist = dist**2
    n = dist.shape[0]
    T1 = np.ones((n,n))*np.sum(dist)/n**2
    T2 = np.sum(dist, axis = 1)/n
    T3 = np.sum(dist, axis = 0)/n

    B = -(T1 - T2 - T3 + dist)/2

    eig_val, eig_vector = np.linalg.eig(B)
    index_ = np.argsort(-eig_val)[:n_dims]
    picked_eig_val = eig_val[index_].real
    picked_eig_vector = eig_vector[:, index_]

    return picked_eig_vector*picked_eig_val**(0.5)

def my_Isomap(data,n=2,n_neighbors=30):
    D = cal_pairwise_dist(data)
    D[D < 0] = 0
    D = D**0.5       #距离矩阵
    D_floyd=floyd(D, n_neighbors)     #k近邻之后得到的距离矩阵
    data_n = my_mds(D_floyd, n_dims=n)
    return data_n

def scatter_3d(X, y):
    fig = plt.figure(figsize=(6, 5))
    ax = fig.add_subplot(111, projection='3d')
    ax.scatter(X[:, 0], X[:, 1], X[:, 2], c=y, cmap=plt.cm.hot)
    ax.view_init(10, -70)
    ax.set_xlabel("$x_1$", fontsize=18)
    ax.set_ylabel("$x_2$", fontsize=18)
    ax.set_zlabel("$x_3$", fontsize=18)
    plt.show()

if __name__ == '__main__':
    X, Y = make_s_curve(n_samples = 500,
                           noise = 0.1,
                           random_state = 42)  #生成S型曲线数据集

    data_1 = my_Isomap(X, 2, 10)    

    data_2 = Isomap(n_neighbors = 10, n_components = 2).fit_transform(X)

    plt.figure(figsize=(8,4))
    plt.subplot(121)
    plt.title("my_Isomap")
    plt.scatter(data_1[:, 0], data_1[:, 1], c = Y)

    plt.subplot(122)
    plt.title("sklearn_Isomap")
    plt.scatter(data_2[:, 0], data_2[:, 1], c = Y)
    plt.savefig("Isomap1.png")
    plt.show()