简介：正交匹配追踪（OMP）是压缩感知与稀疏表示领域的基础算法，以贪心策略从原子库中逐步选取最相关成分。这份资源以MATLAB实现为核心，从信号模型、残差计算到原子选择与最小二乘更新，系统讲解OMP的完整迭代流程和终止条件，面向信号处理、图像重构方向的初学者及需要快速上手稀疏编码的工程师。压缩包共6个文件，包括4个M脚本与2个DOCX文档：前者覆盖OMP主函数、基于迭代次数和稀疏度的重构测试脚本，后者提供算法原理推导、参数设置及常见问题说明；整个压缩包仅293KB，轻量易用，并配有清晰目录结构。目前已有1396人学习浏览，内容实用性得到验证。借助这些代码与文档，读者能够理解支持集扩充、残差更新等关键步骤，直接运行测试脚本观察不同采样率下的重构误差变化，并通过调整迭代次数与容差阈值适配不同信号场景，为后续研究压缩感知、LASSO等进阶方法打下扎实基础。

1. 先说清楚 OMP 算法到底解决什么问题

OMP（Orthogonal Matching Pursuit，正交匹配追踪）是我在压缩感知和稀疏表示里用得最多的算法，没有之一。它的典型场景是：已知一个过完备字典 A 和观测向量 y，想从 y = Ax 里恢复尽可能稀疏的系数 x，同时 x 中非零项越少越好。这个问题在图像压缩感知、雷达成像、阵列信号处理、甚至深度模型剪枝里都会碰到。反直觉的地方在于：方程 y = Ax 本身欠定，解有无穷多个，但 OMP 不去解 L1 范数优化，而是靠“每次挑一个最相关的原子、做一次正交投影、更新一次残差”的贪婪迭代，就能在许多情况下把真正的稀疏解找回来。因为它逻辑简单、实现门槛低，MATLAB 里手写一个核心循环通常不到 30 行，而且调试比大多数凸优化工具更直观。这篇我会把原理、参数、坑一次讲透，并给出一套可以直接跑通的代码。

2. OMP 的原理与 MATLAB 最小实现：从残差内积到正交投影

2.1 把 OMP 的数学过程压缩成“选原子、投影、更新残差”三步

OMP 的输入是测量矩阵或字典 ( A \in \mathbb{R}^{N \times M} )（N 远小于 M 时就是欠定）、观测向量 ( y \in \mathbb{R}^{N \times 1} )、稀疏度 K（或残差容差 tol）。它要解决的问题是：

[ \min_x |x|_0 \quad \text{s.t.} \quad y = Ax ]

这里的 ( |x|_0 ) 是 x 中非零元素的个数。直接穷举所有可能的非零位置是 NP 难，所以 OMP 改用贪婪策略迭代逼近。整个算法的骨架只有三步：

第一步，计算当前残差 r 与字典每一列的内积 ( |\langle r, a_j\rangle| )，取绝对值最大的原子下标 j，加入支撑集 S； 第二步，在支撑集 S 上做最小二乘投影：( x_S = (A_S^T A_S)^{-1} A_S^T y )，这里 A_S 是由已选原子组成的子矩阵； 第三步，用这个投影结果更新残差 ( r = y - A_S x_S )，回到第一步。

“正交”两个字落在第二步：每次求系数不是简单的相关值，而是在选中的原子张成的子空间里做一次正交投影，因此已经被选过的原子不会因为后续选了新原子而过时，迭代过程从原理上是保证残差单调下降的。这也是它比 Matching Pursuit 更稳的核心原因。

Mind the 一步：整个算法能成立的前提是字典的每一列都经过归一化处理。如果某些原子能量特别大，内积绝对值天然就高，会误导选择。最稳妥的顺序是提前把 A 的每一列除以它的 L2 范数，恢复出的稀疏系数再在后续恢复阶段乘回对应范数。下面的实现默认你已经这么做了。

2.2 手写 OMP 的 MATLAB 最小循环版：30 行跑通稀疏恢复

