压缩感知重构算法实战：OMP、CoSaMP、SAMP 选型与调优指南

简介：这份资源面向信号处理、图像重建与无线通信方向的学生、研究人员及工程师，提供压缩感知重构算法的MATLAB实现代码，帮助读者在稀疏信号恢复这一核心问题上快速上手并对比不同策略。压缩包共5个文件，全部为.m脚本，整体约7KB，涵盖CoSaMP、SAMP、SP、ROMP等主流重构算法的具体实现，每个脚本对应一种求解思路，便于逐份阅读与运行验证。已有149人学习关注，说明其在同类学习资料中具备一定参考价值。读者可借助这些代码理解贪婪类算法在支撑集选择、残差更新与迭代终止条件上的差异，并在此基础上调整测量矩阵、稀疏度等参数，观察重构精度与收敛速度的变化；也可将其作为课程实验或课题预研的起点，进一步扩展到MRI成像、遥感数据获取等实际场景中，形成从理论到代码的完整认知链条。

1. 压缩感知重构算法合集：从 CoSaMP 到 SAMP，一套代码怎么选、怎么跑

手里有一份cs-code.zip，里面塞了压缩感知（Compressed Sensing, CS）的多种重构算法，OMP、CoSaMP、SAMP、IRLS 大概率都在。问题从来不是「有没有代码」，而是「我该用哪个、参数怎么设、为什么跑出来跟论文对不上」。压缩感知的核心逻辑就一句话：当信号在某个变换域稀疏时，用远少于奈奎斯特采样率的观测值，通过非线性重构把原信号恢复出来。听起来很美，但落到工程里，测量矩阵怎么生成、稀疏度 K 怎么估、噪声方差怎么设，每一个都能让重构结果从「完美恢复」变成「一团噪声」。这篇笔记面向三类人：刚拿到这套代码想跑通第一个 demo 的新手、在 CoSaMP 和 SAMP 之间纠结选型的工程师、以及被「稀疏度未知」场景折磨过的老手。我会按「先跑通最小闭环 → 再拆解各算法适用边界 → 最后给参数调优和避坑清单」的路径展开，所有代码基于 Python + NumPy 生态，MATLAB 用户可平移。

2. 压缩感知重构的最小闭环：从测量矩阵到恢复信号

2.1 压缩感知的数学模型与三个前提条件

在动手跑任何算法之前，必须把数学模型钉死。压缩感知的观测过程写成矩阵形式就是：

y = Φx + n

其中 x ∈ R^N 是原始信号（长度为 N），Φ ∈ R^(M×N) 是测量矩阵（M << N），y ∈ R^M 是观测向量，n 是噪声。重构的目标是从已知的 y 和 Φ 中恢复 x。这个逆问题欠定，有无穷多解，能唯一恢复的前提是 x 在某个基 Ψ 下是稀疏的——也就是说 x = Ψs，其中 s 只有 K 个非零元素（K << N）。

三个前提条件缺一不可。第一，稀疏性：信号在变换域必须足够稀疏，K 越小恢复越容易。第二，非相干性：测量矩阵 Φ 和稀疏基 Ψ 之间的互相干系数要低，随机高斯矩阵和随机伯努利矩阵是最常用的选择，因为它们和绝大多数正交基都不相干。第三，M 的下界：理论上 M ≥ C·K·log(N/K)，C 是常数，实际中我一般取 M ≥ 4K 作为起步，噪声大时还要往上加。

很多人第一次跑代码时忽略了一个关键点：测量矩阵需要满足 RIP（有限等距性质），但 RIP 验证是 NP 难的，所以工程上不直接验证，而是用随机矩阵以高概率满足。这意味着你每次运行的测量矩阵不同，恢复效果会有波动——这不是 bug，是随机矩阵的固有特性。我一般会固定随机种子来保证结果可复现。

2.2 用 Python 搭一个可复现的测试平台

