简介稀疏表示是现代信号处理与机器学习中的核心概念这份MATLAB配套资料面向希望理解稀疏编码、字典学习及优化求解的学习者与研究者帮助其用少量非零原子表达信号并完成实验验证。包内共6个文件均为m脚本压缩包约4KB虽体积精简却包含主流程演示、数据读取与精度计算等环节的代码实现便于直接阅读、运行和修改。内容覆盖LASSO、BPDN等稀疏求解方法以及基于K-SVD的字典学习并延伸至图像去噪、压缩感知MRI重建等应用场景读者可据此梳理从建模、优化到字典训练的完整脉络观察参数设置对稀疏解与重建效果的影响。已有6567人学习下载适合作为课程学习、论文复现和算法对比的入门参考也为迁移到音频与生物医学信号处理等任务提供可操作的起点。1. 稀疏表示不是玄学一张图里只有少数原子真正起作用从一张灰度图上取 8×8 小块拉成 64 维列向量你会发现它几乎没法用单个基函数说清楚但用字典里三五个原子的线性组合就能逼近到肉眼难辨。这就是稀疏表示最朴素的样子信号 y 写成字典 D 中少数几列的加权和系数向量 x 里绝大多数元素是零。它解决的是一类很具体的问题——候选原子数远多于观测维度、解有无穷多个时怎么挑出那个「用最少原子说清楚」的解。做压缩感知重建、图像去噪、人脸识别、机械故障特征提取的人会反复撞上它新手常把它当成又一个正则化技巧熟手知道真正的分水岭在字典构造和停止条件上。下面按能跑通、能调参、能排错的顺序讲。2. 从 L0 到 L1把组合爆炸换成凸优化的那一步2.1 稀疏性的度量与 L0 为什么算不动稀疏表示的标准写法是 y Dx ε其中 y 是 M 维观测D 是 M×N 的字典且通常 N M过完备x 是 N 维系数。所谓稀疏指 x 中非零元素的个数远小于 N用 L0 范数记作 ||x||_0 ≤ K。约束形式有两种一种是 min ||x||_0 s.t. y Dx另一种是 min ||y - Dx||_2 s.t. ||x||_0 ≤ K。工程上更常用后者因为噪声永远存在硬要求精确重建没有意义。L0 的问题不在于定义不清晰而在于求解代价。要从 N 个原子里挑 K 个候选组合数是 C(N,K)。取 N1024、K5这个数接近 10^13。就算每个组合只花一微秒做一次最小二乘验证也要跑上好几个月。所以 L0 精确解在小规模、低稀疏度时还能用穷举硬扛规模一上去就必须放弃。我见过有人对着一张 256×256 图的几千个块直接调 L0 穷举跑了整晚没出结果这不是算法写错是问题本身不可解。理解这一点很关键后面所有算法——OMP、基追踪、ISTA、K-SVD——都是在绕开这次组合爆炸用不同的近似换可接受的时间。提示判断要不要用稀疏表示先看你的问题有没有「解不唯一 稀疏先验合理」这两个特征。如果 N ≤ M 且矩阵满秩直接最小二乘就行硬上稀疏模型纯属给自己加负担。2.2 L1 松弛的几何解释与 λ 的实际含义把 L0 换成 L1得到的就是最常用的两个形式。基追踪BP写 min ||x||_1 s.t. Dx y带噪版本 Lasso 写 min (1/2)||y - Dx||_2^2 λ||x||_1。L1 是凸的可以转成线性规划或二次规划求解这才是整件事能落地的根本原因。为什么 L1 能促稀疏几何上有个很直观的说法。L2 球的等高面是圆球和约束超平面相切时切点一般落在某个一般位置上分量都不为零。L1 球是菱形二维或多面体高维顶点恰好落在坐标轴上超平面扫过去时绝大多数情况下先碰到某个顶点或棱于是解里出现大量零。维度越高顶点越多碰上的概率越大。这不是严格证明但足以解释为什么实践中 L1 的解稀疏性这么好。λ 的实际含义要讲清楚否则参数没法调。当 λ ≥ ||D^T y||_∞ 时x 的最优解是全零向量对应的就是软阈值公式里 λ 的上界。也就是说λ 从 0 往上扫先经过一段解还不算稀疏的区域然后稀疏度快速上升最后全部归零。我一般会先算出 λ_max ||D^T y||_∞再在 [1e-4, 1e-2] × λ_max 这个区间做对数扫描用交叉验证或者固定的重建误差阈值挑一个。λ 太大解太稀疏重建残差大图像上表现为过度平滑、细节全丢λ 太小解不稀疏等于白做正则噪声全被拟进去。这个区间扫描比拍脑袋设一个 0.1 靠谱得多。2.3 三类求解器的选型贪心、凸优化、迭代阈值同样一个 L1 或 L0 模型求解路线差别很大选错了不是精度差一点而是慢几个数量级。路线代表算法优点明显短板适用规模贪心OMP、CoSaMP、子空间追踪实现简单稀疏度 K 直接可控原子相关性高时选错就回不来M、N 中等K 小凸优化基追踪、内点法、ADMM解质量稳定有理论保证单次求解慢参数敏感中小规模离线场景迭代阈值ISTA、FISTA每次迭代只做矩阵乘和软阈值内存友好收敛慢需调步长大规模、字典结构化选型的判断标准我一般这样定字典列之间相关性低、K 在 10 以内、需要实时性选 OMP要一个可信的基准解、做论文对比实验选凸优化字典是固定结构DCT、小波、规模大、能接受几十上百次迭代选 FISTA。CoSaMP 是 OMP 的加强版每轮选多个原子再剪枝回 K 个重建质量比 OMP 稳代价是要知道稀疏度 K 的准确估计。还有一条常被忽略的字典的互相干性 μ(D) max_{i≠j} |d_i^T d_j| 决定了所有贪心算法和一部分凸优化算法的可靠边界。μ 接近 1 时OMP 会反复选到同一批高相关原子结果全错。这时要么换字典要么改用带相关性惩罚的方法。经验值是 μ 控制在 0.5 以下比较安心DCT 字典的 μ 可以接近 0随机高斯字典在过完备度不太高时也还可以。3. MATLAB 跑通 OMP字典、残差与停止条件3.1 字典矩阵怎么建为什么必须列归一化字典决定了稀疏表示的上限选错字典后面调什么都白搭。常见四种构造方式解析基DCT、小波、Gabor、随机矩阵高斯、伯努利、学习字典K-SVD、在线字典学习、以及这几类的拼接。做图像块处理DCT 基是最省事的起点正交、互相干性低、有快速变换。% 方式一8x8 图像块用的 DCT 字典64 个原子 N 64; D_dct dctmtx(N); % dctmtx 返回的每一行是一个基转置成列 % 方式二过完备 DCT 单位阵拼接保证能精确表示任意块 D_over [D_dct, eye(N)]; % 64 x 128 % 方式三随机高斯字典调试算法用 M 64; N_atom 256; rng(0); % 固定随机种子结果可复现 D_rand randn(M, N_atom); % 无论哪种用之前必须做列归一化 function D normalize_dict(D) nrm sqrt(sum(D.^2, 1)); nrm(nrm eps) 1; % 防止除零 D bsxfun(rdivide, D, nrm); end D normalize_dict(D_over);列归一化这一步很多人嫌麻烦跳过后果非常直接OMP 每轮用内积 |d_i^T r| 挑原子列能量大的原子天然占便宜算法会反复选那几个没归一化的列最终解的能量分布严重失衡重建误差反而比不做稀疏还大。这不是玄学是内积被列尺度污染了。归一化之后每列单位能量内积才真正反映相关性。bsxfun在新版本 MATLAB 里可以用隐式扩展D ./ nrm替代老代码兼容性考虑我一般保留bsxfun。rng(0)这种固定种子在调试阶段一定加上否则每次跑出来结果不同你根本分不清是参数问题还是随机性。3.2 OMP 主循环的 MATLAB 实现OMP 的思路用一句话说完每轮挑和当前残差最相关的原子把已选原子上的最小二乘解出来更新残差重复到满足停止条件。核心就四步写成函数如下。function [x, supp, res_hist] omp_solve(D, y, K, tol) % OMP_SOLVE 正交匹配追踪 % 输入: % D M x N 字典列已归一化 % y M x 1 观测信号 % K 最大稀疏度最多选 K 个原子 % tol 相对残差停止阈值推荐 1e-3 ~ 1e-6 % 输出: % x N x 1 稀疏系数 % supp 被选中的原子下标 % res_hist 每轮残差范数用来画收敛曲线 [M, N] size(D); x zeros(N, 1); supp []; r y; % 初始残差就是观测本身 res_hist zeros(K, 1); x_sub []; for iter 1:K % 步骤 1: 在剩余原子上找和残差相关性最大的一个 proj D * r; proj(supp) 0; % 已选原子屏蔽掉避免重复选 [~, idx] max(abs(proj)); % 相关性接近 0说明残差里没有可解释的成分了 if abs(proj(idx)) 1e-12 break; end supp [supp, idx]; % 步骤 2: 在已选支撑集上做最小二乘重估 D_sub D(:, supp); x_sub D_sub \ y; % 比 pinv(D_sub*D_sub)*... 更稳 % 步骤 3: 更新残差 r y - D_sub * x_sub; res_hist(iter) norm(r) / max(norm(y), eps); % 步骤 4: 相对残差达标就停 if res_hist(iter) tol break; end end if ~isempty(x_sub) x(supp) x_sub; end res_hist res_hist(1:max(iter, 1)); end逐段说逻辑。第 1 步的proj(supp) 0是防止重复选原子的必要动作OMP 理论上不会重复选同一列因为正交化后残差与已选原子正交但数值误差下偶尔会屏蔽掉最省心。第 2 步为什么用\而不是pinv或(D_sub*D_sub)\D_sub*y\在方阵和非方阵上都走 QR 分解条件数不好时比正规方程稳得多正规方程会把条件数平方这是很多「跑出来一堆 NaN」的根源。第 3 步残差用相对值是因为信号绝对幅度差别很大绝对阈值没法通用。res_hist这个输出看起来不起眼但排错时非常有用。正常的 OMP 残差曲线单调下降如果出现平台或者上升基本可以确定是字典线性相关太高或归一化没做。调用示例用一个 8×8 图像块验证重建img double(imread(cameraman.tif)); blk img(101:164, 101:164); % 取一个块区域 y blk(:); % 64x1 D normalize_dict(dctmtx(64)); % 正交 DCT 字典 [x, supp, res] omp_solve(D, y, 10, 1e-3); fprintf(重建相对误差: %.4f%%, 使用原子数: %d\n, ... 100 * norm(y - D*x) / norm(y), numel(supp));用正交 DCT 字典时K10 就能把自然图像块重建到误差百分之几这已经能说明稀疏表示对图像块的有效性。换成过完备字典同样 K 下误差通常会再降一截因为原子更贴合局部结构。3.3 停止条件定了结果才可复现OMP 的停止条件比大多数人想的更重要三种写法各有适用场景。稀疏度 K 硬截断最简单跑满 K 轮就停。优点是输出维度固定做批量处理方便缺点是在噪声大的时候会强行用满 K 个原子去拟合噪声出现过度拟合。图像去噪场景下这个坑特别明显噪点被原子忠实表达去噪变成了保噪。相对残差阈值norm(r)/norm(y) tol更符合物理意义噪声水平高就把 tol 设大一点让算法早停。缺点是不同块的收敛速度差别大有的 3 轮就停有的要 30 轮批量跑的时候时间不可控。我一般把 K 和 tol 同时设上谁先触发谁停形成双保险。支撑集不再变化就停属于补充判据。加在循环里一行if numel(supp) 1 supp(end) supp(end-1)之类没意义正确的是记录上一轮的支撑集比较两轮是否完全一致。这个判据能省掉不少无效迭代。实际参数怎么定做压缩感知重建tol 设 1e-4 到 1e-6K 取理论稀疏度的 1.2 倍做余量做图像去噪tol 直接按噪声标准差 σ 换算经验公式是 tol ≈ σ·sqrt(M)/norm(y)配 K 取 6 到 12做分类特征K 反而不是重点重点是字典本身要有判别性。4. K-SVD 训练字典从固定基到数据驱动4.1 交替优化的两步与初始化策略固定基的问题在于它不认识你的数据。自然图像块、机械振动信号、人脸纹理它们的最优原子形态完全不同靠 DCT 或小波去逼近总归隔了一层。字典学习就是让字典从数据里长出来。目标是 min ||Y - DX||_F^2 s.t. ||x_i||_0 ≤ T其中 Y 是 M×L 的训练样本矩阵每列一个样本X 是 N×L 的系数矩阵。这个目标对 (D, X) 联合非凸但对单独一个变量是凸的所以走交替优化固定 D 求 X稀疏编码固定 X 更新 D字典更新。K-SVD 和 MOD 的区别在更新那一步MOD 用整体最小二乘一次更新所有原子K-SVD 逐个原子更新每次只解一个秩一逼近问题收敛更快、结果更好。初始化我一般用两种。随机挑 L 个训练样本做初始原子简单但可能挑到冗余样本或者用 DCT 字典初始化然后跑几轮迭代收敛更稳尤其在小样本时不容易陷入差局部解。样本量少于原子数时随机初始化有明显风险建议先用 PCA 降维再用前若干个主成分做初始字典。训练样本的预处理常被忽略每个样本块去均值再归一化。不去均值的话字典第一层会浪费大量能量去表达直流分量不归一化的话字典会被高能量样本主导低能量样本几乎学不到东西。4.2 K-SVD 的 MATLAB 实现与 SVD 更新细节function [D, X, err_hist] ksvd_train(Y, n_atom, n_iter, T) % KSVD_TRAIN K-SVD 字典学习 % 输入: % Y M x L 训练样本每列一个已去均值、归一化 % n_atom 字典原子数 N一般取 2M ~ 4M % n_iter 外层交替迭代次数 % T 稀疏编码阶段每个样本的最大非零系数个数 % 输出: % D 训练好的 M x N 字典 % X N x L 稀疏系数矩阵 % err_hist 每轮重建误差 [M, L] size(Y); N n_atom; % --- 初始化随机挑 L 个样本做原子 --- rng(0); idx randperm(L, N); D normalize_dict(Y(:, idx)); X zeros(N, L); err_hist zeros(n_iter, 1); for it 1:n_iter % 阶段一稀疏编码 for n 1:L X(:, n) omp_solve(D, Y(:, n), T, 1e-6); end % 阶段二逐原子更新 for k 1:N % 只保留用到第 k 个原子的样本 use_idx find(X(k, :) ~ 0); if isempty(use_idx) continue; % 该原子无人使用跳过 end % 计算去掉第 k 个原子后的残差矩阵 Ek Y(:, use_idx) - D * X(:, use_idx); % 关键把第 k 个原子自己的贡献加回来否则 SVD 会把它当噪声消掉 Ek Ek D(:, k) * X(k, use_idx); % 对残差做秩一 SVD 逼近 [U, S, V] svd(Ek, econ); D(:, k) U(:, 1); X(k, use_idx) S(1, 1) * V(:, 1); end % 记录当前重建误差 err_hist(it) norm(Y - D * X, fro) / norm(Y, fro); end end代码里最容易写错的是Ek Ek D(:, k) * X(k, use_idx)这一行。K-SVD 的推导里去掉第 k 个原子后的残差要经过 SVD 找最优的一维方向但分子里必须把这个原子原本的贡献减掉再单独处理。漏掉这一步SVD 得到的方向会偏向剩余的成分字典更新方向完全偏了表现出来就是误差曲线振荡、训练不收敛。这是我自己踩过的坑调试了大半天才定位到这一行。svd(Ek, econ)用经济型分解Ek 是 M×|use_idx|完整 SVD 会算出一堆零空间浪费内存。更新系数时X(k, use_idx) S(1,1) * V(:,1)这里把奇异值也乘回去了这样重建能量守恒。另外一个小细节SVD 出来的 U(:,1) 符号不唯一可能整体翻转。对单个原子的符号翻转不影响重建但如果后续要做字典可视化和对比建议统一符号比如把 U(:,1) 里绝对值最大的元素强制为正。4.3 训练参数表与收敛判据参数含义典型取值影响n_atom字典原子数 N2M ~ 4M太大易过完备、μ 上升太小表达力不足n_iter外层迭代次数15 ~ 30超过 30 通常收益很小T稀疏度上限4 ~ 12太小欠拟合太大会把噪声学进字典L训练样本数至少 20N样本太少字典不稳定块大小 M单个样本维度648×8或 25616×16与要处理的信号尺度一致收敛判据我一般看两条err_hist相邻两轮的相对下降小于 1e-4或者跑到预设迭代次数。单看误差曲线的话健康的情况是前 10 轮快速下降之后进入缓慢下降的平台期如果出现上升或者锯齿优先查 4.2 里那行加回残差的代码其次查归一化。训练完的字典最好存下来。一次性训练 20 轮拿到的字典可以反复用于同类型数据避免每次推理都重训。存的时候用save(dict.mat, D)同时把训练时的块大小、T、预处理参数一起存进去换个人接手或者自己过几个月回来看没有这些参数信息基本等于废纸。5. 稀疏表示落地时的常见问题排查五个翻车现场5.1 残差曲线不降反升最后跳出 NaN现象OMP 跑到第 5 轮以后残差开始抖动再往后直接出 NaNx里全是Inf。原因两个可能叠加。一是字典列没归一化内积挑原子的标准被列能量污染选出来的支撑集内部高度相关二是最小二乘用了(D_sub*D_sub)\(D_sub*y)这种正规方程写法支撑集里两列接近共线时条件数爆炸平方之后直接溢出。解决字典用前必做列归一化最小二乘一律改用D_sub \ y内部走 QR 分解条件数只放大一次。如果换完还抖就检查字典的互相干性max(max(abs(D*D - eye(N))))超过 0.9 说明原子太像要么换字典要么在选原子时加相关性惩罚。5.2 稀疏度设大了去噪反而把噪点留下来现象图像去噪任务K 从 6 调到 12PSNR 先升后降视觉上出现一层细碎斑点。原因噪声在字典里总能找到几个原子能表达它K 越大可用的原子越多算法就会把这些噪声成分也精确拟合进去。稀疏表示的「稀疏」约束本身是抗噪的但前提是 K 小于真实信号所需的稀疏度。解决K 的选取要跟噪声水平绑定。经验做法是在干净数据上先统计每个块的真实稀疏度分布取 90% 分位数作为上限再根据噪声标准差适当下调。我一般用 K ceil(0.8 * K_clean) 起步配合相对残差阈值 tol ≈ 3σ让算法在拟合噪声前停下来。5.3 训练字典和测试样本预处理不一致现象字典在训练集上重建误差很低换到测试集上误差翻倍而且不同测试样本之间波动极大。原因训练时每个块做了去均值和归一化测试时只做了其中一步或者归一化用的范数不一样训练用整个块测试用块内非零部分。把测试系数还原成图像时又忘了把均值加回去结果整体偏暗或偏亮。解决预处理写成独立的函数训练和测试都调同一个。去均值时要记录每个块的均值向量重建时y_rec D*x mean_vec。如果对系数也做了归一化还要按范数还原。参数记录进.mat文件跟字典一起存别靠记忆。5.4 K-SVD 训练完字典原子几乎全是方波现象把训练出的 8×8 原子逐个画出来一堆长得像阶跃、方波没有预期的边缘纹理结构。原因训练样本没去均值。图像的直流分量占比很大字典前几轮会把大量原子分配给表达常量块剩下的原子又在剩余误差里互相争夺最终学出来的都是低频形状。T 设得太小也会加剧因为欠拟合时低频分量最容易被优先满足。解决样本零均值化后再训练并把均值作为独立通道单独处理。T 至少设到 6样本块数保证在原子数的 20 倍以上。训练中途把字典画出来看一眼原子有明显纹理结构边缘、角点、条纹才算健康这一步花两分钟能省几小时。5.5 重建结果整体偏移但误差数值看起来很小现象压缩感知重建出来的图整体色调偏灰残差范数不算大可视觉上就是不对劲。原因L1 求解时只求了系数忽略了观测里的直流偏置。或者 OMP 的停止条件只用了相对残差而相对残差对整体偏移不敏感——整体加个常数只让范数变一点点。解决重建后单独检查残差的均值mean(y - D*x)明显非零说明存在系统偏置。要么在观测端去均值、重建后加回要么在字典里显式加入全一向量作为原子让模型自己解释直流。6. 把稀疏表示用到图像去噪换一个惩罚项就能复用的套路6.1 重叠块 全局字典的最小实现图像去噪是稀疏表示最容易见效的场景。思路是把图切成重叠的小块每块单独做稀疏编码重建再按位置平均回去。重叠是为了消除块边界平均是为了压制单块解的不稳定性。function img_den denoise_sparse(img, sigma, D, T, step) % 重叠块稀疏表示去噪 % img: 含噪图像, sigma: 噪声标准差, D: M x N 字典 (M bs^2) % T: 每块稀疏度, step: 滑动步长 bs sqrt(size(D, 1)); % 块边长例如 8 [H, W] size(img); acc zeros(H, W); wgt zeros(H, W); for i 1:step:(H - bs 1) for j 1:step:(W - bs 1) blk img(i:ibs-1, j:jbs-1); blk blk - mean(blk(:)); % 去均值 nrm norm(blk(:)); if nrm 1e-8, continue; end y blk(:) / nrm; % 归一化与训练一致 tol 3 * sigma / (nrm eps); % 噪声相关的残差阈值 x omp_solve(D, y, T, tol); rec reshape(D * x, bs, bs) * nrm; acc(i:ibs-1, j:jbs-1) acc(i:ibs-1, j:jbs-1) rec; wgt(i:ibs-1, j:jbs-1) wgt(i:ibs-1, j:jbs-1) 1; end end img_den acc ./ max(wgt, eps); end两个参数最关键。step取 bs/2 是常见的折中重叠率 75%去块效应明显计算量是块数的 4 倍要更快就把 step 取到 bs代价是会看到明显块边界。tol跟着噪声走噪声大就早停这一步比调稀疏度更有效。6.2 换惩罚项等于换算法同一套框架里把 OMP 换成软阈值迭代就成了经典的 K-SVD 去噪流程本质是 L0 换 L1 的求解方式差异。软阈值的写法只有一行x sign(D*y) .* max(abs(D*y) - lambda, 0)当字典正交时它就是 L1 的精确解过完备时是近端梯度的一步近似跑十几轮也够用。我在资源受限的板子上做图像降噪时更倾向这条路矩阵乘加软阈值实现短、可向量化。参数上的经验8×8 块配 256 原子的学习字典σ ≤ 20 时 T 取 68 就够σ 更大时把 T 提到 10 但必须同步放大 tol否则噪点被忠实表达。图像去噪的 PSNR 提升通常集中在 σ 较小的区间噪声特别大时稀疏表示的边际收益会下降这时候混合 BM3D 那种非局部思路更划算。要说教训最有价值的一条是先检查归一化和预处理再去调稀疏度。我在一个项目上花了整整两天调 K 和 λ最后发现训练脚本里归一化写成了行归一化字典每一列的能量完全是乱的参数怎么调都没用。后来我养成习惯任何稀疏表示代码跑通第一件事就是打印字典列范数的最大值和最小值两者差超过 1e-6 就说明归一化有问题先修这个再说别的。希望帮到你。本文还有配套的精品资源点击获取