企业官网建设流程全解析

简介：面向高维数据降维与可视化需求，等距映射算法是基于流形学习的经典非线性降维方法，能有效挖掘嵌入在复杂数据中的低维流形结构，弥补主成分分析等线性方法在非线性场景下的不足。这份资源提供了该算法的MATLAB实现，适合机器学习、数据挖掘方向的初学者和科研人员使用，可直接用于课程实验或项目验证。压缩包整体仅1KB，共包含3个m源文件，即主算法程序与距离矩阵计算辅助函数，代码紧凑、结构清晰，便于阅读与二次开发。目前该资源已有402人学习浏览，体现了不错的实用性。通过精读这些代码，读者可以完整掌握该算法的实现流程：依据近邻关系构建邻接图，基于最短路径计算测地距离矩阵，再应用多维缩放得到低维嵌入坐标；同时也能深入理解近邻数等关键参数对降维效果的影响，方便后续替换数据、调整算法或扩展到其他流形学习技术。

1. 为什么流形学习选 ISOMAP：从线性降维到等距映射

PCA、MDS 这类线性降维方法，拿到一张瑞士卷形状的数据时，会不分青红皂白地把两个卷层压在一起，投影结果几乎看不出原本的结构。问题出在欧氏距离上：三维空间里隔得很近的两个点，沿着卷面走可能要绕过半圈，这种沿流形量得的距离才是判断样本关系的正确标尺。ISOMAP（Isometric Mapping，等距映射）正是为这类任务设计的流形学习方法。它先构造近邻图，用图上的最短路径逼近测地距离，再用经典 MDS 将这些距离映射到低维空间，从而在降维的同时保持全局几何结构。对形状分析、姿态估计和高维数据可视化来说，ISOMAP 是一个容易在 MATLAB 里落地、结果也足够直观的起点。

2. ISOMAP 的算法拆解：近邻图、测地距离与 MDS 谱分解

2.1 等距映射的核心思想：用图距离逼近流形测地距离

ISOMAP 的思路可以分三步理解。第一步，在原始高维空间里找出每个样本的近邻，形成一张带权无向图；第二步，计算图上任意两点之间的最短路径，作为测地距离的近似；第三步，把最短路径矩阵交给经典 MDS，得到一组低维坐标，使得低维空间里的欧氏距离尽量等于测地距离。整个过程本质上是用图来近似流形，再用谱分解来保距。

与局部线性嵌入 LLE 只保持局部邻域权重不同，ISOMAP 着重保持所有样本对之间的全局距离，因此它对展开型流形效果显著，而在拓扑结构复杂或有大空洞的数据上容易失真。为什么要绕这么大一圈，而不是直接对高维距离做 MDS？因为高维空间里的直线距离会横穿流形，而测地距离才是流形上的真实位移。比如球面上相距半个大圆的两点，三维欧氏距离是直径，但沿球面的最短路径是大圆弧，二者含义完全不同。ISOMAP 用近邻图上的最短路径去近似这个大圆弧，再进入 MDS，等价于把流形“铺平”后再保持点间距离。

2.2 构建邻域图的两种策略：K 近邻与 ε 邻域

构建近邻图是 ISOMAP 成败的第一步。常见做法有两种：固定近邻数 K，选择每个样本最近的 K 个点连边；或者固定半径 ε，两点距离小于 ε 就连边。两者的取舍可以放进一张表里。

在 MATLAB 里构建 K 近邻图，最简单的是用 pdist2 算出距离矩阵，再对每一行排序取前 K+1 个，第一个是自身。也可以用 knnsearch 直接返回邻居索引。我更习惯先保留完整距离矩阵，因为后面 MDS 和残差计算还要用到。下面这段代码演示了带权邻接矩阵的构建，边权取两点间的欧氏距离，随后用 max 做对称化，避免两个点之间出现方向不同的两条边。