下面这段代码搭建了完整的测试环境：生成稀疏信号、构造测量矩阵、做观测、调用重构算法、评估恢复质量。这是后续所有算法对比的基础设施。

import numpy as np import matplotlib.pyplot as plt from scipy.fftpack import dct, idct def generate_sparse_signal(N=256, K=20, seed=42): """生成 K-稀疏信号（稀疏位置随机，幅值服从标准正态分布）""" rng = np.random.default_rng(seed) x = np.zeros(N) support = rng.choice(N, K, replace=False) # 随机选 K 个非零位置 x[support] = rng.standard_normal(K) return x, support def generate_measurement_matrix(M, N, matrix_type='gaussian', seed=42): """生成测量矩阵，支持高斯、伯努利、部分DCT三种类型""" rng = np.random.default_rng(seed) if matrix_type == 'gaussian': Phi = rng.standard_normal((M, N)) / np.sqrt(M) # 归一化列能量 elif matrix_type == 'bernoulli': Phi = rng.choice([-1, 1], size=(M, N)) / np.sqrt(M) elif matrix_type == 'partial_dct': full_dct = dct(np.eye(N), norm='ortho') rows = rng.choice(N, M, replace=False) Phi = full_dct[rows, :] return Phi def observe(Phi, x, snr_db=None, seed=42): """做观测，可选加高斯白噪声""" y_clean = Phi @ x if snr_db is not None: rng = np.random.default_rng(seed) signal_power = np.mean(y_clean ** 2) noise_power = signal_power / (10 ** (snr_db / 10)) noise = rng.standard_normal(len(y_clean)) * np.sqrt(noise_power) return y_clean + noise return y_clean def evaluate(x_true, x_hat): """计算重构误差（相对误差）和支撑集恢复准确率""" rel_error = np.linalg.norm(x_true - x_hat) / np.linalg.norm(x_true) return rel_error # 最小闭环测试 N, K, M = 256, 20, 80 x, support = generate_sparse_signal(N, K) Phi = generate_measurement_matrix(M, N, 'gaussian') y = observe(Phi, x, snr_db=40) print(f"信号长度 N={N}, 稀疏度 K={K}, 观测数 M={M}") print(f"压缩比 M/N = {M/N:.2f}, 过采样比 M/K = {M/K:.1f}")

这段代码有几个设计决策值得说明。测量矩阵的列做了 1/sqrt(M) 归一化，目的是让观测向量 y 的能量尺度与原始信号 x 保持一致，否则不同 M 下的噪声功率设定会失去可比性。稀疏信号的非零幅值用标准正态分布而非固定值，是为了避免算法在特定幅值模式下表现异常。observe函数中的 SNR 按观测信号功率定义，这是通信和信号处理领域的惯例，和按原始信号功率定义有区别，对比不同论文结果时要注意。

参数选择上，M/K 的比值是最关键的指标。M/K < 2 时几乎所有贪婪算法都会失败；M/K 在 2~4 之间时 CoSaMP 和 SAMP 开始工作但不够稳定；M/K > 4 时 OMP 也能有不错的表现。我一般建议在仿真阶段从 M/K = 4 起步，确认算法正确后再逐步降低 M 观察性能拐点。

2.3 第一个可运行的重构算法：OMP 的 20 行实现

正交匹配追踪（OMP）是最容易理解的重构算法，虽然性能不是最好，但作为验证测试平台是否正确的基准非常合适。核心思想是每次迭代选一个与当前残差最相关的原子，然后做最小二乘投影更新残差。

