简介基于压缩感知的二维 OMP 算法 MATLAB 实现包面向图像处理、医学成像、遥感与通信等方向的研究者和工程师用于从低采样率观测中重构二维图像或矩阵信号。压缩包内只有一个核心文件 OMP2D.m资源大小仅 3KB但完整覆盖二维 OMP 的初始化、残差计算、原子选择、系数更新与支撑集迭代流程可适配 DCT 等基矩阵输出重构信号、支撑集与系数向量直接用于 MRI 重建、图像压缩、去噪等典型场景。目前已有 555 人学习下载体现出在稀疏恢复领域的一定参考价值。获取后可快速掌握算法实现思路通过修改基矩阵、阈值或最大迭代次数拓展至自己的应用中也可作为教学和项目开发的基础代码。1. 二维信号直接套一维 OMP为什么总是慢半拍做图像重建、雷达成像或者地震数据处理的人,大概率都遇到过这个场景:手里拿到的是一个 m×n 的观测矩阵,心里清楚它在某个变换域是稀疏的,习惯性操作是把它拉成一个 mn×1 的长向量,然后丢给标准的 OMP(Orthogonal Matching Pursuit)去解。如果信号尺寸是 256×256,字典原子是 65536 维,单次内积计算就要做 65536 次乘法,迭代 K 次就是 K 次全字典扫描——这个计算量对实时处理来说基本是灾难。更麻烦的是,把二维结构强行展平会破坏图像本身的行列相关性,稀疏表示系数变得不再那么“稀疏”,重建质量反而下降。OMP2D 的核心思路就是不做向量化,直接在二维矩阵空间里做匹配追踪。它利用 Kronecker 积字典的可分离结构,把“在 65536 维空间里找原子”拆成“在 256 维行空间和 256 维列空间里各找一次”,计算量从 O((mn)²) 量级降到 O(m²n mn²) 量级。这篇文章就沿着 compressed sensing 的理论框架,把二维 OMP 的实现路径、参数设置和排错经验完整讲透,代码可以直接抄。2. 从向量化到 Kronecker 结构:OMP2D 的运行机制与计算边界2.1 二维压缩感知的数学模型:为什么是双侧投影在 compressed sensing 的标准模型里,一维观测方程写为 y Φx,其中 x ∈ R^N 是稀疏信号,Φ ∈ R^(M×N) 是观测矩阵。到了二维场景,信号本身是矩阵 X ∈ R^(M×N),观测过程变成了Y Φ_l X Φ_r^T E这里 Φ_l ∈ R^(m×M) 是左侧观测矩阵,作用在 X 的行方向;Φ_r ∈ R^(n×N) 是右侧观测矩阵,作用在列方向;E 是噪声。这个双线性观测模型的关键在于:它保持了信号的行列结构信息,而不是像向量化那样把行和列的语义混在一起。如果非要向量化,需要把模型改写成 vec(Y) (Φ_r ⊗ Φ_l) vec(X),其中 ⊗ 是 Kronecker 积。数学上完全等价,但实际操作中 Φ_r ⊗ Φ_l 是一个 (mn)×(MN) 的稠密矩阵,256×256 的图像对应 65536×65536 的矩阵,存储就要 32GB(双精度浮点),根本没法算。OMP2D 的聪明之处在于:它全程不显式构造 Kronecker 积,而是利用恒等式 (A ⊗ B) vec(X) vec(BXA^T) 来做隐式计算。2.2 字典的可分离假设:OMP2D 能成立的前提二维 OMP 能工作的前提是字典具有可分离结构,即原子可以写成两个低维原子的 Kronecker 积:ψ_k,l ψ_row,k ⊗ ψ_col,l换句话说,整个字典矩阵 D D_row ⊗ D_col,其中 D_row ∈ R^(M×K_row),D_col ∈ R^(N×K_col)。很多常见的变换字典天然满足这个性质:二维离散余弦变换(DCT)就是行 DCT 和列 DCT 的 Kronecker 积;二维小波变换同样如此;甚至 learned dictionary 也可以强制施加这种结构约束去做训练。提示:如果信号本身不满足字典可分离假设(比如纹理高度各向异性的图像),OMP2D 的稀疏表示能力会明显下降。这时要么换用全字典的 2D-KSVD,要么老老实实做向量化。2.3 算法流程:每次迭代同时选一组行原子和列原子OMP2D 的每次迭代做两件事:在行方向选一个原子,在列方向选一个原子,然后联合更新系数。用数学语言描述,给定当前残差 R,候选原子对的匹配程度定义为:score(i, j) |(d_row,i)^T R d_col,j|这个 score 同时度量了行原子 i 和列原子 j 对残差的联合贡献。算法维护两个索引集合 I 和 J,分别记录被选中的行原子和列原子,每次迭代挑出使 score 最大的 (i*, j*) 加入集合,然后在张量积子空间 span{d_row,i ⊗ d_col,j : i∈I, j∈J} 上做正交投影更新。这里有个有意思的细节:一次迭代同时选两个方向的原子,支撑集的势每次加 2(Kronecker 积的秩加 1),而不是加 1。所以如果目标稀疏度是 K,迭代次数大约只需要 K/2 次。收敛速度的差异在支撑集规模较大时非常明显。2.4 计算复杂度对比:OMP2D 的三个优势区间用具体数字说话。设观测矩阵维度 mn256,字典行数和列数 MN512,稀疏度 K20,迭代次数 T10。向量化 OMP 的每次迭代,计算所有原子与残差的内积需要 O(mn × MN) 256×256×512×512 ≈ 1.7×10^10 次乘法,这还只是投影匹配那一步,还没算最小二乘更新。OMP2D 的匹配步骤拆成两步:R 与 D_row 的乘积是 O(mM × n) 256×512×256 ≈ 3.4×10^7 次乘法,再与 D_col 乘一次又是同样的量级。总计约 6.8×10^7 次,比向量化版本快了 250 倍。这个差距在大字典场景下只会拉得更大。OMP2D 的优势区间有三类。第一类是图像类信号,行列结构天然有意义;第二类是字典规模大但观测维度小的场景,比如 mn 远小于 mn 时,双侧投影比全空间搜索省太多;第三类是实时性要求高的场景,迭代次数少且每次迭代计算量可控。3. 用 Python 零依赖复现二维 OMP 算法(附可跑代码)3.1 最小可运行版本:30 行实现核心逻辑先给一个不依赖任何第三方稀疏求解库、只用 NumPy 的最小实现。这个版本优先保证逻辑清晰和数据结构的直观性,性能优化留到后面小节。import numpy as np def omp2d(Y, D_row, D_col, K, tol1e-6): 二维正交匹配追踪 Y: (m, n) 观测矩阵 D_row: (m, M) 行方向字典,每列是一个原子 D_col: (n, N) 列方向字典 K: 目标稀疏度(实际支撑集原子对数量) tol: 停止迭代的残差阈值 rows, cols Y.shape # 对字典原子做归一化,避免尺度影响原子选择 D_row D_row / np.linalg.norm(D_row, axis0) D_col D_col / np.linalg.norm(D_col, axis0) # 初始化残差、支撑集索引、系数矩阵 R Y.copy() row_idx, col_idx [], [] # 迭代次数 K(每个原子对贡献一个稀疏系数) for t in range(K): # 双侧投影:R 先左乘 D_row.T,再右乘 D_col Proj D_row.T R D_col i, j np.unravel_index(np.argmax(np.abs(Proj)), Proj.shape) # 检查是否达到收敛容差 if np.abs(Proj[i, j]) tol: break row_idx.append(i) col_idx.append(j) # 用当前支撑集做最小二乘:解 X_hat argmin ||Y - D_row[:,I] X D_col[:,J].T||_F Dr D_row[:, row_idx] Dc D_col[:, col_idx] # 利用 Kronecker 积的伪逆恒等式:(A ⊗ B)^ A^ ⊗ B^ Dr_pinv np.linalg.pinv(Dr) # (|I|, m) Dc_pinv np.linalg.pinv(Dc) # (|J|, n) X_hat Dr_pinv Y Dc_pinv.T # (|I|, |J|) # 重建当前信号并更新残差 Y_recon Dr X_hat Dc.T R Y - Y_recon # 检查残差能量是否低于阈值 if np.linalg.norm(R, fro) tol * np.linalg.norm(Y, fro): break return row_idx, col_idx, X_hat, R这段代码的核心是第 17 行:Proj D_row.T R D_col。这一行同时完成了行方向和列方向的原子匹配,Proj 中第 (i, j) 个元素的绝对值越大,说明行原子 i 加上列原子 j 组成的二维原子对残差的解释能力越强。选出来的原子对索引分别存入 row_idx 和 col_idx,构成支撑集。最小二乘更新部分用了伪逆的 Kronecker 积性质:(A ⊗ B)^ A^ ⊗ B^。这意味着不需要构造完整的二维字典矩阵再求伪逆,只需要分别对 Dr 和 Dc 求伪逆,然后做两次普通矩阵乘法,计算量从 O((mn)³) 直接降到 O(|I|³ |J|³)。当支撑集规模超过 10 时,这个差距是数十倍。3.2 验证实验:合成稀疏矩阵重建下面生成一个已知稀疏的人造信号,用上面的函数验证重建效果和收敛行为。m, n 64, 64 # 观测维度 M, N 128, 128 # 字典原子数 rng np.random.default_rng(42) # 随机字典(满足可分离假设) D_row rng.standard_normal((m, M)) D_col rng.standard_normal((n, N)) # 生成稀疏系数矩阵:只有 4 个非零元素 X_true np.zeros((M, N)) inds [(10, 20), (35, 5), (80, 90), (120, 50)] for (u, v) in inds: X_true[u, v] rng.standard_normal() * 5 # 幅度放大便于观察 # 合成观测 Y D_row X_true D_col.T # 运行 OMP2D row_idx, col_idx, X_hat, R omp2d(Y, D_row, D_col, K4) print(选择行原子:, row_idx) print(选择列原子:, col_idx) print(重建相对误差:, np.linalg.norm(Y - D_row X_hat D_col.T, fro) / np.linalg.norm(Y))这段实验的预期结果是重建相对误差在 1e-8 量级(无噪声情况下),且选出的行原子和列原子与真实支撑集完全一致。输出应该看到 row_idx 35, 120, 80, 10 ,col_idx [5, 50, 90, 20]。这里需要注意:OMP2D 选原子是成对选的,所以支撑集是 (row_idx, col_idx) 的笛卡尔积而不是一一对应。如果真实信号只有 4 个独立的非零系数,理论上迭代 4 次就能精确重建。但如果在 K 次迭代中某些系数恰好被分成两对选出来,支撑集会多出一些“伪原子对”,最小二乘会把多余的系数压到近似零,不影响重建精度。3.3 和向量化 OMP 的公平对比:同一信号、同一字典为了确认 OMP2D 不只是“能跑”,还要和标准一维 OMP 做一次公平对比。做法是把同样的二维观测展平,字典换成 Kronecker 积的显式形式,然后用 sklearn 的 OrthogonalMatchingPursuit 或手写一维 OMP 来跑。from sklearn.linear_model import OrthogonalMatchingPursuit def omp1d_vec(Y, D_row, D_col, K): m, n Y.shape y Y.flatten() # (m*n,) D np.kron(D_col, D_row) # (m*n, M*N) 显式 Kronecker 积 omp OrthogonalMatchingPursuit(n_nonzero_coefsK, fit_interceptFalse) omp.fit(D, y) x omp.coef_ # (M*N,) return x.reshape(M, N) # 还原为 2D 系数矩阵 # 重建并对比 X_omp_vec omp1d_vec(Y, D_row, D_col, K4) err_vec np.linalg.norm(Y - D_row X_omp_vec D_col.T, fro) err_mat np.linalg.norm(Y - D_row X_hat D_col.T, fro) print(f向量化 OMP 误差: {err_vec:.2e}, OMP2D 误差: {err_mat:.2e})从结果看,两种方法在无噪声情况下应该达到同等的重建精度,但运行时间的差距会非常直观。在 64×64 观测、128×128 字典的配置下,向量化 OMP 构造 D 矩阵就需要 64×64×128×128×8 字节 ≈ 536MB 内存,而 OMP2D 的 D_row 和 D_col 加起来只有 128×64×2×8 ≈ 131KB。内存差 4000 倍,时间差至少一个数量级。这就是不在显式空间里做 Kronecker 展开的根本原因——不是数值精度的问题,是物理内存和计算资源根本不允许。4. 稀疏度与停止条件怎么设:OMP2D 的三个必调参数与一组对照表4.1 稀疏度 K:比一维 OMP 更容易误判的参数OMP2D 的 K 参数含义和一维 OMP 有微妙差别。一维 OMP 里 K 直接对应非零系数的个数,而 OMP2D 里每次迭代选一对原子,K 选的是原子对数量。举例来说,如果真实信号的系数矩阵 X 是满秩的(比如一个 3×3 的块状支撑,9 个非零系数),向量化 OMP 需要 K9 才能完全重建,而 OMP2D 里这个 3×3 的块支撑集恰好对应 3 个行原子和 3 个列原子的笛卡尔积。理论上只需要 3 次迭代就能精确重建,因为 3 对原子生成的张量积子空间已经覆盖了全部 9 个系数。但真实信号往往是“部分填充”的:支撑集是若干个原子对的笛卡尔积,但里面有些元素是空的。这种情况下把 K 设在列支撑大小和实际非零系数数量之间,最小二乘会把空位压成接近零,不会影响重建精度。实际调参经验是:K 设为目标稀疏度除以支撑集的平均秩,通常取 2 到 4。# 自适应 K 的一种实用策略:设一个较大的迭代上限,让停止条件来兜底 row_idx, col_idx, X_hat, R omp2d(Y, D_row, D_col, K20, tol1e-8) print(实际迭代次数:, len(row_idx)) # 可能远小于 20这种设法的好处是让算法根据残差自动停止,避免 K 猜错导致过拟合。代价是迭代次数可能比最优多几次,但 OMP2D 每次迭代的计算量本身不大,多几次迭代完全可以接受。4.2 停止容差 tol:决定“够好”的标尺tol 参数控制的是残差能量的相对下降,在代码里以两种方式生效:一是原子匹配分数阈值,二是残差范数阈值。大多数实现只用了其中一个,OMP2D 两个都用,因为二维场景下残差能量下降不是单调的——某些迭代一步能砍掉 90% 的残差,另一步可能只砍掉 2%。单一阈值容易提前停或者不停。场景推荐 tol理由无噪声合成数据1e-8 ~ 1e-10精确重建是可达的,容差设严不亏实际观测含 1% 噪声1e-3 ~ 1e-4再往下追的都是噪声分量实时处理(视频帧)1e-2 ~ 5e-3视觉效果够用即可,省迭代病态字典(原子高度相关)1e-4 ~ 1e-5容差太松支撑集选错,太紧会过拟合噪声4.3 字典归一化:最容易被忽略的前置参数代码里做了 D_row 和 D_col 的列归一化,这一步在二维场景下比一维更关键。原因在于:如果行原子和列原子的范数不统一,Proj D_row.T R D_col 里的分数 |d_row,i^T R d_col,j| 会受两个方向原子范数乘积的影响,导致选出的原子对偏向“范数大的行 范数大的列”,而不是“真正匹配残差方向”的原子。在 MATLAB 的 OMP2D 原始实现里,归一化是在字典传入前做的,但 Python 移植版本经常漏掉这一步。一个容易踩的坑是:归一化后如果再做字典学习或者在线更新,要保证新加入的原子也做相同处理,否则新旧原子的尺度不一致,匹配结果会突然漂移。4.4 一组对照实验:参数变化如何影响结果下面用合成数据展示参数对重建 PSNR 的影响,方便抄作业时做基准。import time def benchmark_omp2d(K_candidates, tol_candidates): m, n 64, 64 M, N 128, 128 rng np.random.default_rng(0) D_row rng.standard_normal((m, M)) D_col rng.standard_normal((n, N)) D_row / np.linalg.norm(D_row, axis0) D_col / np.linalg.norm(D_col, axis0) X_true np.zeros((M, N)) for _ in range(8): # 8 个稀疏系数 u, v rng.integers(0, M), rng.integers(0, N) X_true[u, v] rng.standard_normal() Y D_row X_true D_col.T Y_noisy Y 0.01 * rng.standard_normal(Y.shape) for K in K_candidates: for tol in tol_candidates: t0 time.time() row_idx, col_idx, X_hat, _ omp2d(Y_noisy, D_row, D_col, K, tol) elapsed time.time() - t0 recon D_row X_hat D_col.T psnr 20 * np.log10(np.max(Y) / np.linalg.norm(Y - recon, fro) * np.sqrt(Y.size)) print(fK{K:2d}, tol{tol:.0e}, 迭代{len(row_idx):2d}, PSNR{psnr:5.2f} dB, 耗时{elapsed*1000:6.2f} ms) benchmark_omp2d([4, 8, 12], [1e-2, 1e-5, 1e-8])从输出中能清楚看到:当 K 从 4 增加到 8 时,PSNR 显著上升;超过真实稀疏度后,K 继续增大 PSNR 提升变缓甚至略微下降(开始拟合噪声);而 tol 对结果的影响主要体现在迭代提前终止上,当 tol1e-2 时可能只迭代 2~3 次就停了,重建质量明显不足。最佳配置是 K8(真实稀疏度)、tol1e-5,兼顾质量和速度。5. 常见发散、慢收敛与支撑集漂移的排错路径5.1 迭代发散:残差范数不减反增的元凶一个高频问题:残差能量在迭代到一半时突然反弹。最常见的原因有两个。第一个是字典的列之间近似线性相关,导致伪逆求解时出现巨大的数值解——最小二乘系数虽然拟合了当前支撑集,但系数之间互相抵消,重建信号能量远超原始信号,残差自然反弹。排查方法是检查观测矩阵的条件数:# 检查支撑集对应字典子块的条件数 Dr_sub D_row[:, row_idx] Dc_sub D_col[:, col_idx] # 转为 Kronecker 积并求条件数 D_sub np.kron(D_sub.Dc, D_sub.Dr) # 注意:这里 D_sub 只是示例,正式代码需要正确索引 cond np.linalg.cond(D_sub) print(f支撑集字典条件数: {cond:.2e}) if cond 1e12: print(警告:支撑集字典接近奇异,最小二乘结果不可靠)第二个原因是残差更新时没有把当前支撑集的正交投影完全减去。在 OMP2D 的实现里,X_hat 是在完整支撑集张量积空间上的最小二乘解,重建信号 D_row[:, row_idx] X_hat D_col[:, col_idx].T 理论上一定在支撑空间内,但如果实现里不小心用了逐次累加的方式(每次迭代只更新新增原子对的系数,而不是全部重估),残差就会带上前面的估计偏差。正确做法是每次都重新计算全部支撑原子的系数,代码里的X_hat Dr_pinv Y Dc_pinv.T正是这一点。5.2 支撑集漂移:选错原子后怎么拉回来支撑集漂移是指前几轮选出的原子在后续迭代中被证实是次优的,但由于 OMP 类算法不允许删除原子,错误会被一直带到最后。二维场景更严重的是:行方向和列方向的错误会互相放大。行选错一个原子,列方向的匹配分数全部失真;列再选错,残差里混入的错误成分就更多,行方向的新选择也会被带偏。缓解支撑集漂移的一个实用做法是“预检”——在正式迭代前,先用低秩近似给一个热启动:# 用 SVD 做初步的行列原子筛选 U, S, Vt np.linalg.svd(Y, full_matricesFalse) # 取前 r 个主成分,用主成分方向与字典原子做一次粗匹配 r min(3, len(S)) U_r, Vt_r U[:, :r], Vt[:r, :] # 对行方向:计算 U_r 与 D_row 的余弦相似度 row_scores np.abs(U_r.T D_row).max(axis0) # 每个原子的最大匹配分数 col_scores np.abs(D_col.T Vt_r.T).max(axis1) # 只保留分数前 30% 的原子作为候选集,缩小 OMP2D 的搜索空间 row_candidates np.argsort(row_scores)[-int(0.3 * M):] col_candidates np.argsort(col_scores)[-int(0.3 * N):]这个热启动能滤掉大量明显不相关的原子,让 OMP2D 在搜索时不会选到那些“局部最优”的错误原子。代价是如果真实支撑集的主成分方向太弱,可能被误杀。一般建议 r 不超过 5,候选集保留比例不低于 20%。5.3 收敛缓慢:残差下降曲线出现长尾另一种常见现象是前几步残差降得很快,后面逐步“磨蹭”,每步只降一点点。这可能意味着信号的稀疏表示系数不是集中在少数大系数上,而是分布在大量小系数上。此时 OMP2D 会花很多次迭代去捡那些小系数,效率远不如阈值法。实际工程中我一般这样处理:设置最大迭代次数上限(默认 2×K),同时监控残差下降率,如果连续三步的下降率都小于 5%,就提前终止。具体实现:residual_energies [] for t in range(2 * K): # ... 核心迭代逻辑 ... residual_energies.append(np.linalg.norm(R, fro)) if t 3: drops [(residual_energies[i-1] - residual_energies[i]) / (residual_energies[i-1] 1e-12) for i in range(-3, 0)] if all(d 0.05 for d in drops): break这种基于下降率的停止策略比固定迭代次数更鲁棒,特别是在混合稀疏信号(少量大系数加大量小系数)的场景下,能避免后面那几十次“捡芝麻”的无效迭代。5.4 精度对标:如何验证你的 OMP2D 实现正确移植代码后第一件事,是拿已知支撑集的合成信号验证。具体做法:固定 D_row 和 D_col,随机选 5 个行原子和 5 个列原子构成支撑集,系数矩阵完全随机填充,合成 Y。然后跑 OMP2D,检查两件事:选出的行原子和列原子是否与真实支撑集一致(允许顺序不同);在无噪声情况下,重建误差是否在 1e-12 以下。标准的 OMP 实现应该满足“稀疏度等于真实非零系数数量时精确重建”的保证。如果这个基本测试都过不了,说明实现里有 bug,而不是参数问题。此时从三个地方排查:投影计算是否正确使用了转置;伪逆更新是否每次重算了全部支撑系数;残差更新用的 Y_recon 是否用了最新的完整支撑集。最常见的 bug 是把残差更新写成 R R - d_i * coef,而不是 R Y - D_support coef,前者在频率域实现没问题,在时域会累积误差。6. 利用矩阵结构再快一步:预计算留一学习与低秩系数更新建完可跑的实现,下一步是榨干 Kronecker 结构带来的所有性能红利。有三个技巧在工程中实测定有用。第一个技巧是预计算伪逆的留一更新。上面的代码里,每次迭代都np.linalg.pinv(Dr)和np.linalg.pinv(Dc),当支撑集规模从 1 长到 15 时,总计算量是 Σk³ 的量级。如果改用增量更新——在已有伪逆基础上加一个新列——复杂度只有 O(k²)。具体做法是维护支撑集字典的 QR 分解,每加一个原子做一次 Givens 旋转更新。对于输出行和列的伪逆,使用scipy.linalg.qr_insert可以做到真正 O(k²) 的更新。第二个技巧是系数矩阵的低秩近似。在很多实际信号里,支撑集虽然理论上是一个完整笛卡尔积,但系数矩阵的数值秩很低(例如图像中的平滑区域,OMP2D 选出的行原子和列原子很多,但有效秩只有 3~4)。在每次迭代后,对 X_hat 做一次薄 SVD,令其秩不超过预设值 r,可以抑制噪声放大。这个操作在字典原子高度相关的病态场景下特别有用——它实际上给最小二乘加了一个隐式的正则化。第三个技巧是在频域做投影预计算。如果 D_row 和 D_col 是 DCT 或傅里叶字典,那么 D_row.T R D_col 本质上就是二维变换本身,可以直接用 FFT 计算,复杂度从 O(mM × n nN × m) 降到 O(mn log mn)。对 1024×1024 的图像级重建,这一步能省掉 90% 的匹配时间。这就是为什么很多 compressed sensing 的实时实现偏爱变换字典而不是学习字典——学习字典的原子千变万化,没法用快速变换加速。最后提一个验证预期:当观测维度从 64 升到 256,字典原子数从 128 升到 1024 时,OMP2D 的运行时间应该近似线性增长(每次迭代 O(mM × n nN × m)),而不是像向量化 OMP 那样二次方增长。如果发现你的实现随规模增长明显超出线性,大概率是某个地方不小心构造了完整 Kronecker 矩阵,回头检查内存占用就能定位问题。把这个规模测试写进 CI,每次改动字典或观测参数后跑一遍,能防止后续优化把可分离结构弄丢。本文还有配套的精品资源点击获取