ISOMAP 降维算法的 MATLAB 实现

1. 核心函数 isomap_embedding.m

function [Y, D_geo] = isomap_embedding(X, k, d)
% ISOMAP_EMBEDDING  等距特征映射降维
%   [Y, D_geo] = isomap_embedding(X, k, d)
%
% 输入:
%   X : N×D 矩阵，每行一个样本
%   k : 近邻个数
%   d : 目标维度 (d <= D)
%
% 输出:
%   Y     : N×d 降维后的坐标
%   D_geo : N×N 测地距离矩阵
%
% 依赖: Statistics Toolbox (pdist, knnsearch) 
%       Graph Toolbox (graph, distances)

    N = size(X, 1);
    
    % ---------- Step 1: 计算欧氏距离矩阵 ----------
    fprintf('Step 1: Computing Euclidean distances...\n');
    D_euc = squareform(pdist(X, 'euclidean'));   % N×N
    
    % ---------- Step 2: 构建 k-NN 邻域图 ----------
    fprintf('Step 2: Building %d-NN graph...\n', k);
    [idx, ~] = knnsearch(X, X, 'K', k+1);        % 最近邻包含自身
    idx = idx(:, 2:end);                         % 去掉自身
    
    % 建立稀疏邻接矩阵 (边权 = 欧氏距离)
    A = sparse(N, N);
    for i = 1:N
        neighbors = idx(i, :);
        for j = neighbors
            if A(i,j) == 0 && i ~= j
                A(i,j) = D_euc(i,j);
                A(j,i) = D_euc(i,j);
            end
        end
    end
    
    % 检查图是否连通
    G = graph(A);
    bins = conncomp(G);
    if max(bins) > 1
        warning('图不连通！共有 %d 个连通分量。请增大 k 或使用 Landmark ISOMAP。', max(bins));
    end
    
    % ---------- Step 3: 计算测地距离 (最短路径) ----------
    fprintf('Step 3: Computing geodesic distances via shortest paths...\n');
    D_geo = distances(G);                        % 返回 N×N 矩阵，不可达为 Inf
    
    % 如果有无穷大，将其设为一个大数（但最好保证连通）
    D_geo(isinf(D_geo)) = 10 * max(D_geo(isfinite(D_geo)));
    
    % ---------- Step 4: 经典 MDS ----------
    fprintf('Step 4: Performing MDS to %d dimensions...\n', d);
    Y = cmdscale(D_geo, d);                      % cmdscale 直接完成双中心化和特征分解
    
    % 注：cmdscale 返回的 Y 已经是 d 维坐标，且默认按特征值降序排列
end

---

2. 演示脚本：瑞士卷数据降维

%% ISOMAP 演示 — 瑞士卷展开
clear; clc; close all;

% ---------- 生成瑞士卷数据 ----------
N = 1500;                     % 样本数
t = (3*pi/2)*(1+2*rand(N,1)); % 角度
height = rand(N,1);           % 高度

X = [t.*cos(t), height, t.*sin(t)];   % 三维瑞士卷
X = X - mean(X);                       % 中心化

% ---------- ISOMAP 降维 ----------
k = 12;              % 近邻数（需保证连通）
d = 2;               % 目标维度

[Y, D_geo] = isomap_embedding(X, k, d);

% ---------- 可视化 ----------
figure('Position', [100 100 1200 500]);

subplot(1,2,1);
scatter3(X(:,1), X(:,2), X(:,3), 12, t, 'filled');
title('原始瑞士卷 (3D)');
colormap jet; axis equal; view([45 30]);
xlabel('x'); ylabel('y'); zlabel('z');

subplot(1,2,2);
scatter(Y(:,1), Y(:,2), 12, t, 'filled');
title(sprintf('ISOMAP 降维至 2D (k=%d)', k));
colormap jet; axis equal;
xlabel('Y_1'); ylabel('Y_2');

sgtitle('ISOMAP 等距特征映射');

运行结果：原始瑞士卷被“展开”成一个扇形平面，颜色根据原始角度渐变，说明流形结构被正确恢复。

---

参考代码 isomap-数据降维算法 www.youwenfan.com/contentcnv/81399.html

3. 关键注意事项

3.1 图连通性

• 若提示图不连通，请逐步增大 k 直至连通（或使用 Landmark ISOMAP 处理大规模数据）。
• 可通过 conncomp(graph(A)) 检查分量数量。

3.2 新样本嵌入（Out-of-Sample）

ISOMAP 本身不支持新数据点直接嵌入。常用近似方法：

• Landmark ISOMAP：只对少量锚点（landmarks）做全 ISOMAP，其余点通过三角测量嵌入。
• 训练一个回归模型：用降维后的坐标作为标签，训练从高维到低维的映射（如神经网络、核岭回归）。

3.3 大规模数据

• 当 \(N > 3000\) 时，全距离矩阵 \(O(N^2)\) 和最短路径 \(O(N^2 \log N)\) 可能很慢。
• 改用 Landmark ISOMAP（cmdscale 也支持只传入部分距离矩阵）。

3.4 参数选择

• k 的推荐范围：\(2 \times d_{intrinsic} \sim N/10\)，具体需尝试。
• 目标维度 d 可通过观察 cmdscale 输出的特征值衰减曲线决定（[Y,eigvals] = cmdscale(D_geo); plot(eigvals,'o')）。