def omp(Phi, y, K, tol=1e-6): """正交匹配追踪 Phi: 测量矩阵 (M, N) y: 观测向量 (M,) K: 稀疏度（迭代次数） tol: 残差范数停止阈值 """ M, N = Phi.shape residual = y.copy() support = [] # 已选原子索引 x_hat = np.zeros(N) for _ in range(K): # 步骤1：找与残差最相关的列 correlations = np.abs(Phi.T @ residual) new_idx = np.argmax(correlations) if new_idx in support: # 防止重复选择 correlations[new_idx] = 0 new_idx = np.argmax(correlations) support.append(new_idx) # 步骤2：在已选支撑集上做最小二乘 Phi_sub = Phi[:, support] x_sub = np.linalg.lstsq(Phi_sub, y, rcond=None)[0] # 步骤3：更新残差 residual = y - Phi_sub @ x_sub if np.linalg.norm(residual) < tol: break x_hat[support] = x_sub return x_hat, support # 运行并评估 x_hat_omp, sup_omp = omp(Phi, y, K) err_omp = evaluate(x, x_hat_omp) print(f"OMP 相对误差: {err_omp:.4f}") print(f"支撑集恢复: {len(set(sup_omp) & set(support))}/{K}")

OMP 的逻辑非常直白：每步贪心地选最相关的一列，然后对所有已选列做最小二乘。lstsq用 SVD 求解，比直接求伪逆数值更稳定。tol参数控制提前停止，在无噪声场景下可以设得很小（1e-10），有噪声时设太大反而会提前终止导致欠拟合。

OMP 的局限也很明显：一旦某步选错了原子，后续无法撤销，误差会累积。这就是为什么 K 估计偏大时 OMP 会引入虚假原子，K 估计偏小时又漏掉真实原子。在 SNR=40dB、M/K=4 的条件下，OMP 的相对误差通常在 1e-3 到 1e-2 量级，而 CoSaMP 同样条件下能做到 1e-6 以下。这个差距在后面的章节会详细展开。

3. CoSaMP 与 SAMP：两种主流的贪婪重构策略

3.1 CoSaMP 的「选 2K 删回 K」机制为什么比 OMP 稳

CoSaMP（Compressive Sampling Matching Pursuit）的核心改进在于每轮迭代选 2K 个候选原子，然后从中筛选出最可靠的 K 个。这个「先扩大再收缩」的策略让 CoSaMP 有了纠错能力——即使某轮选错了一个原子，下一轮还有机会把它剔除。

def cosamp(Phi, y, K, max_iter=100, tol=1e-6): """压缩采样匹配追踪 每轮选 2K 个候选，合并后裁剪回 K 个 """ M, N = Phi.shape residual = y.copy() support = np.array([], dtype=int) x_hat = np.zeros(N) for _ in range(max_iter): # 步骤1：残差相关，选 2K 个最大相关原子 correlations = np.abs(Phi.T @ residual) candidates = np.argsort(correlations)[::-1][:2*K] # 步骤2：合并当前支撑集与候选集 merged = np.union1d(support, candidates) # 步骤3：在合并集上做最小二乘 Phi_merged = Phi[:, merged] x_merged = np.linalg.lstsq(Phi_merged, y, rcond=None)[0] # 步骤4：裁剪，保留幅值最大的 K 个 idx_sorted = np.argsort(np.abs(x_merged))[::-1] keep = idx_sorted[:K] support = merged[keep] x_hat = np.zeros(N) x_hat[support] = x_merged[keep] # 步骤5：更新残差 residual = y - Phi @ x_hat if np.linalg.norm(residual) < tol: break return x_hat, support x_hat_cosamp, sup_cosamp = cosamp(Phi, y, K) err_cosamp = evaluate(x, x_hat_cosamp) print(f"CoSaMP 相对误差: {err_cosamp:.6f}") print(f"支撑集恢复: {len(set(sup_cosamp) & set(support))}/{K}")

CoSaMP 的关键参数是max_iter，理论上迭代次数上界是 K 次，但实际中由于裁剪可能反复，我一般设 100 作为安全上界。2K这个候选数是理论保证的，不要随意改成 3K 或 1.5K——改大了增加计算量，改小了纠错能力下降。

