1. ISOMAP建模

MDS的基本思想是在降维后的低维空间中保持原始空间中样本之间的距离。低维空间中的距离使用欧氏距离衡量；通常原始空间中的距离也用欧氏距离描述，但是如果原始空间中的数据分布在非线性的流形(manifold)上，则欧氏距离无法反映两点之间的真实距离。

如下图(A)所示的数据分布在空间中的一个“瑞士卷”流形上，图中两点之间的欧氏距离用蓝色虚线表示；两点间的实际距离为测地线(geodesic)距离，用蓝色实线表示。从图中可以看出，流形空间中的直线距离是不可达的。

通常流形上的测地线距离不好直接计算，但是流形在局部与欧氏空间同胚，这意味着流形在局部具有欧氏空间的性质。对于流形空间中的每个点，可以根据欧氏距离计算其近邻点，进而建立一个近邻连接图。图中近邻点之间存在连接，非近邻点之间则不存在连接。因此计算两点之间的测地线距离问题，就转换成计算近邻连接图上对应两点之间的最短路径问题。

如图(B)所示采用近邻距离代替测地距离，近邻点之间使用欧氏距离度量，近似路径用红色实线表示。将该流形展开后得到图(C)，可以看出近邻距离可以近似低维流形上的测地线距离。

近邻图的构建有两种方法。可以选择与每个点欧氏距离最近的$k$个点作为邻近点，此时得到的近邻图称为$k$近邻图；也可以选择与每个点欧氏距离小于$\epsilon$的点作为邻近点，此时得到的近邻图称为$\epsilon$近邻图。

如果近邻范围设置得比较大，则可能将欧氏距离很近但测地距离很远的点划分为近邻点，造成“短路”问题；如果近邻范围设置得比较小，则可能有些区域与其他区域不存在连接，造成“断路”问题。

2. 求解ISOMAP

ISOMAP算法的一般流程如下：

给定输入样本$X \in \Bbb{R}^{N \times d}$、近邻参数$n$和降维维度$k$;
对于每个样本点计算与其距离最近的$n$个样本点(不包括自身)的距离，其他距离视为无穷大；
使用最短路径算法更新每两个样本点之间的最短距离；
对上述距离矩阵应用MDS算法，得到降维结果。

在计算最短路径时，可以使用Dijkstra算法或Floyd-Warshall算法。下面介绍Floyd-Warshall算法。

计算任意两点之间的最短路径，如果两点$i$和$j$之间存在一点$k$，使得$d_{ij}>d_{ik}+d_{kj}$，则更新两点之间的距离。最终得到重新填充的距离矩阵$D_1$。

def Floyd(D, n_neighbors = 10):
    N = D.shape[0]
    # 填充距离矩阵
    Max = np.max(D)*1000
    D1 = np.ones([N,N])*Max
    index = np.argsort(D, axis=1)
    for i in range(N):
        D1[i, index[i,0:n_neighbors+1]]=D[i, index[i,0:n_neighbors+1]]
    # 更新距离矩阵
    for k in range(N):
        for i in range(N):
            for j in range(N):
                if D1[i,k]+D1[k,j]<D1[i,j]:
                    D1[i,j] = D1[i,k]+D1[k,j]
    return D1

3. 实现ISOMAP

① ISOMAP from scratch

由上述介绍，ISOMAP的一般步骤如下：

给定输入样本$X \in \Bbb{R}^{N \times d}$、近邻参数$n$和降维维度$k$;
对于每个样本点计算与其距离最近的$n$个样本点(不包括自身)的距离，其他距离视为无穷大；
使用最短路径算法更新每两个样本点之间的最短距离；
对上述距离矩阵应用MDS算法，得到降维结果。

def MDS(D, k):
    N = D.shape[0]
    # 计算内积矩阵B
    T1 = np.sum(D**2, axis=0, keepdims=True)/N
    T2 = np.sum(D**2, axis=1, keepdims=True)/N
    T3 = np.sum(D**2, keepdims=True)/N**2
    B = (T1+T2-T3-D**2)/2
    # 计算特征值和特征向量
    eig_values, eig_vectors = np.linalg.eig(B)
    # 选择前k个最大的特征值标号
    index = np.argsort(-eig_values)[:k]
    # 选择对应的特征值和特征向量
    MDS_values = eig_values[index]
    MDS_vectors = eig_vectors[:, index] # (N, k)
    # 降维
    reduced_data = np.diagflat(MDS_values**0.5).dot(MDS_vectors.T) # (k, N)
    return reduced_data

def ISOMAP(data, k):
    N, d = data.shape
    D = np.zeros([N, N])
    # 计算距离矩阵D
    for i in range(N):
        for j in range(N):
            D[i, j] = np.sqrt(np.sum((data[i]-data[j])**2))
    D1 = Floyd(D, n_neighbors = 10)
    reduced_data = MDS(D1.T, k=2)
    return reduced_data.T

② ISOMAP from sklearn

ISOMAP也可以通过sklearn库快速实现：

from sklearn.datasets import make_s_curve
from sklearn.manifold import Isomap
import matplotlib.pyplot as plt

data, color = make_s_curve(n_samples = 500,
                           noise = 0.1,
                           random_state = 42)
X_isomap = Isomap(n_neighbors = 10, n_components = 2).fit_transform(data)

plt.figure(figsize=(10, 5))
plt.scatter(X_isomap[:, 0], X_isomap[:, 1], c=color, label="ISOMAP")
plt.legend()
plt.show()