直接给一个最经典的实现，它强调“看得懂逻辑”，而不是“性能最优”：

function x = omp_solve(A, y, K, tol) % A : N×M 字典矩阵，调用前请先做列归一化 % y : N×1 观测向量 % K : 最大稀疏度（迭代上限） % tol : 残差能量阈值，不限制可设 Inf [N, M] = size(A); r = y; % 初始残差 = 观测 selected = false(M, 1); % 记录哪些原子被选过 support = []; % 支撑集下标 x = zeros(M, 1); for iter = 1:K % 1) 计算残差与所有原子的相关性 corr = A' * r; [~, j] = max(abs(corr)); % 防止重复选择同一个原子 if selected(j) corr(j) = -inf; [~, j] = max(abs(corr)); end selected(j) = true; support = [support, j]; %#ok<AGROW> % 2) 在支撑集上做最小二乘（正交投影） As = A(:, support); xs = As \ y; % 等价于 pinv(As)*y，但数值上更稳 r = y - As * xs; % 3) 残差能量达标则提前停止 if norm(r) <= tol break; end end x(support) = xs; end

这段代码对应上面三步，核心只有三行：corr = A' * r做相关检测，As \ y做正交投影，r = y - As * xs更新残差。几个参数要说明：K给的是最大迭代次数，即人为假设的稀疏度上限；tol给的是残差能量阈值，适合噪声已知的场合。二者是“或”的关系，先触发哪个就停在哪个。典型做法是头几次运行只给K不给tol，用来观察残差下降曲线；确认拐点后再把tol加进去。

这个版本的问题也很明显：每次循环都要重新计算A' * r（一个 N×M 的矩阵乘法），还要对As重新做一次最小二乘，循环 K 次就是 K 次完整的大矩阵运算。当 N 和 M 都到几千、K 到几百时，运行时间会很难看。所以实际工程里我不会用这个版本上生产，它可以作为功能验证和调试的基础。

2.3 用预计算 Gram 矩阵做加速：避免重复做最小二乘

上面最小版本里最耗费时间的部分，是每轮循环都重新对 A_S 求解一次最小二乘。如果我们提前把 Gram 矩阵算好，就能在循环中只操作支撑集对应的子矩阵。原理是：

[ x_S = (A_S^T A_S)^{-1} A_S^T y = (A_S^T A_S)^{-1} (A^T y)_S ]

其中 ( A^T y ) 是一次性算好的，支撑集变化时只需要从里面取对应行。于是每次迭代只需做一个小的线性方程组求解，而不是重新构造大矩阵。实现如下：

function x = omp_gram(A, y, K, tol) % 预计算版本：适合原子数多、每轮支撑集较小时的加速 G = A' * A; % Gram 矩阵，M×M Aty = A' * y; % 预计算观测与所有原子的内积 r = y; nAtoms = size(A, 2); selected = false(nAtoms, 1); support = []; x = zeros(nAtoms, 1); for iter = 1:K corr = A' * r; [~, j] = max(abs(corr)); if selected(j) corr(j) = -inf; [~, j] = max(abs(corr)); end selected(j) = true; support = [support, j]; %#ok<AGROW> % 从 Gram 子矩阵求解，避免重复构造 As Gs = G(support, support); xs = Gs \ Aty(support); % 等价于 A(:,support) \ y r = y - A(:, support) * xs; if norm(r) <= tol break; end end x(support) = xs; end

这个版本把 OMP 循环核心从“大矩阵乘 + 最小二乘”降成了“小 Gram 子矩阵求逆 + 稀疏矩阵乘”。当 M 在几千、K 在几十这个量级时，提速通常两三倍以上。但需要留意，Gs \ Aty(support)的数值稳定性和直接用A(:, support) \ y不完全一样。如果 Gram 矩阵有条件数很大的列子集，后者的 QR 分解会自动做主元缩放，前者的普通\也有可能内部做同样处理，但总体上前者更容易受到原子共线性的影响。所以如果字典相干性比较高，我建议回到直接最小二乘版本，或者用下一章讲到的 QR 更新。

