简介字典学习是稀疏表示与压缩感知的关键技术旨在通过原子库和稀疏编码高效表达信号在图像去噪、分类、特征提取、压缩感知及通信等场景中应用广泛。该MATLAB仿真资源包以多模态字典学习为主线面向图像处理、机器学习与通信领域的研究者和工程师提供从理论到代码的完整实践入口。资源共9个文件包含6个.m脚本、1个mexw64加速模块、1个.mat多模态样例数据集和1个txt说明文档整体大小仅3.8MB。代码覆盖K-SVD、MOD、OMP、Elastic Net等经典字典学习算法并融合ADMM优化、在线任务驱动学习、多类决策融合、字典投影因子化等模块各类函数封装清晰配合样例数据可直接运行训练与测试流程可用于图像复原、特征提取、跨模态融合等实验也便于在readme指引下调整参数、对比不同算法性能。目前已有1614人学习使用适合具备一定MATLAB基础的中高级学习者作为算法研习、论文复现或工程预研的起点。1. 稀疏表示与字典学习先搞清楚它在解决什么问题两年前我第一次在MATLAB里跑字典学习仿真的时候卡了整整一周。不是K-SVD的公式看不懂而是没搞清楚稀疏编码和字典更新这两步到底怎么配合。后来我把整个流程手写了一遍才算真正啃下这块硬骨头。这篇分享就沿着我的实现思路来从稀疏表示的核心思想到K-SVD的MATLAB实现再到图像去噪仿真实验和几个逃不掉的坑一次说清楚。内容适合刚接触稀疏表示、想用MATLAB做仿真验证的同学也适合已经把工具包跑通但想弄懂底层原理的人。1.1 为什么偏偏是字典字典学习Dictionary Learning的核心目标是从一堆训练信号中自动学出一组基使得每个信号都能用这组基里尽可能少的原子atom线性组合出来。这个尽可能少就是稀疏性也是整个方法的灵魂。一个直观的类比假设你要用乐高积木拼各种造型。固定变换DCT、小波相当于只给一套固定规格的积木无论拼什么只能用这套而字典学习是根据你手上所有造型反过来设计一套最匹配的积木。信号处理里这套积木就是字典D每个造型的拼法就是稀疏系数a。数学上写成X D * A其中X是n×N的训练信号矩阵D是n×K的字典K通常大于n称为过完备字典A是K×N的系数矩阵要求每一列a_i的非零元素个数不超过T。写成优化问题就是min ||X - DA||_F²s.t. ||a_i||₀ ≤ T刚开始看到这个式子的同学都会问||·||₀是非凸的甚至求解它是NP难的这怎么优化K-SVD的聪明之处就在于没有直接解这个难题而是把问题拆成稀疏编码和字典更新两步交替迭代每一步都有成熟的高效解法。1.2 和固定字典相比学出来的好在哪我用一张表把差异列清楚方便理解对比项固定字典DCT/小波学习字典K-SVD原子来源预先定义的正交基函数从训练数据中自动学出数据匹配度一般假设信号结构固定高专门适配当前数据稀疏能力需要较多系数才能逼近同等逼近程度下系数更少适用范围对自然图像、音频有普适性训练集与测试集同分布时最佳计算成本低一次变换即可高需要离线训练一句话总结固定字典是以不变应万变学习字典是看菜下碟。对自然图像这类结构复杂的信号学习字典通常能用更少的原子达到更低的逼近误差这是它能在去噪、压缩感知、人脸识别等任务里刷出更好成绩的根本原因。1.3 做仿真前需要准备的MATLAB基础这篇文章默认你会基本的MATLAB脚本和函数编写了解矩阵切片、向量范数、SVD分解svd函数和最小二乘反斜杠运算。如果这些还不熟建议先跑一遍MATLAB自带demo或者找份速查手册。另外提醒一点下面代码里用了vecnorm它是R2017b之后才有的函数老版本MATLAB会直接报错可以用sqrt(sum(D.^2,1))替代。这个细节我是在帮同学调代码时发现的他的R2016a一运行就崩查了半天才定位到这一行。2. K-SVD算法拆解为什么交替优化能work2.1 两个子问题固定一个优化另一个K-SVD的全称是K-Singular Value Decomposition它的核心思想是把联合优化问题拆成两个子问题固定字典D求系数A也就是稀疏编码固定系数A更新字典D也就是字典更新。这就像装修房子你不可能同时决定墙怎么刷、家具怎么摆总是先固定一个再调整另一个反复迭代逼近。数学上这种交替优化不能保证收敛到全局最优但大量实验表明K-SVD在信号处理任务里都能落到一个很好的局部最优解。原因是每一步都让目标函数值严格下降OMP求出的系数是当前字典下的最优稀疏表示SVD更新又让字典在Frobenius范数意义下最优地逼近误差矩阵两步叠加重构误差单调递减实践中很少出现震荡。2.2 稀疏编码阶段OMP在做什么OMPOrthogonal Matching Pursuit正交匹配追踪是使用最广泛的稀疏编码算法。思路非常直接每次从字典里选一个与当前残差最匹配的原子再用最小二乘把已选原子上的系数整体重估一遍。注意是整体重估这是OMP与MP匹配追踪的关键区别——OMP保证每一步得到的系数都是在已选原子张成的子空间中的最优解而不会因为新选原子导致前面的系数失效。步骤拆解如下残差初始化为信号本身 r y计算所有原子与残差的内积取绝对值最大的那个原子加入支撑集在支撑集上解最小二乘 a_s D_s \ y更新残差 r y - D_s * a_s重复2~4直到选了T个原子或残差足够小。整个过程计算量主要集中在第3步的最小二乘K和T都不大时速度很快这也是K-SVD能在普通PC上跑起来的原因。2.3 字典更新阶段SVD在做什么字典更新是逐原子进行的。更新第k个原子时先找出所有用到了这个原子的训练样本计算去掉第k个原子后的误差矩阵E_k然后对E_k做SVD分解取最大奇异值对应的左奇异向量作为新原子同时用奇异值和右奇异向量更新对应的系数行。为什么用SVD因为SVD的秩1近似是Frobenius范数意义下的最优近似这保证了单步更新能把重构误差压到最小每更新一个原子误差都会下降一块。有个容易忽略的细节更新第k个原子时其他原子的系数完全不动。因为E_k只在用到第k个原子的样本列上计算SVD更新只影响第k列原子和第k行系数。这种局部更新设计是K-SVD最巧妙的地方也是它能比MODMethod of Optimal Directions收敛更快的原因。3. MATLAB从零实现手写OMP与K-SVD3.1 OMP函数关键是用反斜杠而不是inv先给完整的OMP实现注释写在关键行function a omp(D, y, T) % OMP 正交匹配追踪 % D: n x K 字典 % y: n x 1 目标信号 % T: 稀疏度上限 [~, K] size(D); r y; % 残差初始化为信号 idxSet []; % 支撑集原子索引 a zeros(K, 1); for iter 1:T proj abs(D * r); % 所有原子与残差的内积 [~, j] max(proj); % 找最匹配的原子 if ismember(j, idxSet) break; % 防止重复选同一个原子 end idxSet [idxSet, j]; Ds D(:, idxSet); aTemp Ds \ y; % 最小二乘必须用反斜杠 r y - Ds * aTemp; % 更新残差 if norm(r) 1e-8 break; end end a(idxSet) aTemp; end代码里几个故意这么写的点proj abs(D * r)本质上是在做投票把n维残差投影到K个原子方向上谁和残差越像内积绝对值越大ismember(j, idxSet)防止同一个原子被选两次虽然T循环次数有限但一旦重复选择会陷入无效迭代必须提前break。最值得注意的是第3步的最小二乘。很多初学者习惯写inv(Ds*Ds)*Ds*y这是给自己挖坑。当两个原子高度相似时Ds*Ds接近奇异inv会严重放大数值误差而反斜杠运算符内部走QR分解或Cholesky分解数值稳定性好得多速度也更快。这个习惯在MATLAB里几乎所有涉及最小二乘的地方都适用。3.2 K-SVD主循环K-SVD主循环的完整代码如下function [D, A] ksvd(Y, K, T, numIter) % K-SVD 字典学习 % Y: n x N 训练信号每列一个样本 % K: 字典原子个数 % T: 稀疏度 % numIter: 迭代次数 [n, N] size(Y); % 随机选K个训练样本做初始字典逐列归一化 rng(42); initIdx randperm(N, K); D Y(:, initIdx); D D ./ vecnorm(D, 2, 1); % 老版本用 sqrt(sum(D.^2,1)) for iter 1:numIter % 稀疏编码阶段 A zeros(K, N); for i 1:N A(:, i) omp(D, Y(:, i), T); end % 字典更新阶段 for k 1:K useIdx find(abs(A(k, :)) 1e-6); if isempty(useIdx) continue; % 这个原子没人用稍后处理 end % 计算移除第k个原子后的误差矩阵 E Y(:, useIdx) - D * A(:, useIdx); Ek E D(:, k) * A(k, useIdx); % SVD取主分量作为新原子并更新系数行 [U, S, V] svd(Ek, econ); D(:, k) U(:, 1); A(k, useIdx) S(1, 1) * V(:, 1); end err sum(sum((Y - D * A).^2)) / N; fprintf(iter %2d, 平均重构误差: %.6f\n, iter, err); end end初始化时我从训练样本里随机抽K列作为初始字典然后逐列归一化。这里有个细节如果初始字典不归一化第一轮OMP选原子时内积大小会被原子能量干扰导致选出来的原子偏向能量大的方向整个训练过程都会受影响。归一化后每个原子等价选核完全由方向决定。svd(Ek, econ)里的econ参数也很关键经济型分解在矩阵不是方阵时能省掉大量计算。如果写svd(Ek)默认会返回完整的U、S、V矩阵维度一大内存和时间都翻倍纯属浪费。3.3 图像块提取与重构平均字典学习在图像上一般不直接对整个图像做而是切成小块patch常用的块大小是8×8或16×16每个块拉成64维或256维的列向量。我没有用im2col因为Image Processing Toolbox不是人人都有自己写个循环函数更通用function patches extractPatches(img, patchSize, step) % 提取图像块每个块按列拉直 % img: H x W 灰度图 % patchSize: 块边长正方形 % step: 滑动步长 [H, W] size(img); patches []; for i 1:step:H-patchSize1 for j 1:step:W-patchSize1 p img(i:ipatchSize-1, j:jpatchSize-1); patches [patches, p(:)]; end end end重构时要把稀疏系数还原成图像块再按重叠位置平均回去。这里必须维护一个权重矩阵每个像素被多少个块覆盖就除以多少否则重叠区域的像素会过亮function recon reconstructImage(D, A, imgSize, patchSize, step) [H, W] imgSize; recon zeros(H, W); weight zeros(H, W); [~, N] size(A); idx 1; for i 1:step:H-patchSize1 for j 1:step:W-patchSize1 patch reshape(D * A(:, idx), patchSize, patchSize); recon(i:ipatchSize-1, j:jpatchSize-1) ... recon(i:ipatchSize-1, j:jpatchSize-1) patch; weight(i:ipatchSize-1, j:jpatchSize-1) ... weight(i:ipatchSize-1, j:jpatchSize-1) 1; idx idx 1; end end recon recon ./ weight; end步长step1时每个像素被64个块覆盖权重矩阵就是均值滤波器step3时图像边缘会有覆盖不均匀的情况除以weight能自动处理。这个重叠块平均技巧在图像去噪、超分辨率里是标配操作。4. 图像去噪仿真从取块到字典可视化4.1 完整实验脚本我拿经典的Lena灰度图256×256像素值归一化到0~1做实验。流程是给干净图叠加标准差为25/255的高斯白噪声用含噪图切成8×8图像块、步长取1跑K-SVD学字典再用OMP对每个块做稀疏重构最后平均回去imgClean im2double(imread(lena_gray.png)); sigma 25 / 255; rng(0); imgNoisy imgClean sigma * randn(size(imgClean)); patchSize 8; step 1; patches extractPatches(imgNoisy, patchSize, step); K 256; T 8; numIter p a hrefhttps://download.csdn.net/download/weiba_lu/10160896 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p 稀疏表示与字典学习：K-SVD算法原理及MATLAB图像去噪实现
简介字典学习是稀疏表示与压缩感知的关键技术旨在通过原子库和稀疏编码高效表达信号在图像去噪、分类、特征提取、压缩感知及通信等场景中应用广泛。该MATLAB仿真资源包以多模态字典学习为主线面向图像处理、机器学习与通信领域的研究者和工程师提供从理论到代码的完整实践入口。资源共9个文件包含6个.m脚本、1个mexw64加速模块、1个.mat多模态样例数据集和1个txt说明文档整体大小仅3.8MB。代码覆盖K-SVD、MOD、OMP、Elastic Net等经典字典学习算法并融合ADMM优化、在线任务驱动学习、多类决策融合、字典投影因子化等模块各类函数封装清晰配合样例数据可直接运行训练与测试流程可用于图像复原、特征提取、跨模态融合等实验也便于在readme指引下调整参数、对比不同算法性能。目前已有1614人学习使用适合具备一定MATLAB基础的中高级学习者作为算法研习、论文复现或工程预研的起点。1. 稀疏表示与字典学习先搞清楚它在解决什么问题两年前我第一次在MATLAB里跑字典学习仿真的时候卡了整整一周。不是K-SVD的公式看不懂而是没搞清楚稀疏编码和字典更新这两步到底怎么配合。后来我把整个流程手写了一遍才算真正啃下这块硬骨头。这篇分享就沿着我的实现思路来从稀疏表示的核心思想到K-SVD的MATLAB实现再到图像去噪仿真实验和几个逃不掉的坑一次说清楚。内容适合刚接触稀疏表示、想用MATLAB做仿真验证的同学也适合已经把工具包跑通但想弄懂底层原理的人。1.1 为什么偏偏是字典字典学习Dictionary Learning的核心目标是从一堆训练信号中自动学出一组基使得每个信号都能用这组基里尽可能少的原子atom线性组合出来。这个尽可能少就是稀疏性也是整个方法的灵魂。一个直观的类比假设你要用乐高积木拼各种造型。固定变换DCT、小波相当于只给一套固定规格的积木无论拼什么只能用这套而字典学习是根据你手上所有造型反过来设计一套最匹配的积木。信号处理里这套积木就是字典D每个造型的拼法就是稀疏系数a。数学上写成X D * A其中X是n×N的训练信号矩阵D是n×K的字典K通常大于n称为过完备字典A是K×N的系数矩阵要求每一列a_i的非零元素个数不超过T。写成优化问题就是min ||X - DA||_F²s.t. ||a_i||₀ ≤ T刚开始看到这个式子的同学都会问||·||₀是非凸的甚至求解它是NP难的这怎么优化K-SVD的聪明之处就在于没有直接解这个难题而是把问题拆成稀疏编码和字典更新两步交替迭代每一步都有成熟的高效解法。1.2 和固定字典相比学出来的好在哪我用一张表把差异列清楚方便理解对比项固定字典DCT/小波学习字典K-SVD原子来源预先定义的正交基函数从训练数据中自动学出数据匹配度一般假设信号结构固定高专门适配当前数据稀疏能力需要较多系数才能逼近同等逼近程度下系数更少适用范围对自然图像、音频有普适性训练集与测试集同分布时最佳计算成本低一次变换即可高需要离线训练一句话总结固定字典是以不变应万变学习字典是看菜下碟。对自然图像这类结构复杂的信号学习字典通常能用更少的原子达到更低的逼近误差这是它能在去噪、压缩感知、人脸识别等任务里刷出更好成绩的根本原因。1.3 做仿真前需要准备的MATLAB基础这篇文章默认你会基本的MATLAB脚本和函数编写了解矩阵切片、向量范数、SVD分解svd函数和最小二乘反斜杠运算。如果这些还不熟建议先跑一遍MATLAB自带demo或者找份速查手册。另外提醒一点下面代码里用了vecnorm它是R2017b之后才有的函数老版本MATLAB会直接报错可以用sqrt(sum(D.^2,1))替代。这个细节我是在帮同学调代码时发现的他的R2016a一运行就崩查了半天才定位到这一行。2. K-SVD算法拆解为什么交替优化能work2.1 两个子问题固定一个优化另一个K-SVD的全称是K-Singular Value Decomposition它的核心思想是把联合优化问题拆成两个子问题固定字典D求系数A也就是稀疏编码固定系数A更新字典D也就是字典更新。这就像装修房子你不可能同时决定墙怎么刷、家具怎么摆总是先固定一个再调整另一个反复迭代逼近。数学上这种交替优化不能保证收敛到全局最优但大量实验表明K-SVD在信号处理任务里都能落到一个很好的局部最优解。原因是每一步都让目标函数值严格下降OMP求出的系数是当前字典下的最优稀疏表示SVD更新又让字典在Frobenius范数意义下最优地逼近误差矩阵两步叠加重构误差单调递减实践中很少出现震荡。2.2 稀疏编码阶段OMP在做什么OMPOrthogonal Matching Pursuit正交匹配追踪是使用最广泛的稀疏编码算法。思路非常直接每次从字典里选一个与当前残差最匹配的原子再用最小二乘把已选原子上的系数整体重估一遍。注意是整体重估这是OMP与MP匹配追踪的关键区别——OMP保证每一步得到的系数都是在已选原子张成的子空间中的最优解而不会因为新选原子导致前面的系数失效。步骤拆解如下残差初始化为信号本身 r y计算所有原子与残差的内积取绝对值最大的那个原子加入支撑集在支撑集上解最小二乘 a_s D_s \ y更新残差 r y - D_s * a_s重复2~4直到选了T个原子或残差足够小。整个过程计算量主要集中在第3步的最小二乘K和T都不大时速度很快这也是K-SVD能在普通PC上跑起来的原因。2.3 字典更新阶段SVD在做什么字典更新是逐原子进行的。更新第k个原子时先找出所有用到了这个原子的训练样本计算去掉第k个原子后的误差矩阵E_k然后对E_k做SVD分解取最大奇异值对应的左奇异向量作为新原子同时用奇异值和右奇异向量更新对应的系数行。为什么用SVD因为SVD的秩1近似是Frobenius范数意义下的最优近似这保证了单步更新能把重构误差压到最小每更新一个原子误差都会下降一块。有个容易忽略的细节更新第k个原子时其他原子的系数完全不动。因为E_k只在用到第k个原子的样本列上计算SVD更新只影响第k列原子和第k行系数。这种局部更新设计是K-SVD最巧妙的地方也是它能比MODMethod of Optimal Directions收敛更快的原因。3. MATLAB从零实现手写OMP与K-SVD3.1 OMP函数关键是用反斜杠而不是inv先给完整的OMP实现注释写在关键行function a omp(D, y, T) % OMP 正交匹配追踪 % D: n x K 字典 % y: n x 1 目标信号 % T: 稀疏度上限 [~, K] size(D); r y; % 残差初始化为信号 idxSet []; % 支撑集原子索引 a zeros(K, 1); for iter 1:T proj abs(D * r); % 所有原子与残差的内积 [~, j] max(proj); % 找最匹配的原子 if ismember(j, idxSet) break; % 防止重复选同一个原子 end idxSet [idxSet, j]; Ds D(:, idxSet); aTemp Ds \ y; % 最小二乘必须用反斜杠 r y - Ds * aTemp; % 更新残差 if norm(r) 1e-8 break; end end a(idxSet) aTemp; end代码里几个故意这么写的点proj abs(D * r)本质上是在做投票把n维残差投影到K个原子方向上谁和残差越像内积绝对值越大ismember(j, idxSet)防止同一个原子被选两次虽然T循环次数有限但一旦重复选择会陷入无效迭代必须提前break。最值得注意的是第3步的最小二乘。很多初学者习惯写inv(Ds*Ds)*Ds*y这是给自己挖坑。当两个原子高度相似时Ds*Ds接近奇异inv会严重放大数值误差而反斜杠运算符内部走QR分解或Cholesky分解数值稳定性好得多速度也更快。这个习惯在MATLAB里几乎所有涉及最小二乘的地方都适用。3.2 K-SVD主循环K-SVD主循环的完整代码如下function [D, A] ksvd(Y, K, T, numIter) % K-SVD 字典学习 % Y: n x N 训练信号每列一个样本 % K: 字典原子个数 % T: 稀疏度 % numIter: 迭代次数 [n, N] size(Y); % 随机选K个训练样本做初始字典逐列归一化 rng(42); initIdx randperm(N, K); D Y(:, initIdx); D D ./ vecnorm(D, 2, 1); % 老版本用 sqrt(sum(D.^2,1)) for iter 1:numIter % 稀疏编码阶段 A zeros(K, N); for i 1:N A(:, i) omp(D, Y(:, i), T); end % 字典更新阶段 for k 1:K useIdx find(abs(A(k, :)) 1e-6); if isempty(useIdx) continue; % 这个原子没人用稍后处理 end % 计算移除第k个原子后的误差矩阵 E Y(:, useIdx) - D * A(:, useIdx); Ek E D(:, k) * A(k, useIdx); % SVD取主分量作为新原子并更新系数行 [U, S, V] svd(Ek, econ); D(:, k) U(:, 1); A(k, useIdx) S(1, 1) * V(:, 1); end err sum(sum((Y - D * A).^2)) / N; fprintf(iter %2d, 平均重构误差: %.6f\n, iter, err); end end初始化时我从训练样本里随机抽K列作为初始字典然后逐列归一化。这里有个细节如果初始字典不归一化第一轮OMP选原子时内积大小会被原子能量干扰导致选出来的原子偏向能量大的方向整个训练过程都会受影响。归一化后每个原子等价选核完全由方向决定。svd(Ek, econ)里的econ参数也很关键经济型分解在矩阵不是方阵时能省掉大量计算。如果写svd(Ek)默认会返回完整的U、S、V矩阵维度一大内存和时间都翻倍纯属浪费。3.3 图像块提取与重构平均字典学习在图像上一般不直接对整个图像做而是切成小块patch常用的块大小是8×8或16×16每个块拉成64维或256维的列向量。我没有用im2col因为Image Processing Toolbox不是人人都有自己写个循环函数更通用function patches extractPatches(img, patchSize, step) % 提取图像块每个块按列拉直 % img: H x W 灰度图 % patchSize: 块边长正方形 % step: 滑动步长 [H, W] size(img); patches []; for i 1:step:H-patchSize1 for j 1:step:W-patchSize1 p img(i:ipatchSize-1, j:jpatchSize-1); patches [patches, p(:)]; end end end重构时要把稀疏系数还原成图像块再按重叠位置平均回去。这里必须维护一个权重矩阵每个像素被多少个块覆盖就除以多少否则重叠区域的像素会过亮function recon reconstructImage(D, A, imgSize, patchSize, step) [H, W] imgSize; recon zeros(H, W); weight zeros(H, W); [~, N] size(A); idx 1; for i 1:step:H-patchSize1 for j 1:step:W-patchSize1 patch reshape(D * A(:, idx), patchSize, patchSize); recon(i:ipatchSize-1, j:jpatchSize-1) ... recon(i:ipatchSize-1, j:jpatchSize-1) patch; weight(i:ipatchSize-1, j:jpatchSize-1) ... weight(i:ipatchSize-1, j:jpatchSize-1) 1; idx idx 1; end end recon recon ./ weight; end步长step1时每个像素被64个块覆盖权重矩阵就是均值滤波器step3时图像边缘会有覆盖不均匀的情况除以weight能自动处理。这个重叠块平均技巧在图像去噪、超分辨率里是标配操作。4. 图像去噪仿真从取块到字典可视化4.1 完整实验脚本我拿经典的Lena灰度图256×256像素值归一化到0~1做实验。流程是给干净图叠加标准差为25/255的高斯白噪声用含噪图切成8×8图像块、步长取1跑K-SVD学字典再用OMP对每个块做稀疏重构最后平均回去imgClean im2double(imread(lena_gray.png)); sigma 25 / 255; rng(0); imgNoisy imgClean sigma * randn(size(imgClean)); patchSize 8; step 1; patches extractPatches(imgNoisy, patchSize, step); K 256; T 8; numIter p a hrefhttps://download.csdn.net/download/weiba_lu/10160896 stylecolor:#ec7500;font-size:14px; 本文还有配套的精品资源点击获取 /a img altmenu-r.4af5f7ec.gif srchttps://csdnimg.cn/release/wenkucmsfe/public/img/menu-r.4af5f7ec.gif stylewidth:16px;margin-left:4px;vertical-align:text-bottom;cursor:text; /p