哈喽，大家好~
今儿和大家来聊聊降维技术中，非常重要的一环，Isomap。
大家可以想象，在一张弯曲的纸上画了一条曲线，再把纸折叠成一个“瑞士卷”（Swiss roll）。
在三维空间里，两点看起来可能很近（欧氏距离小），但沿着纸面的真实“道路”却很远。
那么，Isomap 的目标就是：利用样本点之间沿流形的“测地线距离”，把数据嵌入到低维空间，同时尽量保持这些测地线距离。这样即便原始数据在高维空间中弯曲复杂，Isomap 也能把内在的低维结构“铺平”呈现出来。
核心逻辑
局部近邻假设：流形是“局部平坦”的。
也就是说，任意一个点附近的一小块区域，点与点之间的直线距离近似测地线距离。于是我们可以用 k 近邻或 ε-邻域把局部连接建立起来。
建邻接图：把每对近邻点用一条边连接，边权等于这两点的欧氏距离。这样局部距离保留，但全局结构不直接体现。
近似测地线距离：在邻接图上计算最短路径距离（例如 Dijkstra 或 Floyd），得到点对间通过图的最短路径长度，作为测地线距离的近似值。这样就能把“沿流形的远近”信息捕捉到。
MDS（经典多维尺度）把距离矩阵转换为内积矩阵，然后做特征分解（eigendecomposition），取前 d 个特征向量按特征值的平方根缩放，得到 d 维坐标。
数学上：
-
给定测地线距离矩阵 ，先取平方得 ；
-
做中心化：令 ；
-
计算内积矩阵（Gram 矩阵）：
-
特征分解 ，按降序取前 个正特征值与对应特征向量，低维嵌入为
这就是 Isomap 的数学核心。
参数与注意事项：
-
k（邻居数）或 ε：太小会把图切分成多个连通分量（无法估计某些点间的测地线）；太大则会把非流形上的“捷径”连上，破坏拓扑。通常选择一个既保证连通又保持局部性的 k（例如 5~15）；
-
数据噪声：邻接图可能会连上噪声点，影响测地线估计；
-
计算复杂度：计算全对距离和求最短路径对大样本比较耗时，可以用近似邻居或稀疏化策略。
完整案例
我们使用一个常见的流形数据集 Swiss roll（瑞士卷）做示例，展示整个 Isomap 流程~
import torch
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.csgraph import shortest_path
from scipy.spatial.distance import pdist, squareform
seed = 42
torch.manual_seed(seed)
np.random.seed(seed)
# 1）生成 Swiss Roll 数据
def make_swiss_roll(n_samples=2000, noise=0.0, random_state=None):
# 生成 t（一维参数）和高度 h
if random_state is not None:
np.random.seed(random_state)
t = 1.5 * np.pi * (1 + 2 * np.random.rand(n_samples))
h = 21 * np.random.rand(n_samples)
x = t * np.cos(t)
z = t * np.sin(t)
X = np.vstack((x, h, z)).T
X = X + noise * np.random.randn(*X.shape)
return torch.from_numpy(X).float(), t # 返回 torch 张量与 t
X, t = make_swiss_roll(n_samples=2000, noise=0.5, random_state=0)
# 2）构造 k-NN 邻接图（使用欧氏距离）
def construct_knn_graph(X_tensor, k=10, symmetric=True):
# X_tensor: (n, d)
n = X_tensor.shape[0]
# 计算成对距离（使用 PyTorch）
with torch.no_grad():
D = torch.cdist(X_tensor, X_tensor, p=2) # shape (n, n)
D_np = D.numpy()
# 初始化带无穷的矩阵
INF = 1e12
G = np.full((n, n), INF, dtype=float)
for i in range(n):
idx = np.argsort(D_np[i]) # 包含自己（0）
neighbors = idx[1:k+1] # 排除自己
G[i, neighbors] = D_np[i, neighbors]
if symmetric:
# 对称化（取最小或平均），这里取最小
G = np.minimum(G, G.T)
# 保证对角为 0
np.fill_diagonal(G, 0.0)
return G, D_np
k = 10
G, D_euclid = construct_knn_graph(X, k=k, symmetric=True)
# 3）在图上计算最短路径距离（测地线近似）
dist_matrix = shortest_path(G, method='D', directed=False, unweighted=False)
# 若图不连通，某些距离会是 inf；需要处理（这里我们假设连通）
# 4）经典 MDS（用测地线距离矩阵做 MDS）
def classical_mds_from_dist(D, n_components=2):
# D: (n, n) 距离矩阵（非负），可以是测地线距离
n = D.shape[0]
# 如果存在 inf，需要先处理（这里简单替换为很大值）
D = np.array(D, dtype=float)
inf_mask = np.isinf(D)
if inf_mask.any():
# 使用一个很大的数替代（或更好地先保证图连通）
D[inf_mask] = np.nanmax(D[~inf_mask]) * 10.0
# 1) 平方化
D2 = D ** 2
# 2) 中心化矩阵 H
J = np.eye(n) - np.ones((n, n)) / n
# 3) Gram 矩阵 B
B = -0.5 * J.dot(D2).dot(J)
# 4) 特征分解
# 使用 numpy eigh（B 是对称）
eigvals, eigvecs = np.linalg.eigh(B)
# 降序排序
idx = np.argsort(eigvals)[::-1]
eigvals = eigvals[idx]
eigvecs = eigvecs[:, idx]
# 保留正特征值
positive = eigvals > 1e-10
L = np.diag(np.sqrt(eigvals[positive][:n_components]))
V = eigvecs[:, positive][:, :n_components]
Y = V.dot(L)
return Y, eigvals
Y_2d, eigvals = classical_mds_from_dist(dist_matrix, n_components=2)
# 5）计算 residual variance（残差方差），用于评价低维表示保留原始距离信息的程度
def residual_variance(D_orig, D_embedded):
# 计算 1 - R^2，其中 R 是距阵向量化的相关系数
# D_orig, D_embedded: 对称距离矩阵
v1 = D_orig[np.triu_indices_from(D_orig, k=1)]
v2 = D_embedded[np.triu_indices_from(D_embedded, k=1)]
# 只保留有限值
mask = np.isfinite(v1) & np.isfinite(v2)
v1 = v1[mask]
v2 = v2[mask]
# 相关系数 R
R = np.corrcoef(v1, v2)[0, 1]
return 1 - R**2
# 计算 2D 嵌入后的距离矩阵
D_emb_2d = squareform(pdist(Y_2d))
res_var_2d = residual_variance(dist_matrix, D_emb_2d)
print("Residual variance (2D):", res_var_2d)
# 6）绘图：准备多张图来分析 Isomap 的效果
import matplotlib as mpl
mpl.rcParams['figure.dpi'] = 120
# 图 1：原始瑞士卷（3D），按参数 t 着色
from mpl_toolkits.mplot3d import Axes3D
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection='3d')
p = ax.scatter(X[:,0].numpy(), X[:,1].numpy(), X[:,2].numpy(),
c=t, cmap='rainbow', s=8)
ax.set_title("Swiss Roll 原始数据（按 t 着色）")
plt.colorbar(p, ax=ax, shrink=0.6, label='t parameter')
plt.tight_layout()
plt.show()
# 图 2：邻接图（只画一部分点的边以免图太乱），在 3D 空间中显示 k-NN 边
def plot_knn_graph_3d(X_np, G, sample_for_edges=500):
n = X_np.shape[0]
idx_sample = np.random.choice(n, size=sample_for_edges, replace=False)
pts = X_np[idx_sample]
fig = plt.figure(figsize=(8, 5))
ax = fig.add_subplot(111, projection='3d')
ax.scatter(X_np[:,0], X_np[:,1], X_np[:,2], c='lightgray', s=6, alpha=0.6)
# draw edges for sampled nodes
for i, ii in enumerate(idx_sample):
neighbors = np.where(G[ii] < 1e11)[0]
for j in neighbors:
if j <= ii: continue # 避免重复画线
ax.plot([X_np[ii,0], X_np[j,0]],
[X_np[ii,1], X_np[j,1]],
[X_np[ii,2], X_np[j,2]],
c=plt.cm.viridis((i % 256)/256), alpha=0.6, linewidth=0.6)
ax.set_title(f"k-NN Graph (k={k}) 部分边展示")
plt.tight_layout()
plt.show()
plot_knn_graph_3d(X.numpy(), G, sample_for_edges=400)
# 图 3：测地线距离矩阵的热力图（可视化距离结构）
plt.figure(figsize=(6,5))
plt.imshow(dist_matrix, cmap='plasma', aspect='auto')
plt.colorbar(label='geodesic distance')
plt.title("测地线距离矩阵（heatmap）")
plt.tight_layout()
plt.show()
# 图 4：Isomap 结果（2D 嵌入），按 t 着色
plt.figure(figsize=(7,5))
plt.scatter(Y_2d[:,0], Y_2d[:,1], c=t, cmap='rainbow', s=8)
plt.title("Isomap 2D 嵌入结果（按 t 着色）")
plt.colorbar(label='t parameter')
plt.xlabel('dim 1')
plt.ylabel('dim 2')
plt.tight_layout()
plt.show()
# 图 5：残差方差 vs 嵌入维度（1..6），以及特征值谱
res_vars = []
dims = list(range(1,7))
for d in dims:
Yd, eigs = classical_mds_from_dist(dist_matrix, n_components=d)
D_emb = squareform(pdist(Yd))
rv = residual_variance(dist_matrix, D_emb)
res_vars.append(rv)
plt.figure(figsize=(7,4))
plt.plot(dims, res_vars, marker='o', color='crimson')
plt.xticks(dims)
plt.xlabel('嵌入维度 d')
plt.ylabel('Residual Variance (1 - R^2)')
plt.title('Residual Variance vs Embedding Dimension')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
# 额外绘制特征值谱（scree plot）
plt.figure(figsize=(7,4))
eigvals_sorted = eigvals[:50] # 前 50 个
plt.plot(np.arange(1, len(eigvals_sorted)+1), eigvals_sorted, marker='o', color='tab:blue')
plt.xlabel('特征值索引 (降序)')
plt.ylabel('特征值')
plt.title('Classical MDS 特征值谱（前若干）')
plt.grid(alpha=0.3)
plt.tight_layout()
plt.show()
核心步骤
生成数据：
-
我们用参数 以及高度 生成瑞士卷： 。这里 可以看作流形上的“坐标”，用它来给点着色，便于观察能否被 Isomap 展开成一维（或二维）结构。
构造 k-NN 邻接图：
-
先用 PyTorch 的 torch.cdist计算所有点间的欧氏距离矩阵 ； -
对每个点找到最近的 个邻居，将这些边加入邻接矩阵 ，边权为欧氏距离。对称化确保无向图（取最小值或平均都可以，这里取最小）。
计算测地线距离：
-
在邻接图上求最短路径（Dijkstra），得到一个完整的点对测地线距离矩阵 。这一步是核心：将“局部欧氏距离”通过图的连通性拼成“全局测地线距离”。
经典 MDS：
-
将距离矩阵平方化： ； -
中心化矩阵 ； -
构造 Gram 矩阵 ； -
对 做特征分解，取前 个正特征值与对应特征向量，嵌入为 。 -
这是经典 MDS 的标准公式，能够把“距离信息”转为“点坐标”。
评估指标：
-
用来评估低维嵌入保持原始测地线距离的程度。 -
定义为 ，其中 是原始测地线距离向量与嵌入后欧氏距离向量之间的 Pearson 相关系数。值越小，说明嵌入越好保留距离关系。
可视化分析
Swiss Roll 原始 3D:
观察原始数据的形状与颜色连续性。若 Isomap 有效，后续 2D 图中颜色会呈现从一端到另一端的平滑变化（表示 t 顺序被恢复）。
k-NN 邻接图:
直观检查邻接图是否合理。若 k 过小，图可能断开（孤立的小块）；若 k 过大，图可能出现将不同“层”连接的长边（捷径），破坏测地线估计。
测地线距离矩阵热力图:
检查点与点之间测地线距离的全局结构。若流形结构被很好的保留，距离矩阵的近邻会形成对角附近的“暖色带”。
Isomap 2D 嵌入:
评估 Isomap 是否把瑞士卷“展开”。如果成功，颜色应该沿低维空间形成一条连续的带或路径。
如果颜色混乱或折叠，说明参数 k 选择不当或数据噪声影响。
残差方差 & 特征值谱:
残差方差随维度的下降说明低维表示在保留原始测地线信息方面的改进。通常我们在残差下降趋于平缓的拐点处选择维度。
特征值谱展示了前若干维度的重要性。特征值快速衰减意味着一个小维度能解释大部分结构。
总结
Isomap 的魅力在于简单直观：先把局部欧氏信息连成图，再通过图最短路得到近似的测地线距离，最后用经典 MDS 把距离“转化”为低维坐标”；
但 Isomap 对参数 k 和数据噪声敏感，且对大规模数据计算开销大，需要结合近似技术（比如 Landmark）或预处理来改进性能；
最后
