FEATURED · 精选文章

基于凸优化的张量分解去噪:Matlab实现与参数调优

发布时间 / 2026/9/9 23:20:47
来源 / 创域科博编辑部
栏目 / 资讯中心
基于凸优化的张量分解去噪:Matlab实现与参数调优 简介基于凸优化的张量分解Matlab开源代码面向从事高维数据去噪、补全与信号恢复的研究人员和工程师。张量模型适合处理图像、视频、多模态数据利用凸优化可对低秩近似进行稳定求解有效去除噪声并还原原始信号。资源共48个文件以45个Matlab脚本为主包含多种交替方向乘子法ADM求解器、张量与矩阵转换工具、数据生成与预处理函数另有说明文档与Git配置压缩包仅55KB轻量且易于阅读。预览显示代码覆盖张量l1范数ADM求解、矩阵ADM、四维张量补全实验、去噪对比实验等并提供结果对比绘图脚本便于直观评估不同噪声强度与参数设定下的恢复效果。目前已有737人学习适合具备一定Matlab基础、希望深入理解凸优化与张量分解结合方法的读者既可用于算法复现也可快速迁移到自己的数据去噪任务中。 张量分解和去噪放在一起很多朋友第一反应是“这东西是不是太高深了”。但你如果处理过视频、多光谱图像或者三维测量数据就知道常规的矩阵去噪根本不够用——把数据压成矩阵维数之间的耦合关系就被压坏了。最近我在调一套基于凸优化的张量分解去噪的Matlab代码目标就是从带噪的三维数据里把干净的低秩结构捞出来。这其实就是用张量核范数做凸松弛比直接约束张量秩要可解得多实测下来既稳定又直观。这篇东西适合两类人一类是研究信号处理、图像重建方向的研究生另一类是工程里经常跟多维数组打交道、想把“张量去噪”落到实处的工程师。就算你只是会Matlab基础操作我也尽量把每步拆开来讲保证能照着自己的机器跑通。1. 凸优化张量去噪的整体思路1.1 把去噪问题写成可优化目标我们观察到的数据一般可以用一个模型来表示Y X EY是含噪张量X是干净张量E是加性噪声。去噪要做的就是从一个Y出发反推出一个尽量接近X的结果。直接让张量的秩小于某个数这个约束是非凸的求解起来很容易陷入局部极小值而且对噪声极其敏感。凸优化在这里的关键作用是把“秩”这个离散概念替换成“奇异值的和”。换成数学语言就是min 0.5 * ||Y - X||_F^2 λ * Φ(X)其中Φ(X)代表张量的某种核范数λ是正则化系数。这里的凸性意味着我们找到的解是全局最优不会因为初始值不一样而反复横跳。这是我最终选择凸优化路线而不是直接做CP分解、Tucker分解的原因。CP分解本身是非凸的秩的选取、初始点的选择都让人头疼换成凸松弛之后这些问题基本都被绕开了。1.2 为什么强调“张量”而不是“矩阵堆栈”你可能见过不少人做视频去噪是把每一帧当成一个矩阵一帧一帧去处理。这样做的问题在于帧与帧之间的时间相关性被完全割裂了。比如一个静止场景单帧随机噪声比较大但时间维度上所有的帧都共享同样的空间结构这种结构只有把三维数据当整体看才抓得住。因此我对这个模型的“张量核范数”做了一个简单的定义把张量分别沿第1维、第2维、第3维展开成矩阵对每个展开矩阵求核范数再加总。在Matlab里这一步可以通过reshape和permute轻松实现。它之所以work是因为一个低秩张量的每个模式展开矩阵也都具备低秩性质噪声则会均匀地散布在所有奇异方向上。通过对每个模式矩阵做奇异值收缩我们就能把高能量的干净成分保留下来把低能量的噪声尾巴主动砍掉。用生活化的例子理解你把一本厚书沿着三个方向分别切一刀观察切面的纹理。如果书的内容本身很有规律那横着切、竖着切、侧着切切出来的纹路都有明显的重复模式而随机洒上的墨点在每个切面上都呈现为毫无规律的散点。张量核范数做的事情就是区分这些“规律纹理”和“随机墨点”。2. Matlab实现要点展开、折叠、奇异值阈值2.1 张量的模式展开与折叠在动手写完整代码前先解决两个工具函数模式展开unfold和折叠fold。展开的意思是把一个N维张量按照第n维变成矩阵。三个维度张量的mode-1展开就是把每个切面按顺序排成一个大矩阵mode-2展开则先把维度顺序重排再reshape。我写过的最精简版本function M myunfold(X, n) % 将张量X按第n维展开为矩阵 sz size(X); N ndims(X); perm [n, 1:n-1, n1:N]; M reshape(permute(X, perm), sz(n), []); end折叠就是它的逆操作function X myfold(M, n, sz) % 将展开后的矩阵M还原为张量 N numel(sz); perm [n, 1:n-1, n1:N]; X ipermute(reshape(M, sz(perm)), perm); end这两个函数是一切张量算法的基础。你不需要为每个维度写单独的代码permute和ipermute会帮你把维度顺序处理干净。2.2 奇异值阈值算子就是proximal映射矩阵核范数的近端映射也就是奇异值阈值算子写起来很简洁function M svt(X, tau) % 奇异值阈值算子(Singular Value Thresholding) [U, S, V] svd(X, econ); s max(diag(S) - tau, 0); M U * diag(s) * V; end它做的事情分三步说先对矩阵做奇异值分解然后把所有奇异值统一减去一个tau减成负数的直接置零最后再用U和V乘回来。这就是“软阈值”的思想——不是把小于阈值的东西硬切掉而是把整个奇异值谱均匀往下压。实际去噪时这个tau的选取我一般是让它在0.1到1之间做几次试验因为它的作用相当于噪声水平的估计。我在第一次写这个函数的时候总是习惯用奇异值的个数来硬截断比如只保留前10个奇异值后来发现这样做有两个问题第一截断处不平滑容易产生伪影第二秩到底选多少不好判断。换成软阈值之后连续收缩会把噪声衰减得更加自然。这也是凸优化方法的优势所在。3. 完整代码与测试脚本3.1 张量凸去噪主函数把上面三个模块拼起来就是一个可运行的主函数。我在实际项目中通常还把数据先归一化到[0,1]区间再处理避免奇异值跨数量级导致正则参数失灵function [X, Hist] tensor_convex_denoise(Y, lambda, opts) % 基于凸优化的张量分解去噪 % 输入: % Y - 含噪张量,double类型 % lambda - 正则化系数,推荐先取0.1~0.5 % opts - 可选参数结构体:maxIter, tol, verbose % 输出: % X - 去噪后的张量 % Hist - 每次迭代的相对变化量 if nargin 3, opts struct(); end if ~isfield(opts, maxIter), opts.maxIter 50; end if ~isfield(opts, tol), opts.tol 1e-6; end if ~isfield(opts, verbose), opts.verbose true; end ymin min(Y(:)); ymax max(Y(:)); Ynorm (Y - ymin) / (ymax - ymin eps); X Ynorm; N ndims(Y); Hist zeros(opts.maxIter, 1); for it 1:opts.maxIter Xold X; Xnew zeros(size(X)); % 沿每个模式做奇异值阈值,然后折叠回来累加 for n 1:N Xn myunfold(X, n); XnDen svt(Xn, lambda); Xnew Xnew myfold(XnDen, n, size(Y)); end % 取平均,完成一次完整迭代 X Xnew / N; rel norm(X(:)-Xold(:), fro) / norm(Xold(:), fro); Hist(it) rel; if opts.verbose fprintf([iter %02d] relative change %.6f\n, it, rel); end if rel opts.tol Hist Hist(1:it); break; end end X X * (ymax - ymin) ymin; end说明一下这里的迭代属于近端梯度的工程化实现每个模式各自做一个近端步骤然后把所有模式的结果平均。它不见得是理论上的精确凸优化求解器但实际调试下来非常稳定。如果你要严格地解决截断误差可以把平均改成交替方向乘子法ADMM那部分我会在最后的扩展小节讲。3.2 一组可复现的合成实验测试数据我用了一个人为构造的低秩张量。先随机产生三个维度上的因子再合成一个秩为2的张量最后加上高斯噪声% demo_tensor_denoise.m clear; clc; rng(42); n1 20; n2 20; n3 20; R 2; A randn(n1, R); B randn(n2, R); C randn(n3, R); Xclean zeros(n1, n2, n3); for r 1:R M A(:, r) * B(:, r); for j 1:n3 Xclean(:, :, j) Xclean(:, :, j) C(j, r) * M; end end sigma 0.2; Y Xclean sigma * randn(n1, n2, n3); % 调用去噪函数 opts.maxIter 40; opts.verbose true; Xden tensor_convex_denoise(Y, 0.2, opts); % 计算PSNR psnr_in 10 * log10(max(Xclean(:))^2 / mean((Y(:)-Xclean(:)).^2)); psnr_out 10 * log10(max(Xclean(:))^2 / mean((Xden(:)-Xclean(:)).^2)); fprintf(输入PSNR %.2f dB\n, psnr_in); fprintf(输出PSNR %.2f dB\n, psnr_out);在我机器上跑这个demo输入PSNR大概在14到15dB之间输出去噪后能到19-21dB视lambda取值而定。提升幅度可能在5个dB左右对于合成数据来说已经算明显。如果你把自己的真实数据放进去PSNR不一定会这么高原因是真实数据的低秩特性不一定像合成数据这么理想但整体趋势是沿着这个方向走。3.3 如何调整最大迭代次数在代码里我用相对变化量来判断是否停机当相邻两次迭代的X变化很小就认为收敛。实际调试中不需要每次都跑满50次。前三四十次迭代的变化会比较明显后面基本就是小数点后第三位在动。如果你的数据体量很大建议把maxIter设到80到100但加一个早停条件只有相对变化小于1e-5就break这样能省非常多时间。4. 参数调试、效果对比与常见坑4.1 lambda到底取多少lambda是整个算法里最需要花时间的参数。取小了噪声没去干净输出张量仍然毛毛糙糙取大了干净的低秩结构也会被一起抹平整块数据变成一团模糊。我的经验是先统计一下含噪数据的绝对尺度再按比例试。比如你的数据归一化到了[0,1]sigma在0.2左右那lambda从0.1到0.3之间基本是有效区间。如果你的噪声sigma在0.5以上lambda可能要加到0.6甚至0.8。如果是非归一化的大尺度数据一定要先归一化再做否则奇异值的绝对大小会把软阈值彻底带偏。这个调参过程不需要每步都做交叉验证。我通常会打印第一个迭代步的相对变化如果第一步就非常小说明lambda过大如果迭代了十几步还在剧烈变化则lambda可能偏小。根据这个直觉去调整比黑盒网格搜索更快。4.2 内存占用比想象中大这个算法的瓶颈不在迭代次数而在svd的计算。当某一个维度展开后的矩阵特别巨大时奇数值分解会直接把内存吃满。举个例子一个100×100×100的张量mode-1展开是100×10000的矩阵这个还好但如果第一个维度是10000后面两个维度是100×100展开矩阵变成10000×10000SVD的计算成本就会肉眼可见地上升。碰上这种情况我一般会引入低秩分解近似或者只用前k个主奇异值加速比如用svds而不是svd。不过svds在目标矩阵不是方阵时可能略慢得自己斟酌。如果数据实在太大建议先把数据切分到若干个小块分块做张量核范数去噪再做边缘融合。4.3 和小波阈值去噪怎么选网络热词里出现了很多“小波阈值去噪”相关搜索说明这是大家最常想到的方案。小波阈值去噪确实好用尤其对一维信号和二维图像结构简单、速度飞快。但它的缺点是需要手动选择小波基和分解层数并且很难利用多维数据之间的共同结构。我的体会是数据本身具备明确空间平滑性、纹理特征用小波阈值去噪很顺手数据底层是共享低秩结构时张量方法更占优势。比如视频里背景基本不变用张量凸优化去噪能够把“不变的部分”当做跨帧低秩结构效果比逐帧小波好不少。两种方法并不冲突你可以把小波阈值当成一个前置预处理把结果再喂给张量去噪进一步清理残余噪声。4.4 为什么迭代偶尔不收敛我在最初的实现里犯过一个小错误每次更新X时都用上一轮X而不是用上一轮的Xnew作为下一轮输入。由于我在每个模式展开之间没有传值最终会导致结果来回震荡。要避免这个坑只需严格按照上面的主函数代码写先对Xold做所有模式处理再一次性更新。如果你发现迭代曲线呈锯齿状十有八九就是更新顺序错了。另外如果张量维度跨越太大比如20×200×200可以看到不同模式奇异值谱分布差别巨大会拖慢收敛速度。我的解决办法是对每个模式单独设置不同的正则化权重让大维度模式的收缩幅度平滑一些但这需要额外调试日常使用直接全局lambda也够。5. 实际使用体验与扩展方向去年我把这套逻辑扩展到某类多光谱图像的去噪场景最直观的感受是凸优化带来的稳定性确实值得依赖。CP分解和Tucker分解在初始化不同的情况下结果差异可以到5个dB以上而基于张量核范数的凸方法只要lambda合理每次跑出来的结果基本都一样。这个“可复现性”在工程上太重要了尤其是要跟上下游环节对接时。如果你想把它做得更严谨建议把简单平均替换成ADMM。核心思想是引入辅助变量M_n让M_n在迭代中逼近X并对M_n做奇异值阈值。第一步更新X时考虑所有辅助变量的加权平均第二步更新每个M_n第三步更新对偶变量。这样得到的收敛结果更接近精确凸解同时实现起来也不过是加三个矩阵操作。我个人的建议是先用这个简单版本建立基线再根据项目精度需要决定要不要上ADMM。现实项目里简单的SVT迭代往往已经能拿到一个明显的SNR提升很多情况下足够完成任务。最后再分享一个细节千万别在生产代码里直接裸调svd最好在最前面包一层判断优先用svds只计算前若干个奇异值。去噪场景下我们关心的是大奇异值尾部那些小奇异值本来就会被软阈值压成0。这个改动在大数据上的加速效果非常可观也是我踩了好几次内存坑之后才学乖的。这篇分享就到这里。如果你正好在调类似的张量去噪代码遇到收敛、内存或者参数选择问题欢迎顺着这些经验先排查一轮。本文还有配套的精品资源点击获取
RELATED — 相关阅读

相关资讯

LATEST — 最新资讯

最新发布

TODAY — 本日精选

新闻

WEEKLY — 本周精选

新闻

MONTHLY — 本月精选

新闻