简介:字典学习的MATLAB仿真资源包面向信号处理、图像处理与机器学习方向的初学者和研究人员,聚焦稀疏表示与字典更新方法,涵盖K-SVD、OMP/MOD、Elastic Net等经典算法。资源共9个文件,压缩包仅3.8MB,包括6个MATLAB脚本、1个编译好的mexw64加速函数、1个多模态样本数据mat文件及1份说明文档,代码模块划分清晰。实现围绕多模态字典学习与联合分类展开,覆盖字典投影、因子分解、随机梯度下降、联合ADMM求解等关键环节,并提供可直接运行的示例数据,方便观察稀疏编码残差与字典收敛过程。使用者可借助注释和文档梳理从字典初始化到更新迭代的完整流程,也可在图像压缩、特征提取、信号恢复等场景中替换数据验证算法。已有1614人学习下载,适合需要边读代码边复现实验的读者作为上手与二次开发的参考。
1. 字典学习整体设计与仿真思路
1.1 字典学习到底是解决什么问题的
前几年我一直在做图像去噪和压缩感知相关的课题,起初用的是固定变换基,比如DCT、小波,效果怎么说呢,常规场景下还行,但一旦遇到纹理复杂或者结构特殊的图像,稀疏表示系数不够稀疏,重建质量就明显掉下去。后来把目光转向了字典学习——说白了,就是别再拿着预先设计好的基函数去硬套数据,而是让算法自己从训练样本中“学”出一组更贴合数据结构的原子来。
字典学习的核心思想并不复杂:给定一个训练样本矩阵 Y,我们要找一个字典 D 和一个稀疏系数矩阵 X,使得 Y ≈ D * X,同时 X 的每一列都尽可能稀疏。这里的 D 就是我们要学习的“字典”,它的每一列叫做一个原子(atom)。形象点理解,固定基就像用普通话的拼音去拼各地方言,总有些不合适;字典学习则是先听一段方言录音,再归纳出适合这套方言的拼音规则,后面再用这个规则去转写效果自然会好很多。
这套方法在MATLAB里做仿真有一个天然优势:MATLAB的矩阵运算接口非常高效,无论是OMP(正交匹配追踪)这种逐列稀疏编码,还是SVD分解更新字典,都能直接用原生函数完成,不太需要自己费力去写底层线性代数实现。
1.2 为什么选择K-SVD作为主算法
字典学习的算法路线有好几条:有最朴素的MOD(Method of Optimal Directions),有在线字典学习(Online Dictionary Learning),还有后来工程上用得最多的K-SVD。我在仿真中最终选择K-SVD,主要看中三个点。
第一,K-SVD把“稀疏编码”和“字典更新”两个阶段拆得非常清楚。稀疏编码阶段先用OMP求每个样本在当前字典下的稀疏表示;字典更新阶段再逐个原子做SVD分解,找到能最大程度降低残差的更新方向。这两个阶段交替迭代,逻辑清楚,MATLAB代码写起来结构也明白,对于想入门字典学习的同学来说是个特别好的切入点。
第二,K-SVD对字典初始化的敏感度相对低。最早做MOD的时候,字典初始化的好坏会直接影响最终收敛结果,而K-SVD在更新字典时采用的是逐原子更新的策略,每个原子都会吸收和自己相关性最高的样本残差进行分析,整体鲁棒性好不少。我试过用DCT字典、随机高斯矩阵、甚至直接从训练样本中抽列向量做初始字典,K-SVD迭代到后面基本都能收敛到可用的字典。
第三,K-SVD的扩展性很好。后面很多改进算法,比如带标签信息的判别式字典学习D-KSVD、用于图像分类的LC-KSVD,都是在K-SVD骨架上做的修改。你先把基础版K-SVD仿真吃透了,后面读这些改进方法的论文、复现实验,会感觉顺畅很多,不用重新啃一遍地基。
1.3 仿真方案的总体框架
我的MATLAB仿真框架整体分四个模块,后面会逐一展开讲:
- 训练样本准备模块:从原始图像中按滑窗方式提取图像块,组合成训练矩阵。
- 稀疏编码模块:用OMP算法对每个训练样本求稀疏系数。
- 字典更新模块:用SVD分解逐原子更新字典。
- 应用验证模块:用学到的字典做图像去噪或者重建,量化PSNR/SSIM指标。
2. 字典学习核心细节与MATLAB关键函数解析
2.1 字典学习方程怎么用矩阵形式表达
在做MATLAB仿真之前,首先得把数学模型写在纸上。标准K-SVD的目标函数如下:
argmin_{D,X} ||Y - D*X||_F^2,s.t. 对每个 i,||x_i||_0 ≤ T0
这里的 ||·||_F 是Frobenius范数,就是矩阵所有元素平方和再开根号,用来度量重建误差;||x_i||_0 是零范数,指的是向量中非零元素个数;T0 就是稀疏度约束——每个样本最多用 T0 个原子来线性表示。
这个优化问题本身是一个NP难问题,所以K-SVD用的策略是启发式迭代优化:先固定字典D去求解稀疏系数X(这就是稀疏编码阶段),再固定稀疏系数X去更新字典D(这是字典更新阶段),交替进行。虽然不能保证全局最优,但在实际工程中,这种交替优化策略得到的解质量已经很高了。
2.2 OMP稀疏编码的MATLAB实现
在K-SVD里,稀疏编码阶段的工作可以理解为“拿着当前的字典,给每个样本找最精简的线性表示”。我这边用的是OMP(正交匹配追踪),它的数学原理是每一步都从字典中挑一个与当前残差相关性最强的原子,再通过最小二乘求解当前原子集合下的系数,更新残差,直到达到稀疏度约束或残差阈值。
MATLAB里可以直接用函数,也可以自己写一个精简版本。我是建议初学阶段一定自己写一遍OMP,因为后面你要给OMP加约束、改停止条件、甚至做批量加速,都得吃透这个算法内部流程。我贴一下自己项目里用的核心循环代码:
function X = omp(D, Y, T0) % D: 字典矩阵,大小为 n x K % Y: 训练样本矩阵,大小为 n x N % T0: 稀疏度约束,每个样本的非零系数个数上限 % X: 稀疏系数矩阵,大小为 K x N [~, N] = size(Y); [~, K] = size(D); X = zeros(K, N); for col = 1:N y = Y(:, col); residual = y; support = []; % 原子索引集合 coeffs = []; % 对应系数 for iter = 1:T0 % 计算当前残差与所有字典原子的内积绝对值 correlations = abs(D' * residual); % 排除已经选过的原子 correlations(support) = 0; [~, idx] = max(correlations); support = [support, idx]; % 用最小二乘更新当前支撑集下的系数 D_s = D(:, support); coeffs = D_s \ y; % 更新残差 residual = y - D_s * coeffs; % 如果残差足够小,提前终止 if norm(residual) < 1e-6 break; end end X(support, col) = coeffs; end end这个实现里有个细节需要注意:每次迭代选原子时只依据内积绝对值大小,不考虑原子之间的相关性。所以OMP的精髓在于,每选入一个新原子后都会用最小二乘重新计算一遍全部系数,把前面对已经选入原子系数的估计也修正过来。这就是它和MP(匹配追踪)最大的区别——OMP不需要反复调整已选原子的系数,收敛速度也快得多。
注意:如果训练样本量很大,或者字典原子数K很大,这个循环写法的效率会存在问题。批量加速的标准做法是预计算 D'*D 和 D'*Y,这样每次内积计算就不需要重新算整个矩阵乘法。后面我会在常见问题里展开讲。
2.3 字典更新阶段的SVD分解逻辑
字典更新是K-SVD和MOD最大的差异点。MOD在更新字典时是把所有原子一次性整体做最小二乘更新;K-SVD则是一个一个原子轮流更新,每次只处理当前原子,用SVD去找到最优的更新方向和对应系数。
具体流程是这样的:假设现在要更新字典D中的第j个原子 d_j,先记录到所有使用过这个原子的训练样本索引集合,记为 w_j = { i | X的第i行第j个元素不为0 }。然后计算去掉该原子贡献后的整体误差矩阵:
E_j = Y - sum_{k≠j} d_k * x_k^T
然后把E_j限制在w_j对应的列上,得到 E_j^R。对这个矩阵做SVD分解:
[U, S, V] = svd(E_j^R)
用U的第一列去更新原子 d_j,用S(1,1)乘以V的第一列去更新对应的稀疏系数行。为什么要这么做?因为SVD分解中,秩一近似是矩阵在Frobenius范数意义下的最优近似,所以用第一主成分去更新原子,能最大程度减小重建误差。
核心代码片段如下:
function D = ksvd_dict_update(D, X, Y) % X: 当前稀疏系数矩阵 % 逐原子更新 [~, K] = size(D); for j = 1:K % 找到使用过第j个原子的样本索引 wj = find(X(j, :)); if isempty(wj) continue; end % 限制X、D的行列范围 X_j = X(:, wj); X_j(j, :) = 0; % 计算误差矩阵 E_j = Y(:, wj) - D * X_j; % SVD分解 [U, S, V] = svd(E_j, 'econ'); % 更新原子和系数 D(:, j) = U(:, 1); X(j, wj) = S(1,1) * V(:, 1)'; end end这里有个地方容易踩坑:svd(E_j, 'econ') 的'econ'选项很关键。如果不加这个选项,当E_j的尺寸不是方形时,MATLAB会返回一个完整的m×m的U矩阵,浪费内存不说,后面取列还可能取错。用'econ'模式后,U只保留与E_j列数对应的前几列,正好满足我们的需求。
3. 完整仿真流程与核心代码实现
3.1 训练样本怎么准备
字典学习的训练样本通常是对原始图像做滑窗采样得到的图像块。我习惯用一个统一的脚本函数来做这个步骤,因为图像块尺寸、滑动步长、是否归一化这些参数会直接影响字典质量。
以一个经典的图像去噪任务为例:输入一张带噪图像,我用的训练块尺寸是8×8,滑动步长设为1个像素。这意味着每个位置都提取一个8×8的块,覆盖整个图像的所有局部区域。块尺寸选得太大,训练样本数会减少,字典学到的是偏全局的结构;选得太小,则很难捕捉到足够的纹理细节。实践中8×8是一个不错的平衡点,也是论文里最常见的配置。
样本提取后,需要把每个8×8的块拉成一个64维的列向量,拼成一个64×N的训练矩阵。所有样本还需要做均值移除——就是每个样本向量减去自身的均值。这一步非常关键,因为图像块中的直流分量(均值)会占据字典原子的一个重要维度,不去除的话,字典会花大量原子去编码亮度均值,而不是纹理结构。去噪时再把均值加回去就行。
function patches = im2col_patches(img, block_size, step) % img: 灰度图像,double类型,范围[0,1] % block_size: 块尺寸,比如8 % step: 滑动步长,设为1时重叠最大 [h, w] = size(img); patches = []; for i = 1:step:(h - block_size + 1) for j = 1:step:(w - block_size + 1) block = img(i:i+block_size-1, j:j+block_size-1); col = block(:); % 均值移除 col = col - mean(col); patches = [patches, col]; end end end这段代码为了清晰写成了双重循环,但实际大规模训练时建议用MATLAB自带的im2col函数,速度会快很多:
patches = im2col(img, [block_size, block_size], 'sliding'); patches = patches - mean(patches, 1);3.2 主循环:K-SVD迭代训练
有了前面的稀疏编码和字典更新模块,主循环就非常简单清晰了。下面给出一个完整的K-SVD训练脚本,字典大小设为256个原子,迭代次数设为30轮,稀疏度T0设为6:
% ==== 参数设置 ==== block_size = 8; % 图像块尺寸 dict_size = 256; % 字典原子数 T0 = 6; % 稀疏度 max_iter = 30; % 迭代轮数 % ==== 读取图像并提取训练样本 ==== img = double(imread('lena.png')) / 255; if size(img, 3) == 3 img = rgb2gray(img); end Y = im2col(img, [block_size, block_size], 'sliding'); Y = Y - mean(Y, 1); % ==== 字典初始化 ==== % 使用DCT过完备字典初始化,也可以用随机高斯矩阵 D = dct_dict(block_size * block_size, dict_size); % ==== K-SVD迭代 ==== X = zeros(dict_size, size(Y, 2)); for iter = 1:max_iter % 稀疏编码阶段 X = omp(D, Y, T0); % 字典更新阶段 D = ksvd_dict_update(D, X, Y); % 计算重建误差 recon_err = norm(Y - D * X, 'fro') / norm(Y, 'fro'); fprintf('Iter %d, relative error: %.4f\n', iter, recon_err); enddct_dict这个函数负责生成一个过完备DCT字典作为初始字典。从信号处理的角度来说,DCT基函数本身就能较好地表示平滑图像块,用它做初始值能显著减少前期迭代次数。我常用的实现:
function D = dct_dict(n, K) % n: 原子维度(比如64) % K: 字典原子数量 D = zeros(n, K); for k = 1:K freq = (k - 1) * pi / K; for i = 1:n D(i, k) = cos((i - 1) * freq); end end % 每列归一化 D = D ./ sqrt(sum(D.^2, 1)); end每轮迭代打印重建误差是个好习惯。以我的实验经验,前5轮误差下降会非常明显,从0.3左右快速降到0.1以下;到第15轮以后基本进入平台期,此时继续增加迭代轮数对最终效果的影响已经很小。如果你跑完30轮发现误差还在明显下降,说明训练样本量偏大或字典尺寸偏大,可以适当增加迭代轮数。
3.3 用训练好的字典做图像去噪
字典训练完成之后,应用阶段的核心思想是:对于一幅待去噪图像,先按同样的滑窗方式提取图像块,然后用OMP在当前字典下对每个块做稀疏编码,再用D * X重建出去噪后的图像块,最后加权平均还原成完整图像。这里重建误差大的成分绝大多数是噪声,而字典本身学到的是图像的结构信息,所以重建结果能够有效抑制噪声。
图像块重叠区域的处理需要额外注意。由于我采用的滑动步长为1,相邻图像块之间有大量重叠,每个像素会被多个重建块覆盖。最简单有效的合并方式是对每个像素位置累加所有覆盖它的重建值,再除以覆盖次数——这就是所谓的平均融合。
function denoised = denoise_image(img, D, T0) [h, w] = size(img); block_size = sqrt(size(D, 1)); block_size = round(block_size); Y = im2col(img, [block_size, block_size], 'sliding'); means = mean(Y, 1); Y = Y - means; % 稀疏编码 X = omp(D, Y, T0); % 重建 rec = D * X + means; % 融合回完整图像 % 使用accumarray思路,也可以直接用col2im的'average'方式 denoised = col2im(rec, [block_size, block_size], [h, w], 'sliding'); endcol2im函数在这里做了两件正确的事:把每个图像块向量重新排列回块矩阵,然后在重叠位置取平均。用'average'选项时,重叠区域是简单平均;用'sliding'选项时,MATLAB默认会在重叠区域做加权平均(实际上等价平均)。对于大部分场景,这个平均策略已经够用。
4. 常见问题、调参与效率优化实录
4.1 OMP求解慢,训练时间过长怎么办
我最早测试时,用256×256的Lena图、字典256原子、30轮迭代,直接跑原生循环OMP,单轮迭代就要约40秒,30轮下来20多分钟,这个速度做实验确实难受。
优化思路有三个方向。第一个是预处理:因为字典D在稀疏编码阶段是不变的,所以D'*D可以提前算好缓存,OMP每轮计算内积时直接查表而不是重新算整个矩阵乘法。第二个是减少不必要的计算:OMP内循环里,D_s = D(:, support) 这个子矩阵的乘法每次都在重建,可以预先计算残差与所有原子的内积,然后每次更新时只修正新加入原子带来的那部分影响。第三个是整体向量化:对多个样本同时做OMP,把内层循环改为矩阵操作。不过这个方法实现复杂度高,我建议新手先做前两个优化,通常就能提速3到5倍。
还有一种通用的备选方案是直接用批量稀疏编码函数,例如用l1-magic工具箱或SPAMS库。但如果学习K-SVD原理,还是建议先自己把OMP写出来,跑通了理解了,再考虑外部库。
4.2 字典训练不收敛或重建误差震荡
我遇到过一种情况:迭代到中段重建误差不仅不降反而小幅上升,这通常是两个原因造成的。
第一个原因是稀疏编码阶段的OMP没有严格按照稀疏度约束。比如某个样本在字典中确实找不到T0个足够好的原子,而OMP强制选了T0个,这些额外的原子反而引入了更大的重建误差。解决办法是在OMP里增加一个残差阈值判断,如果某一步残差已经小于设定阈值就提前停止,不要强行选满T0个原子。
第二个原因则是数值稳定性问题。当字典中有两个原子高度相关时,最小二乘求解D_s \ y可能会因为矩阵病态产生很大的系数值,进而导致字典更新时出现异常。解决办法是在每次原子更新后做一次归一化,并检查原子之间的互相关,遇到高度相关的原子组可以剔除其中一个后重新初始化。
4.3 字典可视化后观察到大量噪声型原子
每次训练完,我都会写一段可视化代码把字典原子打印出来看一眼:
function display_dict(D) n = sqrt(size(D, 1)); K = size(D, 2); margin = 2; grid_size = ceil(sqrt(K)); canvas = ones((n + margin) * grid_size, (n + margin) * grid_size); for k = 1:K atom = reshape(D(:, k), [n, n]); % 归一化到[0,1]便于显示 atom = (atom - min(atom(:))) / (max(atom(:)) - min(atom(:)) + 1e-8); r = floor((k - 1) / grid_size); c = mod(k - 1, grid_size); canvas(r*(n+margin)+1 : r*(n+margin)+n, ... c*(n+margin)+1 : c*(n+margin)+n) = atom; end imshow(canvas, []); end如果训练完的字典里出现大量这样的原子——频率分布杂乱无规律、和噪声图像块差不多的原子——基本可以断定是训练样本本身质量问题。最常见的错误是忘记做均值移除。当图像块均值没有被移除时,字典会分配相当一部分原子专门表示不同亮度级别的直流分量,这些原子从可视化角度看就是一块块均匀的色块,压缩了真正用于表示纹理结构的原子数量,字典表达能力会明显下降。
另一个容易导致噪声型原子出现的情况是初始字典中某些原子在训练过程中彻底丢失了更新机会,也就是所谓的“死亡原子”。在字典更新阶段,某些原子在整个训练集中没有被任何样本使用过,这时wj为空,原子就一直停留在初始值。后续如果初始值是随机的,这些死亡原子自然会显示成噪声状。处理方法是每次迭代结束后检查各原子的使用频率,对连续多轮没有参与编码的原子,用当前残差最大的训练样本块重新初始化。
4.4 参数调节的实战经验
参数调节方面,我把自己的经验和踩坑记录整理成了一张速查表,在做仿真实验之前不妨先对照一下:
| 参数 | 推荐范围 | 说明与心得 |
|---|---|---|
| 图像块尺寸 | 6×6到12×12 | 8×8最通用;块越大越偏全局结构,块越小越偏细节纹理 |
| 滑动步长 | 1到4 | 步长越小重叠越多,样本量越大,但训练耗时也越长;步长1效果最佳 |
| 字典原子数 | 128到512 | 256是均衡点;原子数太多容易过拟合,太少则表达能力不足 |
| 稀疏度T0 | 4到10 | 与原子数正相关;256原子时推荐6,128原子时推荐4到5 |
| 迭代轮数 | 20到50 | 用相对误差变化来判定收敛,不要盲目加大轮数 |
| 训练样本数 | 原子数的5到10倍 | 样本太少字典学不透,样本太多训练时间无法接受 |
单独强调一个容易被忽视的参数:训练样本量。我刚入门时认为样本越多越好,直到一次实验发现,用10000个样本训练的字典效果反而比用50000个样本训练的效果好。原因是图像块之间重叠度高,大量样本高度冗余,冗余样本并不会提供新的信息结构,反而拉偏了字典对重点结构的拟合。现在我的习惯是用im2col提取全部块后,随机均匀抽取其中的一部分作为训练集,通常取10000到30000个块就足够了。
5. 字典学习在MATLAB中的扩展应用
5.1 从图像去噪到图像超分辨率
K-SVD训练出的字典不仅能做去噪,还能直接扩展到超分辨率重建。超分辨率的基本思路是把字典学习拆成两个字典——低分辨率字典D_l和高分辨率字典D_h,两者共享同一个稀疏系数。训练时,把高分辨率图像块下采样成低分辨率图像块,然后分别对这两组块做联合字典训练,约束条件是它们在高维稀疏空间里的表示系数一致。
应用时,输入一幅低分辨率图像,提取图像块,用D_l求稀疏系数,再直接用D_h和这个系数重建出高分辨率图像块。我在MATLAB里做过一组对比实验,在3倍下采样场景下,比双三次插值的PSNR高出1.5到2.5dB,比传统的稀疏编码超分辨率方法(ScSR)效果略好,但计算时间增加了不少。这类实验很适合作为字典学习进阶方向,代码结构和K-SVD几乎一样,只是训练数据多了一路。
5.2 基于学习字典的图像分类
如果训练样本带类别标签,还可以做判别式字典学习。核心思想是在K-SVD的目标函数里加一个分类误差项:
argmin_{D,X,W} ||Y - DX||_F^2 + α * ||H - WX||_F^2,s.t. ||x_i||_0 ≤ T0
这里的H是样本的标签矩阵,W是一个线性分类器参数。这个优化仍然可以用类似K-SVD的交替迭代解法求解,区别在于字典更新后还要多加一步更新分类器W。在我做的一个手写数字识别实验中,用D-KSVD学到的字典在MNIST子集上能达到96%左右的准确率,比直接用原始像素训练SVM略高,而且字典具有很好的可解释性——每个原子大致对应某种笔画模式的组合。
6. 仿真实验过程中记得养成的几个习惯
代码写多了之后,我总结出几个对实验效率影响很大的习惯,这里一并分享给你。
第一,每次实验必须记录参数配置和随机种子。字典学习中用到了随机初始化(如果不用DCT初始化)、随机采样训练样本,这些随机因素会导致每次实验结果有细微差异。记录随机种子之后,任何实验结果都可以被精确复现,这在写论文或者做实验对比时太重要了。MATLAB里用rng(seed)这一行命令就能固定随机流。
第二,训练过程的误差曲线一定可视化出来。别只在命令行里打印数字,把每轮的相对重建误差存成数组,训练结束后plot出来。可视化比命令行打印更能反映训练健康状况:正常曲线是快速下降后平滑收敛,如果曲线有锯齿或者平台期极其漫长,说明参数设置有隐患,可以提前停止实验并调整参数,不用傻等全部训练完成。
第三,实验文件夹里固定保存一份字典可视化图片。我通常每次训练完成后立即执行display_dict函数,把最终的字典结构保存成png。从这些可视化图片中可以直观地比较不同参数下学到字典的差异,写论文时也需要这样的插图,所以保持这个习惯非常有益。
字典学习这套方法我在MATLAB里前前后后复现过很多次,从最开始的照搬论文公式、调不通参数,到后来能根据自己的数据特点灵活调整算法结构,中间踩过不少坑。希望这篇文章里写到的代码和调试心得,能帮你少走一些弯路。如果你在我描述的某个环节遇到了不一样的问题,欢迎在评论区把具体的报错信息、参数设置和图像尺寸贴出来,我看到了会尽量回复。最后再提醒一句:代码跑通只是第一步,真正理解每个矩阵维度变化背后的含义,才是做仿真实验最大的收获。
本文还有配套的精品资源,点击获取