另外再补充一个很实用的变体：当字典是正交基或近似正交基时，A(:, support) \ y可以直接换成A(:, support)' * y，此时正交投影退化为内积，OMP 就成了纯贪婪选择。这种写法在做正交基字典（如 DCT、小波基）时非常快，但千万不要把它用在过完备字典上，否则就是 MP 算法而不是 OMP，精度会明显下降。

3. MATLAB 里 OMP 的四个必调参数与字典设计：稀疏度、容差和相干性

3.1 稀疏度 K 与残差容差 tol：两个终止条件怎么配合

OMP 的终止条件有两个：最大迭代次数 K 和残差容差 tol。新手最常见的困惑是到底用哪个。我的经验是：如果对信号的稀疏程度有先验，K 优先；如果已知噪声方差，tol 优先；现实中两者一起用，谁先满足谁停。

稀疏度 K 的估计没有万能公式，但有两条常用路径。第一，如果信号来自某个已知的变换域（比如图像 DCT 系数），统计训练集里保留 95% 能量所需的系数个数，取个上界当 K。第二，如果完全没有先验，就跑一遍不带 tol 的 OMP，记录每一轮的残差能量变化。你会发现残差通常先快速下降，后进入平缓区，平缓区的起点就是合理的 K。超过这个点再迭代，拟合法就会把噪声也拟合进去，恢复结果反而变差。

残差容差 tol 的设定和噪声水平直接挂钩。对于 y = Ax + n，若 n 是方差 ( \sigma_n^2 ) 的高斯白噪声，残差 ( |r|_2^2 ) 的期望约等于 ( N \sigma_n^2 )。我一般把 tol 设在 ( \sqrt{N} \sigma_n ) 的 1.05 到 1.2 倍之间，太低会让算法多迭代几次去拟合噪声，太高又会提前停止导致欠拟合。MATLAB 里如果知道信噪比，可以直接转换：

sigma_n = std(y) / (10^(SNR_dB / 20)); tol = 1.1 * sqrt(length(y)) * sigma_n;

注意这只是一个经验起点，具体还要看字典的相干性。相干性越高，同样的 tol 下恢复误差会越大，这个下面展开。

3.2 字典 A 的归一化与非相干性：不归一化一定会翻车

字典设计是 OMP 里最容易翻车的环节。第一个坑是没做列归一化。比如你有两个等能量原子 a1 和 a2，但 a1 的 L2 范数是 2，a2 是 0.5，那么即使残差 r 与 a2 的方向完全一致，内积绝对值也可能比 a1 小。OMP 依据的就是内积绝对值最大，最终会错误地选成 a1。所以无论什么场景，我拿到字典的第一件事就是A = A ./ vecnorm(A)。

第二个坑是原子相干性过高。两个原子方向接近，残差投影在它们上的内积会同时很高，OMP 可能先选了其中一个，下一轮又把另一个当“新信息”选进来，造成支撑集里出现一对近似重复的原子。描述这个现象的量是互相干系数：

[ \mu = \max_{i \neq j} |\langle a_i, a_j\rangle| ]

MATLAB 里检查非常方便：

% 假设 A 已经按列归一化 G = A' * A; G(1:size(A, 2) + 1:end) = 0; % 把对角线置零 mu = max(abs(G(:))); fprintf('mutual coherence = %.3f\n', mu);

当 mu 大于 0.8 时，OMP 的支撑集可靠性会明显下降；大于 0.95 时基本可以预测恢复会乱选。如果 mu 偏高，常见处理是把字典改造成正交化字典（比如用 QR 对字典做一次预处理，但这个会改变原子的物理含义，只适合对表示系数做中间计算），或者改用更鲁棒的原子选择策略——比如同时看内积和残差投影增益，而不是只看内积绝对值。

3.3 用 OMP 做图像分块重建：块尺寸与字典的张成关系