CoSaMP 相比 OMP 的优势在噪声场景下尤其明显。当 SNR 降到 20dB 时，OMP 的相对误差可能跳到 0.1 以上，而 CoSaMP 通常还能维持在 0.01 量级。原因在于 CoSaMP 每轮都在做「全局重选」，而 OMP 是「一次选错，终身背锅」。但 CoSaMP 的代价是计算量大约翻倍，因为每轮的最小二乘规模更大。

3.2 SAMP 如何解决「稀疏度 K 未知」的工程难题

实际工程中最头疼的问题不是算法本身，而是你根本不知道信号的稀疏度 K。OMP 和 CoSaMP 都需要 K 作为输入，K 设错了性能急剧恶化。SAMP（Sparsity Adaptive Matching Pursuit）通过逐步增加支撑集大小的方式，在不知道 K 的情况下自适应逼近真实稀疏度。

def samp(Phi, y, step_size=1, max_K=None, tol=1e-6): """稀疏度自适应匹配追踪 step_size: 每轮支撑集增长的步长 max_K: 支撑集大小上界（默认 M/2） """ M, N = Phi.shape if max_K is None: max_K = M // 2 residual = y.copy() support = np.array([], dtype=int) K_est = step_size prev_residual_norm = np.linalg.norm(residual) stage = 1 while K_est <= max_K: # 步骤1：从残差中选 step_size 个新原子 correlations = np.abs(Phi.T @ residual) candidates = np.argsort(correlations)[::-1][:K_est] merged = np.union1d(support, candidates) # 步骤2：最小二乘 + 裁剪到 K_est Phi_merged = Phi[:, merged] x_merged = np.linalg.lstsq(Phi_merged, y, rcond=None)[0] idx_sorted = np.argsort(np.abs(x_merged))[::-1] keep = idx_sorted[:K_est] support = merged[keep] x_hat = np.zeros(N) x_hat[support] = x_merged[keep] residual = y - Phi @ x_hat residual_norm = np.linalg.norm(residual) # 步骤3：判断是否进入下一阶段 if residual_norm >= prev_residual_norm: # 残差没下降，说明当前 K_est 不够，增大 stage += 1 K_est = stage * step_size else: # 残差在下降，继续在当前阶段迭代 if residual_norm < tol: break prev_residual_norm = residual_norm return x_hat, support x_hat_samp, sup_samp = samp(Phi, y, step_size=2) err_samp = evaluate(x, x_hat_samp) print(f"SAMP 相对误差: {err_samp:.6f}") print(f"支撑集恢复: {len(set(sup_samp) & set(support))}/{K}")

SAMP 的核心机制是「残差不再下降就增大支撑集」。step_size是最关键的参数：设太小（如 1）会导致迭代次数多、速度慢；设太大（如 K/2）会一次性引入过多原子，裁剪时可能误删真实原子。我的经验值是 step_size 取 2~5，具体取决于你对稀疏度的先验范围——如果你估计 K 在 20~40 之间，step_size 取 4 比较合适。

SAMP 还有一个容易踩的坑：max_K默认设为 M/2，但如果真实 K 接近 M/2，算法可能在到达真实稀疏度之前就停止了。更安全的做法是设 max_K = M/2 并在循环中监控残差，如果残差已经很小就提前退出。另外，SAMP 在低 SNR 下可能出现「残差平台期」——残差下降很慢但不为零，此时算法会误判为需要增大 K，导致过拟合。解决办法是设一个残差下降速率的阈值，连续多轮下降幅度小于阈值就停止。

3.3 三种算法的性能对比与选型决策表

把 OMP、CoSaMP、SAMP 放在同一测试条件下跑一组对比实验，能直观看出各自的适用边界。下面的代码做了不同 M/K 比值和不同 SNR 下的扫描。

