简介面向信号处理与机器学习方向的研究者和学生这份资源以稀疏表示为核心提供可直接运行的完整MATLAB代码可用于信号恢复、图像去噪、压缩感知等典型任务的算法验证与实验对比。压缩包共六个文件全部为.m源码整体仅4KB涵盖主程序、PCA预处理、l1eq_pd优化求解、样本读取与准确率计算等模块代码结构清晰方便逐段调试和按需替换字典或基函数。该资源已有6565人学习下载作为轻量级示例具有较高的参考热度。通过阅读并运行这些代码可以直观体会稀疏表示如何以少量非零系数逼近原始信号理解LASSO与BPDN等优化模型的构造和求解思路并进一步掌握字典学习对提升表示效果的作用为后续在图像复原、特征提取或深度学习中的应用提供可复现的编程基础。1. 稀疏表示MATLAB 里从零写一遍比背公式有用得多之前做图像去噪时我拿着现成的稀疏表示工具箱调出来的效果总比论文差一截后来才意识到问题出在“字典”和“稀疏度”这两个参数上。稀疏表示的本质很简单把一段信号 y 写成字典 D 与一个“大多数位置都是 0”的系数向量 x 的乘积y 可以是一维振动波形、二维图像块甚至是 CT 的投影数据。在 MATLAB 里跑它既可以直接调用 lasso也可以用几十行代码自己实现 OMP 算法。这篇文章顺着一条最小闭环讲从理论选型到完整可运行的 MATLAB 代码再到字典学习和 5 个容易翻车的坑适合做毕设、搞图像处理或者刚接触压缩感知的开发者。2. 稀疏表示的理论选型L0/L1 该解哪个OMP/Lasso 怎么选固定字典还是学习字典2.1 稀疏表示的数学问题L0 不可解L1 才是工程解稀疏表示要解决的核心问题可以写成模型y D*x e。y是 m 行观测信号D是 m×n 的字典矩阵x是 n×1 的系数向量e是噪声。当 n 远大于 m 时这个线性系统欠定解不唯一稀疏表示就是在所有可行解里找非零元素个数最少的那一个。这个“最少”用数学语言描述就是||x||_0也就是零范数即统计 x 里有多少个非零元素。严格求解 L0 范数需要遍历所有原子组合这是一个 NP 难问题。真实信号维度稍微上去一点比如字典是 256×512遍历组合就会把计算时间拖到无法接受。工程上最常见的做法是退一步把 L0 范数替换成 L1 范数。L1 是绝对值求和优化问题变成凸优化可以高效求解。MATLAB 自带的lasso函数直接面对这个目标但它在实际稀疏编码中使用时需要额外讨论正则化系数的物理意义反而没有自己写 OMP 来得直观。所以我对初学者的建议是先动手写一个 OMP它把复杂的目标函数拆成“每次挑一个原子、更新一次系数、算一次残差”的循环。OMP 不直接优化全局目标而是用贪婪策略逼近这种“不完美但有效”的思路恰恰是工程落地里最常见的。弄懂 OMP 之后再看 Lasso、ADMM 这些方法你会更容易理解它们各自在优化什么以及为什么有些问题必须换凸优化。目标范数标准问题形式MATLAB 落地方式备注L0min ||x||0s.t. ||y-Dx||2 ε自实现 OMP 类算法NP 难小规模可算L1min ||y-Dx||2^2 λ||x||1lasso、quadprog凸优化系数稳定L2最小二乘解D\y不稀疏只能当基线2.2 算法选型OMP、CoSaMP、Lasso 各适合什么数据规模很多 MATLAB 初学者拿到稀疏表示问题第一反应是搜“稀疏表示 MATLAB 代码”然后复制一个求解函数就开始跑跑不通就归结为“玄学”。其实绝大多数失败都能从算法选型找到原因。我把常用求解算法按适用场景列一张表这张表基本就是我在实际项目里的选型依据。算法优化策略典型适用场景不适用场景OMP贪心逼近 L0字典规模几千以内、稀疏度已知、追求速度字典列高度相关、噪声很强CoSaMP压缩感知专用重加权随机观测矩阵下的信号重建字典结构复杂、相关性高LassoL1 凸优化高维系数稳定性优先、样本量不大需要严格稀疏支撑集时ADMM分裂变量迭代融合 TV、核范数等复杂正则项只求快速验证的场景如果是做图像去噪或振动信号特征提取字典规模一般不超过 1000 列OMP 足够。如果观测矩阵是随机高斯矩阵、观测数远小于信号长度这类压缩感知场景CoSaMP 更稳。如果后端还要做统计推断系数要稳定而不是绝对稀疏我会优先把 Lasso 作为基线因为它有明确的正则化路径。如果问题里除了稀疏项还有全变分或低秩项那 ADMM 这类拆分变量迭代的方法更适合否则目标函数叠在一起没法直接求梯度。在多算法融合的图像处理系统里稀疏表示经常作为前级特征提取模块把信号从原始域映射到稀疏系数域后端再接分类器或重建模块。这种场景我一般会把固定字典先跑通基线再根据指标决定要不要上 K-SVD。选型不是越高级越好而是让每一步的计算代价都在可控范围。2.3 固定字典与学习字典什么时候必须放弃 DCT固定字典指的是直接用数学解析式生成的基函数集合常见有 DCT、小波、傅里叶、Gabor。这类字典构造快、可解释性强MATLAB 里一条dctmtx(8)就能生成 8×8 图像的 DCT 基。对平滑信号和规则纹理固定字典往往已经足够系数稀疏度也满足要求。但它的问题是字典结构固定无法针对特定数据集自适应遇到复杂自然图像时每个块可能需要更多原子才能达到同样的重构精度。学习字典则是从训练样本里“学”出来的。给定一批样本矩阵 Y目标是同时求解字典 D 和稀疏系数 X使得重构误差尽量小这个优化就是字典学习的核心。MOD、K-SVD 是两种最经典的学习算法它们交替执行“稀疏编码”和“字典更新”。K-SVD 的优势在于每次迭代更新一个原子同时更新对应系数行收敛更快字典质量更高。缺点是计算量明显增大尤其是训练样本几千列时MATLAB 里迭代 20 轮通常要跑几十秒到几分钟。字典类型构造速度重构精度适用场景DCT / 小波极快平滑信号够用快速验证、常规压缩过完备 Gabor一次生成即可时频特征好振动信号、语音K-SVD 学习字典慢自然图像去噪更好图像处理、特征学习所以我的选型准则是先用固定字典跑一个最低可行版本如果重构残差始终压不下去再切换到 K-SVD。很多人一上来就训练字典结果把时间耗在迭代上最后稀疏度还是没调明白反而掩盖了真正的参数问题。3. 用 MATLAB 从零实现 OMP 求解稀疏表示完整代码、参数说明与重构验证3.1 最小闭环第一步生成过完备 DCT 字典和稀疏测试信号OMP 算法的输入包含三样东西字典、观测信号、稀疏度或残差阈值。为了让代码可复现我先写一个生成过完备 DCT 字典的函数。这里采用离散余弦基原因是 DCT 与图像压缩关系紧密而且列之间的相关性比普通傅里叶基低很多。% 生成 m x n 的过完备 DCT 字典 % m: 观测维度n: 字典原子数要求 n m 才构成过完备 function D makeDCTDict(m, n) D zeros(m, n); for k 0:n-1 if k 0 % 第一个原子是直流分量归一化处理 D(:, k1) ones(m, 1) / sqrt(m); else % 其余原子是不同频率的余弦基 D(:, k1) sqrt(2/m) * cos(pi * k * (0:m-1) / m); end end % 对每一列做 L2 范数归一化避免 OMP 偏向范数大的原子 D D ./ sqrt(sum(D.^2, 1)); end参数说明m 是观测维度也就是信号的长度或者图像块的像素数n 是字典原子数必须大于 m否则系统相对确定稀疏表示就失去意义。列归一化这一步不能省否则 OMP 在选择原子时范数大的原子天然更容易被选中但选出来的原子并不一定和残差最相关。接下来生成一个稀疏测试信号用于验证整个闭环是否正确。% 生成模拟观测数据 m 64; % 观测维度 n 128; % 字典原子数 D makeDCTDict(m, n); K_true 6; % 真实稀疏度也就是非零系数个数 rng(42); % 固定随机种子保证可复现 idx sort(randperm(n, K_true)); % 随机选 6 个原子位置 x_true zeros(n, 1); x_true(idx) randn(K_true, 1); % 这 6 个位置给随机振幅 % 用字典和稀疏系数合成观测再加一点高斯噪声 y D * x_true 0.05 * randn(m, 1);这里 K_true 设成 6代表原始信号只由 6 个基函数组合而成。噪声标准差设为 0.05相对幅值不算大目的是让算法能清楚看到真实支撑集。信号生成完毕之后就可以把这组数据交给 OMP 去恢复 x看看算法能不能找到正确的位置。3.2 OMP 核心代码原子选择、最小二乘更新与残差控制OMP 全称 Orthogonal Matching Pursuit中文叫正交匹配追踪。它的流程非常固定先找与当前残差最相关的原子然后把这些原子联合起来做一次最小二乘投影用投影结果更新残差重复迭代直到满足停止条件。下面是完整实现核心部分不超过 30 行。% 正交匹配追踪 OMP % 输入: % D - 字典矩阵m x n列必须已归一化 % y - 观测向量m x 1 % K - 稀疏度上限即最多允许选择 K 个原子 % tol - 残差相对阈值默认 1e-4 % 输出: % x - 稀疏系数向量n x 1 function x OMP(D, y, K, tol) if nargin 4, tol 1e-4; end [~, n] size(D); x zeros(n, 1); % 系数向量初始化为全零 r y; % 残差初始化为观测信号 selected []; % 已选原子索引集合 for t 1:K % 1. 计算所有原子与残差的内积绝对值 % 内积绝对值越大说明该原子与残差的方向越接近 proj abs(D * r); % 2. 屏蔽已经选过的原子避免重复选择同一个基 proj(selected) -inf; [~, j] max(proj); selected [selected, j]; % 3. 对已选原子构成的子字典做最小二乘 % x_sub 是在当前支撑集下的最优系数 Ds D(:, selected); x_sub Ds \ y; % 4. 用新系数重建信号更新残差 r y - Ds * x_sub; % 5. 如果残差相对能量已经低于阈值提前终止 if norm(r) / norm(y) tol break; end end % 把支撑集上的系数填回完整系数向量 x(selected) x_sub; end每个步骤说明第 1 步的D * r是一个 n 维向量它计算字典所有列与残差的内积物理意义是每个原子与当前残差的“匹配度”。第 3 步的Ds \ y是 MATLAB 里的最小二乘求解因为Ds不一定方阵用反斜杠运算符会比显式写伪逆更稳定。第 5 步的停机条件是残差能量与原始信号能量的比值当它小于 tol 时认为已经收敛。proj(selected) -inf这一段非常关键去掉它你会在迭代日志里看到同一个原子被反复选中的情况。3.3 重构效果怎么验证SNR、残差曲线与稀疏度观察代码写完不算完成必须验证恢复出来的系数是否正确。最直接的办法是把重构信号和原始观测放到一张图里对比同时打印支撑集信息。下面是验证代码。% 调用 OMP 恢复系数 x_hat OMP(D, y, K_true 2, 1e-4); % 计算重构信号和信噪比 y_rec D * x_hat; residual y - y_rec; snr 20 * log10(norm(y) / norm(residual)); fprintf(真实支撑集索引: %s\n, mat2str(idx)); fprintf(OMP 找到的支撑集: %s\n, mat2str(find(x_hat ~ 0))); fprintf(重构信噪比: %.2f dB\n, snr); % 画对比图 figure; subplot(2,1,1); plot(y, linewidth, 1.2); title(带噪观测信号); subplot(2,1,2); plot(y_rec, linewidth, 1.2); title(OMP 重构信号);这里我把稀疏度上限设为 K_true 2也就是 8。这么做是为了模拟真实场景里我们并不知道确切的稀疏度只能给一个略高于真实值的上限。如果恢复成功打印出的支撑集索引应该和真实索引完全一致信噪比通常在几十分贝量级。如果支撑集对不上就需要回头检查字典是否归一化、噪声是否过大、tol 设置是否合理。通过这个最小闭环你就已经理解稀疏表示的核心链路字典构造、稀疏编码、残差判断。下一步是把这个流程升级为“字典也可以学习”的完整方案。4. 把稀疏编码升级成字典学习K-SVD 的 MATLAB 实现与收敛性检查4.1 为什么固定字典会在真实数据上“现原形”固定字典在模拟信号上表现不错但换成真实图像块时问题很快暴露。举个例子cameraman 图像里既有平滑的天空区域又有复杂纹理的人物边缘DCT 字典对平滑块编码效率高对纹理块则需要大量原子稀疏度指标很难好看。更关键的是固定字典没有“记住”训练数据的结构它只是从数学上保证基函数完备而不是针对当前数据最优。K-SVD 的做法是连字典一起训练。输入是一组样本矩阵 Y每个样本是一列输出是学到的字典 D 和对应稀疏系数 X。优化目标是最小化整体重构误差||Y - D*X||_F^2同时约束 X 的每一列稀疏度不超过设定值。这个优化不是一步到位的而是拆成两个阶段反复交替先用 OMP 固定 D 求 X再固定 X 更新 D。这种交替优化在工程里非常常见稀疏表示只是其中一个典型应用。4.2 K-SVD 的 MATLAB 实现稀疏编码阶段与字典更新阶段下面是我常用的 K-SVD 函数实现。为了让代码可读我把前面写的 OMP 直接作为子函数调用。% K-SVD 字典学习 % 输入: % Y - 训练样本矩阵m x L每列一个样本 % K - 想让字典含有的原子数 % iter - 最大迭代次数 % tol - OMP 稀疏编码时的残差阈值 % 输出: % D - 学习得到的字典m x K % X - 稀疏系数矩阵K x L % history - 每轮重构误差曲线 function [D, X, history] KSVD(Y, K, iter, tol) [m, L] size(Y); % 1. 初始化随机从训练样本里抽 K 列作为初始原子 idx randperm(L, K); D Y(:, idx); D D ./ sqrt(sum(D.^2, 1) 1e-12); % 加微小常数防止除零 history zeros(iter, 1); for it 1:iter % --- 稀疏编码阶段 --- % 固定字典 D对每个训练样本单独用 OMP 求稀疏系数 X zeros(K, L); for j 1:L X(:, j) OMP(D, Y(:, j), K, tol); end % --- 字典更新阶段 --- % 逐个原子更新每个原子只影响使用了它的那些样本 for k 1:K % 找到系数矩阵第 k 行非零的样本也就是用过该原子的样本 use_idx find(X(k, :) ~ 0); if isempty(use_idx) continue; % 这个原子没有被任何样本使用先跳过 end % 计算去除第 k 个原子后的残差矩阵 % 注意要先把当前原子贡献加回来否则会重复扣减 E_k Y(:, use_idx) - D * X(:, use_idx) D(:, k) * X(k, use_idx); % 对残差矩阵做 SVD取主奇异向量更新原子和系数 [U, S, V] svd(E_k, econ); D(:, k) U(:, 1); % 新原子取左奇异向量 X(k, use_idx) S(1, 1) * V(:, 1); % 对应系数行更新 end % 记录本轮整体重构误差用于判断是否收敛 history(it) norm(Y - D * X, fro) / norm(Y, fro); if it 1 abs(history(it) - history(it-1)) 1e-6 * history(it-1) break; end end end逻辑说明稀疏编码阶段本质上是把前面写的 OMP 跑 L 遍每一遍对应一个训练样本。字典更新阶段是 K-SVD 区别于其他方法的关键它不是一次更新整本字典而是逐列更新。E_k的含义是“去掉第 k 个原子之后那些使用过它的样本还剩多少残差”。对这个残差矩阵做 SVD取主奇异向量作为新原子这样使得单步更新在 Frobenius 范数意义下最优。S(1,1) * V(:,1)这个写法把奇异值合并到系数行里可以保证更新后重构误差不增。千万别直接用E_k的第一列作为新原子那相当于把一个样本当成了字典原子忽略其他样本信息。K-SVD 的核心恰恰是利用 SVD 从残差矩阵里提取“所有相关样本共享的主成分”。4.3 K-SVD 参数配置字典大小、迭代次数与收敛检查K-SVD 能调的参数不算多但每一个影响都很大。我把常用参数的推荐范围列成表方便你直接当模板使用。参数推荐范围对结果的影响配置注意事项原子数 K64 ~ 256原子越多字典表达能力越强不要超过训练样本数的一半迭代次数 iter10 ~ 30前期误差下降明显超过 30 轮基本没有收益OMP 阈值 tol1e-4 ~ 1e-5控制每个样本用几个原子太大会欠拟合太小会过拟合字典初始化随机抽样 / DCT 基随机抽样会改变每次结果建议固定随机种子便于复现收敛检查我做两件事。第一观察history曲线正常情况下应该单调下降且前 5 轮下降幅度最大后面趋于平缓。第二检查有没有“死原子”也就是在字典更新阶段永远找不到use_idx的原子。如果一个原子从头到尾没有被任何样本使用说明它游离在有效表示空间之外可以直接删掉也可以重新用残差最大的样本去初始化它。如果重构误差峰值一直降不下来先确认稀疏编码阶段用的K是不是太小了。OMP 里传入的稀疏度上限如果低于训练样本真实稀疏度残差就会卡在某个水平后面字典再训也没用。反过来如果tol设得太小每个样本都会被塞入很多原子字典会偏向过拟合泛化能力变差。参数之间存在耦合我最常用的顺序是先把 tol 固定在 1e-4小范围扫 K 和迭代次数找到最佳区间后再回来微调 tol。5. 稀疏表示实战避坑5 条血泪记录与排查方向5.1 字典列没有归一化OMP 提前退出现象OMP 迭代一两轮之后残差就不再下降程序直接退出输出支撑集和真实支撑集完全对不上甚至出现空支撑集。原因字典某些列的 L2 范数远大于其他列内积绝对值排序时大范数原子天然靠前OMP 反复选择它们但真正与残差方向吻合的小范数原子反而没机会被选中。更隐蔽的是如果字典里存在全零列D * r会产生 NaNmax 函数直接失效。解决字典构造完成后统一做一次列归一化并剔除范数接近零的列。MATLAB 里可以用col_norms sqrt(sum(D.^2, 1)); valid col_norms 1e-12; D D(:, valid); D D ./ col_norms(valid);这行代码放在任何字典生成函数末尾是最便宜也最有效的一道保险。5.2 稀疏度 K 设置过大降噪变成“追噪”现象训练集重构误差很低但放在测试集或者真实图像上重构结果明显包含颗粒状噪声图像反而变脏了。原因稀疏度 K 设得过高算法有足够自由度去拟合噪声。噪声在字典原子上的投影通常不为零只要允许更多原子进入支撑集它就一定会被当作有效信号学进去。解决不要拍脑袋给 K。先跑一个扫描K 取 1、2、4、6、8、10、12画出测试集重构误差随 K 变化的曲线找误差下降趋于平缓的拐点。图像块的稀疏度一般控制在每块 4 到 10 个原子一维信号可以根据物理含义估计非零频带数量再乘 1.2 倍作为上限。5.3 字典原子高度相关同一个基被反复换着选现象OMP 打印出来的支撑集索引不重复但残差在第三、第四轮之后就基本不降重构波形轮廓对细节全丢。原因过完备字典如果由 DCT 和小波混合拼接而成不同来源的原子之间可能存在强相关性。OMP 能排除同一个索引再次被选中但无法排除与已选原子高度相关的另一个原子。后选进来的原子不能提供新信息残差自然不降。解决字典生成完成后先算 Gram 矩阵G D*D把相关系数大于 0.9 的原子对找出来只保留其中一个。如果你希望保留完整字典来做对比实验那就改用 Lasso 类的 L1 方法它对相关原子的处理比贪婪算法更平滑。5.4 忽略负系数重构结果相位全乱现象重构信号的幅值范围正常但波形上下颠倒或者某一小段相位和原始信号完全相反。原因很多人看proj abs(D*r)习惯了以为稀疏系数也全是正的重构时给x_hat abs(x_hat)或者手动把负系数清零。实际上字典原子本身没有方向约束负系数代表用原子的反方向去组合信号这种表示是合法的甚至有物理意义。强制取绝对值等于把一半的表示能力扔掉了。解决OMP 代码里只对“选择阶段”取绝对值求解阶段必须保留符号。重构时直接用x_hat(selected) x_sub不要在外面套 abs。验证阶段抽查一下系数的正负分布如果全是正数大概率是你中间哪个环节手动取了绝对值。5.5 中文注释乱码代码没法协作现象用新版 MATLAB 打开别人发来的脚本中文注释全部变成乱码代码逻辑看不出来冷嘲一些注释后保存整个文件都乱了。原因MATLAB 不同版本对脚本文件的默认编码不同老版本常按 GBK 保存新版本默认 UTF-8。文件编码和编辑器解码方式不匹配就出现乱码。这是稀疏表示代码在小组协作里最容易踩的隐性坑因为算法本身没错问题全在文本层。解决统一保存为 UTF-8 编码在 MATLAB 主页的“预设项 - 编辑器/Debugger - 语言”里把字符编码设置为 UTF-8。如果你只需要自己跑建议代码注释尽量用英文避免编码问题干扰算法调试。已经乱码的文件用记事本或其他文本编辑器打开重新“另存为”并选择 UTF-8 编码再回到 MATLAB 打开才会恢复。6. 进阶验证与实用技巧把稀疏表示用在真实图像块上以及两个能复制到最后一步的习惯6.1 用真实图像块检验稀疏表示效果模拟信号验证做完了K-SVD 也训完了最后要看的是这套流程在真实图像上的表现。下面这段代码读取 cameraman 图像把图像切成 8×8 的重叠块取其中一块交给 OMP 重构。% 读取灰度图像并归一化到 0~1 范围 img double(imread(cameraman.tif)) / 255; if size(img, 3) 1 img rgb2gray(img); end % 把图像切成 8x8 的重叠块每个块展开成 64x1 列向量 % im2col 的输出是 64 x NN 为总块数 patches im2col(img, [8 8], sliding); % 取其中一个块作为测试样本 test_patch patches(:, 1000); % 用 64x128 的过完备 DCT 字典做稀疏编码 D makeDCTDict(64, 128); x_hat OMP(D, test_patch, 8, 1e-4); % 重建并计算 PSNR recon D * x_hat; psnr 10 * log10(1 / mean((test_patch - recon).^2)); fprintf(单图像块 PSNR %.2f dB非零系数数 %d\n, psnr, nnz(x_hat));如果你跑出来的 PSNR 偏低先看test_patch的均值是不是很大。图像块均值接近 0.5 时直流原子需要承担很大系数这会挤压其他原子的预算。一个实用的改进是编码前先减均值把均值单独存下来重构后加回去这样稀疏度能被更多用来刻画纹理和边缘。6.2 两个可以直接复制的参数习惯第一个习惯是为每个稀疏表示脚本都保留一个误差历史变量。无论是 OMP 的残差还是 K-SVD 的每轮重构误差都打成一个数组画出来。这个习惯能让你在 30 秒内定位“字典学习是否在收敛”而不是把求解器当黑匣子反复试参数。第二个习惯是固定随机种子。K-SVD 的字典初始化用到了randpermOMP 的测试信号也用到了随机数。如果你不固定种子今天跑出来的结果和明天完全不同参数调节就变成了纯玄学。在脚本开头写一行rng(2024)所有实验就都能复现这比任何调参技巧都优先。我自己的经验是所有稀疏表示参数问题最后都能拆成“字典结构和稀疏度之间的匹配问题”。每次出现异常结果先看支撑集索引是否合理再看残差曲线是否单调下降最后才是调阈值。把这三步固化成肌肉记忆之后稀疏表示从 MATLAB 代码到实际应用的距离会比你想象中短得多。希望帮到你。本文还有配套的精品资源点击获取