图像压缩感知是最常见的 OMP 落地场景之一。做法是先把图像切成一堆小块，比如 8×8 或 16×16 像素，每块拉直成一个 64 或 256 维向量，分别做压缩采样和稀疏恢复，最后拼回整幅图。块尺寸的选择直接影响字典设计：块越小，字典原子维度越低，计算越快，但块与块之间容易产生明显的接缝；块越大，恢复质量越好，但测量矩阵和字典规模都水涨船高。

以 8×8 块为例，一种常见配置是：块拉直成 64 维，用 DCT 基作为稀疏字典 Psi，维度为 64×64；构造一个 ( m \times 64 ) 的高斯随机测量矩阵 Phi，其中 m 是每块采样数，比如 32；把 Phi 和 Psi 相乘得到传感矩阵 A = Phi * Psi，然后用 OMP 恢复每个块的 DCT 系数 theta，再反变换回像素域。核心代码是：

blk = 8; m = 32; % 每块观测数，压缩比 = m / 64 Phi = randn(m, blk^2) / sqrt(m); % 高斯随机投影，行数远小于列数 Psi = dctmtx(blk^2); % 离散余弦变换基 A = Phi * Psi; % 等效传感矩阵 A = A ./ vecnorm(A); % 列归一化（重要） % 对每个图像块做恢复 for bi = 1:blk:size(X, 1) for bj = 1:blk:size(X, 2) patch = X(bi:bi+blk-1, bj:bj+blk-1); y = Phi * patch(:); % 压缩观测 theta = omp_solve(A, y, K, tol); % OMP 恢复 DCT 系数 recons_patch = Psi * theta; % 逆变换回像素域 Xrec(bi:bi+blk-1, bj:bj+blk-1) = reshape(recons_patch, blk, blk); end end

这里面我把 A 的归一化放在了测量矩阵和稀疏字典相乘之后，这是 OMP 的要求。但要注意：归一化后的 A 对应的稀疏系数 theta 已经是“被缩放过”的表示，最后重建时直接乘 Psi 会让结果偏小或偏大。正确做法是记住每列归一化因子，或者让omp_solve返回归一化前的系数。我通常在omp_solve内部只做列归一化，而在返回前用norm_col把系数映射回去：

x = x ./ norm_col'; % norm_col 是归一化时除掉的每个原子范数

3.4 测量矩阵选择：高斯随机阵不是唯一答案

除了字典 A，测量矩阵 Phi 的选择也会对 OMP 恢复质量产生影响。常见的做法是随机高斯矩阵，行数 m 至少是稀疏度 K 的 2~3 倍，条件数接近 1，理论上能和大部分字典配合。但实际工程中你可能会遇到几种情况：观测矩阵是哈达玛矩阵（硬件实现方便但需要满足二进约束），或者是部分傅里叶矩阵。部分傅里叶矩阵和 DCT 字典合用时，传感矩阵 A 的行向量会高度周期化，OMP 的残差更新容易出现振荡。

我的一般策略是：先用高斯矩阵跑通流程，然后换成目标测量矩阵，对比稀疏度-残差曲线。如果换成目标矩阵后同一 K 值下重建误差变大超过 10%，通常不是 OMP 的问题，而是测量矩阵与字典相干性过高，这个信号本身就是难恢复的。

4. OMP 在 MATLAB 里的 5 个高频坑：现象、原因与解决

4.1 原子未归一化，残差降不下去

现象：迭代次数增加了，但残差能量迟迟不降，恢复出的信号到处是尖刺。原因：字典各列范数差异大，OMP 的相关度检测被大范数列带偏。解决：所有 OMP 实现前先做A = A ./ vecnorm(A)；如果字典来自训练数据，归一化因子要存下来，恢复后映射回真实系数。

4.2 稀疏度 K 设置过大，把噪声当信号拟合

现象：K 增大到一定程度后，信噪比反而下降，甚至在无噪声的验证数据上出现重建伪影。原因：真实稀疏度是 8，你设成 20，OMP 在后面十几轮里就是拿残差里的噪声分量做最小二乘。解决：观察残差下降曲线的“膝盖”，在膝盖处截断；或者把 tol 也加上，让残差能量先达到阈值成为停止条件。