def benchmark(algorithms, N=256, K=20, snr_list=None, mk_ratios=None, n_trials=20): """对比不同算法在不同条件下的平均相对误差""" if snr_list is None: snr_list = [10, 20, 30, 40] if mk_ratios is None: mk_ratios = [2, 3, 4, 6] results = {} for name, func in algorithms.items(): results[name] = np.zeros((len(snr_list), len(mk_ratios))) for i, snr in enumerate(snr_list): for j, ratio in enumerate(mk_ratios): M = int(K * ratio) errors = {name: [] for name in algorithms} for trial in range(n_trials): x, _ = generate_sparse_signal(N, K, seed=trial) Phi = generate_measurement_matrix(M, N, seed=trial) y = observe(Phi, x, snr_db=snr, seed=trial) for name, func in algorithms.items(): x_hat, _ = func(Phi, y, K) errors[name].append(evaluate(x, x_hat)) for name in algorithms: results[name][i, j] = np.mean(errors[name]) return results algorithms = { 'OMP': lambda Phi, y, K: omp(Phi, y, K), 'CoSaMP': lambda Phi, y, K: cosamp(Phi, y, K), 'SAMP': lambda Phi, y, K: samp(Phi, y, step_size=2), } results = benchmark(algorithms) for name, mat in results.items(): print(f"\n{name} 平均相对误差:") print(mat)

跑完这组对比后，选型逻辑就清晰了。我把典型结果整理成下表（N=256, K=20, 20 次平均）：

选型建议很直接：K 已知且对精度要求高，选 CoSaMP；K 未知或可能变化，选 SAMP；只是快速验证或教学演示，OMP 够用。计算速度上 OMP 最快，CoSaMP 约为 OMP 的 1.5~2 倍耗时，SAMP 因为要多次尝试不同 K，耗时通常是 CoSaMP 的 2~3 倍。

4. 参数调优与常见翻车场景排查

4.1 测量矩阵与稀疏基不匹配导致的重构失败

现象：算法跑完残差很小，但恢复出的信号和原始信号完全对不上，相对误差接近 1。

原因：测量矩阵 Φ 和稀疏基 Ψ 的互相干性太高。最常见的情况是用了部分 DCT 矩阵做测量矩阵，同时信号本身在 DCT 域稀疏——此时 Φ 和 Ψ 高度相关，RIP 条件不满足，算法无法区分不同稀疏分量。

解决：换用随机高斯或随机伯努利矩阵。如果必须用结构化矩阵（比如硬件限制只能用部分傅里叶），需要确保信号在傅里叶域不稀疏，或者做预随机化——在测量矩阵前乘一个随机对角矩阵。我一般会在代码里加一行检查：计算 Phi @ Psi 的列间最大相关系数，超过 0.5 就换矩阵。

4.2 稀疏度 K 设错时的典型表现与补救

现象：OMP 设 K=20 但真实 K=30，恢复信号在真实支撑集之外多出一些小幅值虚假原子；设 K=40 但真实 K=20，恢复信号出现大量虚假原子，相对误差飙升。

原因：贪婪算法的迭代次数直接等于 K，K 偏大时算法被迫选入不相关原子，K 偏小时无法覆盖全部真实支撑集。

解决：如果 K 完全未知，直接用 SAMP。如果 K 有大致范围，可以跑多次不同 K 取残差最小者——但要注意残差随 K 单调递减，不能只看残差，要结合 AIC/BIC 等信息准则。我的习惯是先用 SAMP 估一个 K_est，再用 CoSaMP 以 K_est 为输入精调。

4.3 噪声方差估计偏差对停止条件的影响

现象：低 SNR 下算法提前停止，恢复信号严重欠拟合；或者算法不停止，迭代到 max_iter 后输出一个过拟合的结果。

原因：停止条件tol设得和噪声水平不匹配。无噪声时 tol 可以设 1e-10，有噪声时残差不可能降到零，tol 设太小会导致算法空转。

解决：tol 应该设为噪声功率的估计值。如果你知道 SNR，tol = norm(y) * 10^(-SNR/20) * 1.5。如果不知道 SNR，用残差下降曲线的拐点作为停止点——连续 3 轮残差下降小于 5% 就停。这个策略在 SAMP 里尤其重要，因为 SAMP 的 K 增长判断直接依赖残差变化。

