1、介绍 先看一些定义 ? 流形是一个具有欧几里得空间性质的局部空间，可以用欧氏距离计算。 ? 同一胚胎，在拓扑学中，如果两个流形可以通过弯曲、延伸和剪切（只要它们最终完全沿着原始剪切的间隙重新粘贴）将其中一个变成另一个，则认为两者是同一胚胎。因此，也可以说，流形是与欧洲空间局部同一胚胎的空间。 ? 测量地线距离，空间中两点的最短或最长路径 
如图所示:蓝色实线表示测地距离，蓝色虚线是我们传统上所知道的欧氏距离 根据欧氏距离找出每个点的近邻点，然后建立近邻连接图。图中的近邻点之间有连接，而不是近邻点之间没有连接。因此，计算两点之间地线距离的问题转化为计算近邻连接图中两点之间最短路径的问题. 在近邻连接图上计算两点间的最短路径，可采用着名的Dikstra算法或Floyd 算法，在得到任意两点的距离后，可以通过 MDS 在低维空间中获得样本点坐标的方法 ISOMAP，使用测地距离使高维空间中任意点对点之间的测地距离关系在低维空间中保持不变，这也是其主要优势，而不是使用原始欧几里得距离，可以更好地控制数据信息的丢失，在低维空间中更全面地显示高维空间的数据 Isomap算法是在MDS以算法为基础的算法，MDS算法是在降维后保持样本间距不变，而ISOMAP还需要实现过程MDS加入算法实现 算法描述过程如下：
2、算法实现 2.1数据的选择 先说数据。这里我的数据是随机生成的。 为什么别无选择？doc2vec文档处理后的数据呢？ 因为这样生成的数据确实是高维的，不使用并不意味着ISOMAP不能处理，因为我解释了ISOMAP该模型的目的是展示该算法的实现过程及其优点。在解释过程中，我们希望在图像的帮助下更方便地显示实现效果。由于高维空间的我们无法想象的，我需要一个三维显示空间来将维度降低到二维，让大家更容易理解过程 2、算法实现 2.1数据的选择 先说数据。这里我的数据是随机生成的。 为什么别无选择？doc2vec文档处理后的数据呢？ 因为这样生成的数据确实是高维的，不使用并不意味着ISOMAP不能处理，因为我解释了ISOMAP该模型的目的是展示该算法的实现过程及其优点。在解释过程中，我们希望在图像的帮助下更方便地显示实现效果。由于高维空间的我们无法想象的，我需要一个三维显示空间来将维度降低到二维，让大家更容易理解过程 2.主要实现过程 ? 生成随机数据 借助于scikit-learn 的 datasets 我们专门使用模块make_s_curve生成流行数据，突出ISOMAP优势
X, Y = make_s_curve(n_samples = 500, #主要数据是X，Y它只是用来生成不同的颜色 noise = 0.1, random_state = 42) X在这里生成的是500*3.每列数据作为三维空间的坐标，相当于我们现在准备的500点。y事实上，这只是我们用来显示不同颜色图形的代表数。在这里，为了显示最初的效果，我们绘制了三维空间图像来显示数据分布
def scatter_3d(X, y): fig = plt.figure() ax = plt.axes(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$", fontsizespan class="token operator">=18)
ax.set_ylabel("$x_2$", fontsize=18)
ax.set_zlabel("$x_3$", fontsize=18)
plt.show(block=False)
结果如下： 生成距离矩阵 上面我们生成了500个点，现在我们先计算两两之间的欧氏距离
1. def cal_pairwise_dist(x):
2. N,D = np.shape(x)
3. dist = np.zeros([N,N])
4. for i in range(N):
5. for j in range(N):
6. # dist[i,j] = np.dot((x[i]-x[j]),(x[i]-x[j]).T)
7. dist[i,j] = np.sqrt(np.dot((x[i]-x[j]),(x[i]-x[j]).T))
8. # dist[i,j] = np.sum(np.abs(x[i]-x[j]))
9.
10. #返回任意两个点之间距离
11. return dist
得到距离矩阵 MDS降维实现 为了凸显ISOMAP的效果，我们这里先实现一下MDS降维处理后的效果，我们可以对比结果 使用MDS算法将数据降维到二维空间
1. def my_MDS(dist, n_dims):
2. # dist (n_samples, n_samples)
3. dist = dist**2
4. n = dist.shape[0]
5. T1 = np.ones((n,n))*np.sum(dist)/n**2
6. T2 = np.sum(dist, axis = 1)/n
7. T3 = np.sum(dist, axis = 0)/n
8.
9. B = -(T1 - T2 - T3 + dist)/2
10.
11. eig_val, eig_vector = np.linalg.eig(B)
12. index_ = np.argsort(-eig_val)[:n_dims]
13. picked_eig_val = eig_val[index_].real
14. picked_eig_vector = eig_vector[:, index_]
15.
16. return picked_eig_vector*picked_eig_val**(0.5)
ISOMAP的实现 我们这里的主要部分就是构建最短路径矩阵，思想上面也说了，我们需要规定临界点的个数k，保留下每个点的k个临界点的欧式距离值，剩余其他的所有距离看作是无穷大，利用相邻点欧氏距离的累加，逐步寻找跨越的点的距离，这个就是测地距离
# 构建最短路径图
def floyd(D,n_neighbors=15):
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)
for i in range(n1):
D1[i,D_arg[i,0:k+1]] = D[i,D_arg[i,0:k+1]]
for k in tqdm(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
计算完距离之后我们得到一个测地距离矩阵 接下来执行我们上面算法描述中的第5步使用MDS对测地距离矩阵进行降维，依旧是降维到二维空间，绘制图像得到我们ISOMAP处理后的结果
其实，sklearn中由关于ISOMAP算法的包，我们可以直接调用，这样比较方便
1. #直接调包
2. data_ISOMAP2 = Isomap(n_neighbors = 10, n_components = 2).fit_transform(X)
这样的结果是 ISOMAP处理后的展现的图像，保留了原来数据的相对位置，他其实是将三维空间的图像展开类似于展开一样展现在二维中，可以看到我门调包处理后的结果与我们自己写的算法实现结果大致相同，说明我们自己手动实现的效果也不错。 当然ISOMAP的作用不仅仅限于这些，我们可以用ISOMAP处理更高维的数据而且保持测地数据不变。
全部代码：
# -*- coding: utf-8 -*-
""" Created on Fri Nov 12 17:52:03 2021 @author: Yangz """
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import make_s_curve
from sklearn.manifold import Isomap
from tqdm import tqdm
from mpl_toolkits.mplot3d import Axes3D
# x 维度 [N,D]
def cal_pairwise_dist(x):
N,D = np.shape(x)
dist = np.zeros([N,N])
for i in range(N):
for j in range(N):
# dist[i,j] = np.dot((x[i]-x[j]),(x[i]-x[j]).T)
dist[i,j] = np.sqrt(np.dot((x[i]-x[j]),(x[i]-x[j]).T)) #欧氏距离的计算 相当于 sqrt((x1-x2)^2+(y1-y2)^2+(z1-z2)^2)
# dist[i,j] = np.sum(np.abs(x[i]-x[j]))
#返回任意两个点之间距离
return dist
# 构建最短路径图
def floyd(D,n_neighbors):#=15
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)
for i in range(n1):
D1[i,D_arg[i,0:k+1]] = D[i,D_arg[i,0:k+1]]
for k in tqdm(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 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 # 计算dist(ij)^2
T2 = np.sum(dist, axis = 1)/n # 计算dist(i.)^2
T3 = np.sum(dist, axis = 0)/n # 计算dist(.j)^2
B = -(T1 - T2 - T3 + dist)/2 #B = Z^T * Z
#对矩阵B做特征分解
eig_val, eig_vector = np.linalg.eig(B) #eig_val矩阵B的特征值，eig_vector归一化的特征向量
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(D,n=2,n_neighbors=30):
D_floyd=floyd(D, n_neighbors)
data_n = my_mds(D_floyd, n_dims=n)
return data_n
def scatter_3d(X, y):
fig = plt.figure()
ax = plt.axes(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(block=False)
if __name__ == '__main__':
X, Y = make_s_curve(n_samples = 500, #主要数据是X，Y只是用来生成不同的颜色
noise = 0.1,
random_state = 42)
scatter_3d(X,Y)
# 计算欧式距离，得到距离矩阵
dist = cal_pairwise_dist(X)
# MDS 降维
# data_MDS = my_mds(dist, 2)
#
# plt.figure()
# plt.title("my_MSD")
# plt.scatter(data_MDS[:, 0], data_MDS[:, 1], c = Y)
# plt.show(block=False)
#
#
# ISOMAP 降维
data_ISOMAP = my_Isomap(dist, 2, 10)
plt.figure()
plt.title("my_Isomap")
plt.scatter(data_ISOMAP[:, 0], data_ISOMAP[:, 1], c = Y)
plt.show(block=False)
#直接调包
data_ISOMAP2 = Isomap(n_neighbors = 10, n_components = 2).fit_transform(X)
plt.figure()
plt.title("sk_Isomap")
plt.scatter(data_ISOMAP2[:, 0], data_ISOMAP2[:, 1], c = Y)
plt.show(block=False)
plt.show()