4.3 重复选择同一个原子导致死循环

现象：程序跑很久不结束，或者support里出现两个相同下标。原因：残差在某个方向上仍有较大分量，而该方向已被选过，但由于字典原子之间非正交，投影后的残差又和该原子保留了一定相关性，于是每轮都挑中它。解决：代码里必须有selected标记位，任何情况下不能重复选同一下标；如果连续两轮选出的原子相关度都大于 0.9，建议检查字典质量，而不是继续调 OMP 参数。

4.4 Gram 子矩阵病态，最小二乘结果数值爆炸

现象：支撑集很小的时候结果正常，支撑集稍大就开始出现 1e5 量级的系数，恢复信号完全失真。原因：支撑集里有两个高度相关的原子，导致 ( A_S^T A_S ) 接近奇异，\求解被条件数放大。解决：优先改用A(:, support) \ y，让 MATLAB 内部走 QR 分解；避免用inv(A_S' * A_S)这种裸写法。另一个有效手段是给 Gram 子矩阵加一个微小的岭正则：(Gs + 1e-6 * eye(k)) \ Aty(support)，牺牲一点残差换稳定性，这在高噪声条件下效果明显。

4.5 把 OMP 当黑匣子，参数和字典不匹配

现象：同样的 K 和 tol，换了一组训练数据后恢复效果大幅波动。原因：固定测量矩阵但更换了信号源，实际稀疏度变化了，原来的 K 不再合适。解决：做一个两步自适应流程：第一次运行用较大的 K 跑完，记录每轮残差；第二次用拐点位置作为 K 重新运行。这段“两遍 OMP”的开销不大，但能明显提升稳定性。

提示：上述所有坑都建议先用合成数据复现一遍，再看真实数据。因为真实数据的“真实稀疏度”你也不知道，出了问题很难判断是算法错还是参数错。

5. 用合成信号验证 OMP 的正确性，再换 QR 更新提速

最后给你一个我反复使用的验证脚本套路，以及一个工程提速技巧。验证永远从合成数据开始，因为只有合成数据你知道真实的稀疏支撑集在哪里。构造方法很直接：随机生成一个稀疏向量 x_true，非零位置随机，非零值用标准的randn；再随机生成高斯字典 A；观测 y = A * x_true，这时无噪声，跑 OMP 后对比支撑集是否正确恢复。然后依次叠加噪声，观察不同信噪比下的恢复误差，画出 SNR 和支撑集恢复率曲线。这一步建议写成函数固化下来，之后每次改字典、改参数都跑一遍，能省大量排错时间。

提速技巧方面，如果字典维度很大，上述的 Gram 版本还不够快，可以考虑在支撑集更新时维护一个 QR 分解。具体做法是：每轮新增一列 a_j，把它和上一轮的残差正交化后追加到 Q 矩阵尾部，系数求解变成xs = R \ (Q' * y)。这样每轮只需更新一次 Q 和 R，无需对 A_S 重新做 QR。MATLAB 里大致写法是：

[Q, R] = qr(A(:, support), 0); % 初始支撑集 % 新增原子 a_new 时 q_new = a_new - Q * (Q' * a_new); q_new = q_new / norm(q_new); Q = [Q, q_new]; R = [R, Q_old' * a_new; zeros(1, size(R,1)), norm(q_new)]; xs = R \ (Q' * y);

代码中需要注意 Gram-Schmidt 的重正交化，否则两次迭代后数值误差会累积。这种写法的提速在小支撑集上并不明显，但当 K 在几百、字典 M 在几万时，能明显缩短每次迭代的线性代数耗时。

我的习惯是：第一版永远是最小循环版，验证正确性；确认逻辑没问题后，才换成 Gram 版或 QR 更新版。直接一上来写加速版，出了问题你会分不清是算法逻辑错还是数值优化引入的 bug。OMP 本身不难，但省去验证、跳过归一化、没有残差监控这三件事，我至少各翻车过两次。希望这篇能帮你少走这些弯路。

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