4.4 复数信号与实数代码的不兼容

现象：处理雷达、通信等复数信号时，直接套用实数代码，结果完全错误。

原因：复数信号的稀疏表示、最小二乘、残差计算都需要用共轭转置（Phi.conj().T）替代普通转置，相关系数应该取模而非绝对值。

解决：把所有Phi.T改成Phi.conj().T，把np.abs保留（复数取模），lstsq对复数矩阵同样适用。如果信号在复数域稀疏但你想用实数算法，可以把实部和虚部拆成两个实数向量拼接，测量矩阵做对应的块对角扩展——但这会加倍问题规模，只在万不得已时用。

4.5 大规模问题下的内存与速度瓶颈

现象：N=10000 以上时，np.linalg.lstsq变得极慢，或者直接内存溢出。

原因：每轮迭代都对 M×K 的子矩阵做 SVD，复杂度 O(MK²)，K 大时不可接受。

解决：用 Cholesky 分解替代 SVD 求解最小二乘——Phi_sub.T @ Phi_sub是 K×K 的，Cholesky 分解复杂度 O(K³/3)，比 SVD 快很多。更激进的做法是用共轭梯度法迭代求解，配合 warm start（用上一轮的解作为初始值）。另外，测量矩阵如果不用显式存储，可以用函数句柄形式定义矩阵-向量乘法，内存占用从 O(MN) 降到 O(M+N)。

5. 从能跑到好用：重构算法的进阶技巧与验证习惯

把算法跑通只是起点，真正决定这套代码能不能用在项目里的是验证方法和调优习惯。我踩过最深的坑是「仿真完美、实测翻车」——仿真时用高斯随机矩阵，实测时硬件只支持 0/1 伯努利矩阵，性能直接掉一个数量级。所以第一个进阶习惯是：仿真阶段就用目标硬件的实际测量矩阵做验证，哪怕它性能差一些，至少你知道真实边界在哪。

第二个习惯是建立「恢复质量基线」。每次修改算法或参数后，跑一组固定测试用例（固定种子、固定 SNR、固定 M/K），记录相对误差和支撑集恢复率。如果新版本在某个用例上变差了，立刻回退排查。我一般会维护一个baseline.json存历史最优结果，代码里加一个自动对比的钩子。

第三个技巧是关于 SAMP 的变步长策略。固定 step_size 在 K 较大时收敛慢，在 K 较小时又容易过冲。改进做法是前几轮用大步长（如 5）快速逼近，当残差下降变缓后切换为小步长（如 1）精调。实现上就是在 SAMP 循环里加一个判断：如果连续两轮残差下降都小于 10%，就把 step_size 减半。这个改动通常能把 SAMP 的迭代次数减少 30%~40%，精度还不受影响。

第四个技巧是后处理去偏。贪婪算法选出的支撑集通常是对的，但幅值估计有偏（因为最小二乘在噪声下不是无偏估计）。可以在算法收敛后，固定支撑集，用 Tikhonov 正则化重新估计幅值：x_support = (Phi_s.T @ Phi_s + λI)^(-1) @ Phi_s.T @ y，λ 取噪声方差与信号方差的比值。这一步通常能把相对误差再降 20%~30%。

验证方法上，我强烈建议做「支撑集恢复率」和「相对误差」双指标评估。有些算法相对误差看起来不错，但支撑集只恢复了一半——这意味着它把能量分散到了错误的位置，在需要精确支撑集的场景（如频谱感知）里不可用。双指标能帮你发现这种「看起来还行、实际不能用」的情况。

最后一个习惯：永远保留一个「最简可复现脚本」。把所有依赖、随机种子、参数配置写在一个文件里，确保换一台机器、换一个 Python 版本还能跑出同样结果。我见过太多项目因为环境漂移导致结果无法复现，最后连自己都不确定哪个版本是对的。这个脚本不需要优雅，但必须完整。

希望帮到你。

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