D = pdist2(X, X, 'euclidean'); n = size(X, 1); A = sparse(n, n); for i = 1:n [~, idx] = sort(D(i, :)); nbrs = idx(2:K+1); % 去掉自身 A(i, nbrs) = D(i, nbrs); end A = max(A, A'); % 对称化

这段代码中，X 是 N×D 的样本矩阵，K 是近邻数。对每一行排序后取第 2 到第 K+1 个索引，因为第 1 个是自身，距离为 0。用 max 做对称化而不是加法平均，是为了防止 A(i,j) 和 A(j,i) 两条边同时存在时最短路径算法把权重算成两倍。如果改用 knnsearch，构建会更快，但需要额外保存邻居距离，整体差别不大。

2.3 测地距离矩阵：用 graph 对象做最短路径

图上最短路径的经典算法是 Floyd-Warshall，直接在稠密距离矩阵上迭代，但复杂度是 O(n^3)。在 MATLAB 中，我一般会改用 graph 对象上的 distances 函数。graph 内部根据稀疏矩阵的密度选择 Dijkstra 或 Johnson 算法，在近邻图这种边数远小于 n^2 的稀疏图上，速度往往比 Floyd 快一到两个数量级。虽然 Floyd 的写法很优雅，但实际数据上跑到一万个样本就会卡到无法忍受。

G = graph(A); % 从稀疏邻接矩阵创建无向图 Dgeo = distances(G); % 返回所有点对最短路径，Inf 表示不连通 Dgeo(isinf(Dgeo)) = max(Dgeo(~isinf(Dgeo))) * 2;

这里把 Inf 替换成一个比最大有效距离大很多的值，是工程上的应急处理，目的是不让后面 MDS 的双中心化产生 NaN。如果断掉的连通分量较大，这种替换本身会引入严重偏差；更严格的做法是先提取最大连通分量，再在分量上运行 ISOMAP。关于连通分量的处理，第 3 章会给出更严谨的替代方案。

2.4 经典 MDS 降维与特征值分解

拿到测地距离矩阵 Dgeo 后，ISOMAP 的收尾工作是经典 MDS。先对平方距离矩阵做双中心化，构造 Gram 矩阵 B = -1/2 * J * D^2 * J，其中 J = I - 1/n * 11^T 是中心化矩阵。然后对 B 做特征值分解，取前 d 个最大特征值对应的特征向量，低维坐标就是特征向量乘以对应特征值的平方根。这样得到的低维空间中，点间欧氏距离在最小二乘意义下最优地逼近测地距离。

J = eye(n) - ones(n) / n; B = -0.5 * J * (Dgeo .^ 2) * J; [V, E] = eig(B); ev = real(diag(E)); [~, idx] = sort(ev, 'descend'); Y = V(:, idx(1:d)) * diag(sqrt(max(ev(idx(1:d)), 0)));

这里使用 eig 而不是 eigs，是因为 B 通常都是满秩矩阵，eig 能直接拿到全部特征向量；当样本数上万时再考虑用 eigs 只求前几个。sqrt(max(...)) 是为了避免负特征值产生复数坐标。测地距离矩阵在噪声下不保证是欧氏距离矩阵，所以 B 中会出现负特征值，这并不影响前 d 个主要特征向量的有效性，但需要截断。

3. MATLAB 实现 ISOMAP：从零写一个可复用的 isomap 函数

3.1 输入输出设计与参数校验

直接写一个函数，而不是只贴散装代码。接口设计为 [Y, R] = isomap(X, k, d)。X 是 N×D 矩阵，每个样本一行；k 是近邻数；d 是目标维度。输出 Y 是 N×d 的低维坐标，R 是保距残差，用来判断降维质量。参数校验放在函数开头，避免后面用到 k 或 d 时产生掩码式错误。

n = size(X, 1); if k < 2 || k >= n error('k 必须在 [2, n-1] 之间'); end if d < 1 || d >= min(n, size(X, 2)) error('d 超出合法范围'); end

这里把 k 的下界定在 2，是因为 k=1 时近邻图只是一条条孤立边，测地距离与欧氏距离几乎没有区别，流形学习失去意义。d 的上界受限于样本数和原始维度，毕竟低维坐标最多只能有 min(n, D)-1 个非零特征值。如果你的 MATLAB 版本比较新，还可以把这段校验放到 arguments 代码块里，但改写成函数后处理报错信息更直观。

3.2 完整函数主体：近邻图、最短路径与经典 MDS

把第 2 章的散装步骤合成一个 isomap.m。完整代码不长，核心就是 pdist2、graph/distances 和 eig 三句话。这里特意保留 for 循环构建近邻图，是希望你能在断点处观察邻居索引；如果追求性能，可以把内层替换成 knnsearch。

function [Y, R] = isomap(X, k, d) n = size(X, 1); D = pdist2(X, X, 'euclidean'); D(1:n+1:end) = Inf; % 排除自身 A = zeros(n, n); for i = 1:n [~, ord] = sort(D(i, :)); nbrs = ord(1:k); A(i, nbrs) = D(i, nbrs); end A = max(A, A'); G = graph(A); Dgeo = distances(G); if any(isinf(Dgeo(:))) Dgeo(isinf(Dgeo)) = max(Dgeo(isfinite(Dgeo))) * 2; end J = eye(n) - ones(n) / n; B = -0.5 * J * (Dgeo .^ 2) * J; [V, E] = eig(B); ev = real(diag(E)); [~, idx] = sort(ev, 'descend'); Y = V(:, idx(1:d)) * diag(sqrt(max(ev(idx(1:d)), 0))); Ydist = pdist(Y, 'euclidean'); dvec = Dgeo(tril(true(n), -1)); R = 1 - corr(dvec(:), Ydist(:)); if R < 0, R = 0; end end

参数说明：X 必须是数值型矩阵，缺失值需要提前处理，不能带入 pdist2。k 的选择直接影响 A 的边数，一般从 min(10, n-1) 附近开始扫描。d 是目标维度，通常先设为 2 或 3 做可视化，再根据残差曲线调整。corr 来自 Statistics Toolbox，如果没有这个工具箱，可以自己算皮尔逊相关系数，公式是 (x-mean(x))'*(y-mean(y)) 除以标准差乘积。用 tril(true(n),-1) 提取 Dgeo 下三角，是为了和 pdist 输出的向量顺序对齐，避免把矩阵上三角重复算进去。

3.3 边界处理：不连通图的替换策略与连通分量检查

3.2 的代码用 max(finiteVals)*2 替换 Inf，能在断图时保住输出不为 NaN，但这属于“尽力而为”。更严谨的做法是先检查连通分量，如果最大连通分量只覆盖了大部分样本，就只在该分量上降维，并把孤立样本的坐标置为 0 或 NaN。下面这段代码可以在进入 MDS 之前使用。

G = graph(A); bins = conncomp(G); counts = accumarray(bins(:), 1); [~, maxBin] = max(counts); mainIdx = find(bins == maxBin);

conncomp 返回每个节点所属分量的编号。max(counts) 定位包含节点最多的分量，mainIdx 就是该分量内的样本索引。后续只需要把 isomap 的输入 X 替换成 X(mainIdx, :)，算完后重新映射到原图位置。如果多个分量体量接近，说明数据本身可以被切成多块独立流形，强行用一个低维坐标表示会失真，这时可以考虑分簇后分别降维。

3.4 特征值分解与残差计算的实现细节

经典 MDS 的特征值分解有一个容易被忽略的点：B = -0.5 * J * D^2 * J 是数值上对称的，但由于浮点误差，eig 返回的特征值可能有微小的虚部，所以代码里用 real 取实部。排序用 sort(ev, 'descend')，取前 d 个之后还要用 max(...,0) 做截断，因为负特征值开根号会得到复数。若发现 Y 中出现大量全零列，多半是 d 超过了正特征值个数，这时需要减小 d。

残差 R = 1 - corr(测地距离, 低维距离) 是 ISOMAP 经典定义的一个变体。它衡量的是降维前后点对距离的单调相关性，R 越接近 0 表示保距效果越好。需要注意，corr 对尺度缩放不敏感，ISOMAP 本身也只要求相对距离一致，因此这个指标比直接算平均绝对误差更适合判断流形展开质量。

4. 参数选择与调优：K 近邻、特征维度与噪声数据的坑

4.1 近邻数 K 对测地线失真的影响

K 的选择是 ISOMAP 最敏感的参数。K 太小，近邻图可能被拆成多个连通分量，测地距离矩阵里出现大量 Inf；K 太大，边缘处原本不相邻的两个卷层会被一条捷径连起来，测地距离被严重低估，展开结果出现重叠。用 MATLAB 调参时，可以做一个扫描，例如 K 从 5 递增到 20，分别运行 isomap，记录残差 R。

ks = 5:20; res = zeros(size(ks)); for i = 1:numel(ks) [~, rr] = isomap(X, ks(i), 2); res(i) = rr; end plot(ks, res, 'o-'); xlabel('K'); ylabel('Residual');

通常残差随 K 先下降后上升，选择平台区的左端点。如果所有 K 下残差都很大，说明数据本身不是单一光滑流形，或者距离度量不合适。另一个辅助指标是断边比例，可以在 distances(G) 后统计 Inf 个数，若 K 增大到某一值后断边比例突然归零，通常说明图已经连通，依然存在较多断边的 K 值不值得尝试。

4.2 本征维度估计与残差曲线

ISOMAP 的目标维度 d 也不是拍脑袋定的。常见方法是在 d = 1:min(10, size(X,2)-1) 内循环，计算残差，画残差曲线。随着维度增加，残差显著下降后进入平台，拐点处可以看作本征维度。下面的代码沿用自写 isomap 的 R 输出，不需要重新写距离计算。

ds = 1:10; res = zeros(size(ds)); for i = 1:numel(ds) [~, res(i)] = isomap(X, 12, ds(i)); end plot(ds, res, 'o-'); xlabel('d'); ylabel('Residual');

除了残差曲线，还可以观察 B 的特征值衰减。特征值排序后，前几个特征值明显大于其余时，拐点同样指示本征维度。但特征值衰减受样本密度影响大，残差曲线更接近“重构误差”语义。我一般两个图一起看：残差曲线负责选 d，特征图负责交叉验证。

4.3 噪声数据与预处理：为什么 ISOMAP 会“短路”

ISOMAP 对噪声敏感的原因在于近邻图只看欧氏距离，不区分“沿流形”和“横穿流形”。一个噪声点可能把两个不相邻的流形片层连接起来，导致很多点对之间的测地距离被低估。处理噪声的常见做法有三类：先降噪再跑 ISOMAP，比如对局部邻域做 PCA 平滑；增大 K，让单个噪声点的影响被周围点稀释；或者改用鲁棒距离，如对距离矩阵做分位数截断。我一般在数据维度很高时先做一个 PCA 预降维，保住 95% 方差，再去跑 ISOMAP。这一步对 MATLAB 里的图像特征特别重要，能避免噪声主导近邻排序。

4.4 ISOMAP 与 PCA、LLE、t-SNE 的适用场景对比

做降维选型时，很多人问 ISOMAP 和 LLE、t-SNE 有什么区别。直接看这张表：

ISOMAP 的目标是恢复低维坐标，而不是像 t-SNE 那样把聚类结构按社区摊开。如果你的数据有明显的球面或圆柱结构，ISOMAP 比 PCA 和 t-SNE 更接近真实内禀坐标。如果数据采样稀疏或者存在多个分量，LLE 和 t-SNE 往往更稳。注意 ISOMAP 需要存储 N×N 的距离和最短路径矩阵，样本数超过两万时内存会吃紧，常见做法是先抽样跑参数，再对全量数据用 Nyström 近似。

5. 验证降维效果：残差、保距误差与 MATLAB 可视化

5.1 生成 Swiss Roll 并运行自写 isomap

Swiss roll 是 ISOMAP 的标准试金石。在 MATLAB 中可以用下面的代码生成一个充分采样的三维瑞士卷，然后调用第 3 章的 isomap 函数。

n = 1500; t = (3 * pi / 2) * (1 + 2 * rand(n, 1)); h = 30 * rand(n, 1); X = [t .* cos(t), h, t .* sin(t)]; X = (X - min(X)) ./ (max(X) - min(X)); [Y, R] = isomap(X, 12, 2); scatter(Y(:,1), Y(:,2), 8, t, 'filled');

这里的 t 变量本身对应瑞士卷展开后的角度方向，用它给散点图着色，能看出展开后的 Y 是否像一张被裁开的扇形。如果 ISOMAP 工作正常，Y 的横轴应该大致沿着 t 的变化方向，纵轴对应 h 的方向。如果图上出现明显的卷曲或叠层，优先怀疑 K 太大。

5.2 残差曲线与类可分性交叉验证

残差是判别降维是否有效的一个定量指标，但只有残差不够。如果降维后还要用于分类，可以用 KNN 交叉验证的准确率来对比不同 K 下的嵌入。MATLAB 里可以用 fitcknn 和 crossval 快速完成，但注意 fitcknn 属于 Statistics Toolbox。操作思路是：对每个候选 K 运行 isomap 得到 Y，把 [Y, labels] 交给 fitcknn，然后 crossval 得到损失。这样能得到“降维后信息损失了多少”的实操答案，比单纯看残差更贴近业务目标。

5.3 常见错误与调试技巧

最后放三个最容易踩的坑。第一，pdist2 在样本数过万时单是矩阵就有几百 MB，建议用分块欧氏距离或先采样一部分做参数探索。第二，graph 对象要求邻接矩阵对称且无自环，构建后可以用 issymmetric(A) 和 G.numedges 检查边数是否符合预期。第三，特征分解后出现复数坐标，基本是 Dgeo 里有 NaN 或负特征值没有截断，回查 Inf 替换策略。调试时我习惯打印 min(Dgeo(:))、max(Dgeo(:)) 和 sum(isinf(Dgeo(:))) 三个量，能快速定位断图和异常距离。

本文还有配套的精品资源，点击获取