简介这份Matlab代码资源面向从事高维数据分析、信号处理与机器学习的研究者和工程师聚焦于用凸优化方法实现张量分解从而完成数据恢复与去噪任务。张量作为矩阵向高维的推广能有效刻画图像、视频、多模态数据的内在结构而凸优化则保证低秩近似求解的稳定性与全局最优性适合具备一定Matlab基础、希望深入张量处理与优化技术的读者。资源包共48个文件以45个m脚本为核心辅以README.md说明文档及git配置项整体约55KB涵盖张量分解、ADM求解、软阈值、矩阵补全、去噪与补全实验及绘图对比等模块目录组织清晰。目前已有738人学习下载。读者可借助其中的算法函数、示例脚本与测试用例快速理解CP、Tucker等分解思路并依据自身数据类型调整迭代次数与阈值完成从数据预处理到去噪效果可视化的完整流程。1. 从一段被噪声淹没的信号说起这套张量分解去噪代码能帮你做什么如果你手头有一段多通道采集的信号——比如阵列传感器输出、高光谱波段数据、或者某个实验平台同时记录的多个通道——大概率会遇到同一个问题单通道看还行多通道拼在一起就全是噪声传统的一维滤波按通道逐个处理处理完通道之间的相关性也被切断了。这套基于凸优化的张量分解 Matlab 代码解决的正是这个场景它把多通道数据组织成高阶张量用张量分解的方式把低秩的信号子空间和稀疏的噪声分离开而不是逐通道去卷滤波核。代码是完整的 Matlab 工程包含分解主流程、凸优化求解和去噪演示脚本适合做信号处理、阵列数据处理、图像去噪方向的研究生和工程师直接拿来跑也适合在此基础上做二次开发。它不依赖深度学习工具箱纯数值计算路线对没有 GPU 的环境很友好。2. 技术选型与整体架构为什么用张量分解而不是逐通道滤波2.1 张量分解去噪的基本思路先把问题说清楚。假设你采集到的数据是一个三阶张量 $\mathcal{X} \in \mathbb{R}^{I_1 \times I_2 \times I_3}$其中两个维度是空间或通道维度第三个维度是时间或样本维度。真实的干净信号 $\mathcal{L}$ 通常具有低秩结构——因为多通道之间的信号是高度相关的而噪声 $\mathcal{E}$ 往往是稀疏的或者满秩的。于是去噪问题可以写成$$\min_{\mathcal{L}, \mathcal{E}} \ \text{rank}(\mathcal{L}) \lambda |\mathcal{E}|_0 \quad \text{s.t.} \quad \mathcal{X} \mathcal{L} \mathcal{E}$$这个问题的难点在于 rank 和 L0 范数都是非凸的直接求解是 NP 难的。工程上的常见做法是松弛成凸问题用张量的核范数tensor nuclear norm替代 rank用 L1 范数替代 L0得到$$\min_{\mathcal{L}} \ |\mathcal{L}|_* \lambda |\mathcal{X} - \mathcal{L}|_1$$这就是这套代码的核心优化目标。它属于凸优化范畴可以用 ADMM交替方向乘子法或者近端梯度法稳定求解。相比逐通道小波阈值去噪这种做法的优势在于分解过程天然保留了通道间的耦合关系去噪后的张量在结构上更一致不会出现某个通道被过度平滑、另一个通道还残留噪声的情况。2.2 代码的模块划分拿到源码包后目录结构大致是这样的不同版本可能略有差异但核心模块是固定的tensor-denoise/ ├── main_denoise.m % 主入口加载数据、调参、输出去噪结果 ├── core/ │ ├── tensor_nuclear_norm.m % 张量核范数计算基于 t-SVD │ ├── tsvd.m % 张量奇异值分解实现 │ ├── soft_threshold.m % 软阈值算子L1 近端算子 │ └── admm_solver.m % ADMM 主循环 ├── utils/ │ ├── tensor2mat.m % 张量展开/折叠工具 │ ├── psnr_ssim.m % 去噪评价指标 │ └── gen_test_data.m % 生成带噪声的测试张量 └── data/ └── sample.mat % 示例数据这个划分很清晰core/放算法核心utils/放工具函数main_denoise.m是唯一需要你手动改的入口。我一般拿到这种结构的项目第一件事就是打开main_denoise.m看它默认加载什么数据、参数怎么设然后直接跑一遍看输出。2.3 凸优化求解器的选择代码里用的是 ADMM而不是 CVX 这类建模工具箱。原因很实际CVX 虽然写起来简洁但每次求解都要重新解析模型对于需要反复迭代调参的去噪任务来说太慢而且 CVX 对三阶以上张量的支持并不好。ADMM 把原问题拆成两个子问题交替求解每个子问题都有闭式解迭代几十次就能收敛速度上完全可接受。ADMM 的增广拉格朗日形式为$$\mathcal{L}\rho(\mathcal{L}, \mathcal{E}, \mathcal{Y}) |\mathcal{L}|* \lambda |\mathcal{E}|_1 \langle \mathcal{Y}, \mathcal{X} - \mathcal{L} - \mathcal{E} \rangle \frac{\rho}{2} |\mathcal{X} - \mathcal{L} - \mathcal{E}|_F^2$$其中 $\mathcal{Y}$ 是拉格朗日乘子张量$\rho$ 是惩罚参数。迭代更新规则在admm_solver.m里实现下面会逐段拆。3. 核心模块拆解从 t-SVD 到 ADMM 迭代的完整链路3.1 张量奇异值分解t-SVD的实现张量核范数的计算依赖 t-SVD。和矩阵 SVD 不同t-SVD 是在傅里叶域沿第三维做矩阵 SVD然后逆变换回来。tsvd.m的核心逻辑如下function [U, S, V] tsvd(X) % 输入 X: n1 x n2 x n3 张量 % 输出 U, S, V: t-SVD 分解结果 [n1, n2, n3] size(X); % 沿第三维做 FFT转到傅里叶域 Xf fft(X, [], 3); U zeros(n1, n1, n3); S zeros(n1, n2, n3); V zeros(n2, n2, n3); for i 1:n3 % 对每个频率切片做矩阵 SVD [Uf, Sf, Vf] svd(Xf(:,:,i), econ); U(:,:,i) Uf; S(:,:,i) Sf; V(:,:,i) Vf; end % 逆 FFT 回到原域 U ifft(U, [], 3); S ifft(S, [], 3); V ifft(V, [], 3); end这里有几个关键点。第一fft(X, [], 3)是沿第三维做变换这个维度通常对应时间或样本变换后每个频率切片是一个独立的矩阵。第二svd(..., econ)用经济型分解因为张量切片往往不是方阵经济型能避免产生大量零奇异值。第三逆变换后得到的 U、S、V 是复数张量但实际使用时通常只取实部或者模值因为原始数据是实数的。注意如果你的 Matlab 版本较老R2016 之前fft对三维数组的支持可能有问题需要手动 reshape 成二维再处理。常见做法是用reshape把n1 x n2 x n3变成n1 x (n2*n3)做完 FFT 再变回来。3.2 软阈值算子与核范数近端算子ADMM 的每个子问题都涉及近端算子。对于 L1 范数近端算子就是软阈值function Y soft_threshold(X, tau) % 软阈值算子: Y sign(X) * max(|X| - tau, 0) % X: 输入张量/矩阵 % tau: 阈值参数对应 ADMM 中的 lambda/rho Y sign(X) .* max(abs(X) - tau, 0); end这个函数看起来简单但tau的取值直接决定去噪强度。tau太大信号细节被一起抹掉tau太小噪声残留明显。代码里默认tau lambda / rho其中lambda是用户在主脚本里设的稀疏惩罚系数rho是 ADMM 惩罚参数。我一般会先固定rho 1然后在一个对数尺度上扫lambda比如 0.001 到 1看哪个值对应的 PSNR 最高。对于核范数近端算子是奇异值软阈值Singular Value Thresholding, SVTfunction L svt_tensor(X, tau) % 张量奇异值软阈值: 对奇异值做软阈值后重构 [U, S, V] tsvd(X); S_thresh soft_threshold(S, tau); L tprod(U, tprod(S_thresh, ttranspose(V))); end这里用到了tprod张量乘积和ttranspose张量转置它们是 t-SVD 体系下的基本运算。tprod的实现本质是在傅里叶域做矩阵乘法再逆变换和 t-SVD 的推导保持一致。3.3 ADMM 主循环的逐段解析admm_solver.m是整个项目的心脏。它的输入是观测张量X、参数lambda、rho、最大迭代次数max_iter和收敛阈值tol输出是去噪后的张量L和噪声张量E。function [L, E, obj_history] admm_solver(X, lambda, rho, max_iter, tol) % ADMM 求解: min ||L||_* lambda*||E||_1 s.t. X L E [n1, n2, n3] size(X); % 初始化 L zeros(n1, n2, n3); E zeros(n1, n2, n3); Y zeros(n1, n2, n3); % 拉格朗日乘子 obj_history zeros(max_iter, 1); for k 1:max_iter % --- L 子问题: 核范数近端算子 --- % L argmin ||L||_* (rho/2)||X - L - E Y/rho||_F^2 temp X - E Y / rho; L svt_tensor(temp, 1 / rho); % --- E 子问题: L1 近端算子 --- % E argmin lambda*||E||_1 (rho/2)||X - L - E Y/rho||_F^2 temp X - L Y / rho; E soft_threshold(temp, lambda / rho); % --- 乘子更新 --- Y Y rho * (X - L - E); % --- 收敛判断 --- primal_res norm(X(:) - L(:) - E(:)); dual_res rho * norm(L(:) - L_old(:)); % 需要保存上一步的 L obj_history(k) tensor_nuclear_norm(L) lambda * sum(abs(E(:))); if primal_res tol dual_res tol fprintf(ADMM 在第 %d 次迭代收敛\n, k); obj_history obj_history(1:k); break; end L_old L; end end逐段说明。L 子问题里temp X - E Y/rho是增广拉格朗日项对 L 的梯度方向然后调用svt_tensor做奇异值软阈值阈值取1/rho。E 子问题同理temp X - L Y/rho然后做元素级软阈值阈值取lambda/rho。乘子更新Y Y rho*(X - L - E)是标准的对偶上升步骤。收敛判断用了原始残差和对偶残差两个指标。原始残差衡量约束满足程度对偶残差衡量 L 的变化量。两个都小于tol才认为收敛。实际跑的时候tol设1e-4到1e-5比较合适太小会迭代很多次但精度提升有限。提示rho的取值对收敛速度影响很大。rho太大L 子问题的阈值1/rho太小去噪不充分rho太小E 子问题的阈值lambda/rho太大噪声被过度剔除。常见做法是设rho 1然后调lambda。如果收敛慢可以试试自适应rho每 10 次迭代根据残差比调整。3.4 主脚本的调用流程main_denoise.m把上面这些模块串起来% 加载数据 load(data/sample.mat); % 假设里面有变量 X_clean 和 X_noisy % 参数设置 lambda 0.05; % 稀疏惩罚系数 rho 1; % ADMM 惩罚参数 max_iter 200; % 最大迭代次数 tol 1e-5; % 收敛阈值 % 调用 ADMM 求解 [L_denoised, E_noise, obj_hist] admm_solver(X_noisy, lambda, rho, max_iter, tol); % 评价指标 psnr_val psnr_ssim(X_clean, L_denoised); fprintf(去噪后 PSNR: %.2f dB\n, psnr_val); % 可视化 figure; subplot(1,3,1); imshow(X_noisy(:,:,1), []); title(含噪); subplot(1,3,2); imshow(L_denoised(:,:,1), []); title(去噪); subplot(1,3,3); imshow(X_clean(:,:,1), []); title(干净);这段脚本可以直接跑。如果你的数据不是图像格式把imshow换成plot或者imagesc就行。psnr_ssim是工具函数输入干净张量和去噪张量输出 PSNR 和 SSIM 两个指标。4. 跑起来环境配置、数据准备与参数调优4.1 环境要求与依赖检查这套代码对 Matlab 版本的要求不高R2015b 以上基本都能跑。不需要任何额外工具箱——这是它比很多深度学习去噪代码方便的地方。如果你用的是 Octave大部分函数也能兼容但fft对三维数组的处理和 Matlab 略有差异可能需要手动调整tsvd.m里的维度参数。跑之前先检查两件事一是core/和utils/是否都在 Matlab 的搜索路径里二是data/sample.mat是否存在。如果路径不对main_denoise.m会在load那一步报错。常见做法是在脚本开头加一行addpath(genpath(pwd))把当前目录及子目录都加进路径。4.2 用自己的数据替换示例数据示例数据sample.mat里通常包含两个变量X_clean干净张量和X_noisy加噪张量。你要用自己的数据只需要保证变量名一致或者改main_denoise.m里的加载逻辑。数据格式要求是三维数组n1 x n2 x n3。如果你手头是二维矩阵比如单张图像需要先扩展成三维。常见做法是用repmat复制几份或者用reshape把多个样本拼成第三维。比如你有 100 张 256x256 的图像可以拼成256 x 256 x 100的张量。% 把多个二维样本拼成三维张量 data zeros(256, 256, 100); for i 1:100 data(:,:,i) load(sprintf(sample_%03d.mat, i)); end X_noisy data;如果你的数据维度更高比如四阶张量这套代码需要改tsvd.m和tprod.m的维度处理逻辑。三阶是这套代码的原生支持范围四阶以上需要额外做 matricization。4.3 参数调优的实操方法lambda是最关键的参数。它控制稀疏项和低秩项之间的权衡。lambda越大越多能量被归到噪声项 E 里去噪越激进lambda越小越多噪声被保留在 L 里。我一般用网格搜索先设lambda从 0.001 到 0.5 取 10 个对数等距值每个值跑一遍 ADMM记录 PSNR。然后在这个最优值附近再细化搜索。整个过程大概几分钟比盲目试参数靠谱得多。lambda_list logspace(-3, log10(0.5), 10); psnr_list zeros(size(lambda_list)); for i 1:length(lambda_list) [L, ~, ~] admm_solver(X_noisy, lambda_list(i), 1, 200, 1e-5); psnr_list(i) psnr_ssim(X_clean, L); end [~, best_idx] max(psnr_list); fprintf(最优 lambda: %.4f, PSNR: %.2f dB\n, lambda_list(best_idx), psnr_list(best_idx));rho一般不用调固定 1 就行。但如果你的数据量很大比如n1*n2*n3 1e6ADMM 收敛会变慢可以把rho设大一点比如 5 或 10同时把max_iter加到 500。4.4 运行时间与内存估算三阶张量的大小直接决定内存占用。tsvd需要对n3个切片做矩阵 SVD每个切片是n1 x n2所以内存峰值大约是O(n1*n2*n3)的几倍。对于256 x 256 x 100的数据单精度下大约 25 MB双精度 50 MB普通笔记本完全能扛住。但如果你的数据是512 x 512 x 500内存会到 GB 级别需要改成分块处理或者用single精度。运行时间方面一次 ADMM 迭代的主要开销在tsvd和tprod。256 x 256 x 100的数据一次迭代大约 0.5 到 1 秒收敛通常需要 50 到 150 次迭代总时间 1 到 3 分钟。如果开了parfor并行化tsvd里的频率切片循环能快 2 到 4 倍。5. 避坑与排查跑这套代码时最容易翻车的几个地方5.1 现象去噪结果全黑或全白PSNR 极低原因通常是lambda设得太大或太小。lambda太大时E 子问题把几乎所有能量都当成噪声剔除L 接近零张量lambda太小时E 接近零L 几乎等于原始含噪数据去噪没效果。解决方法是先跑一遍lambda扫描看 PSNR 曲线的峰值位置。如果曲线单调下降说明lambda的搜索范围整体偏大把上限调小如果单调上升说明范围偏小把上限调大。5.2 现象ADMM 不收敛迭代到 max_iter 还在震荡常见原因是rho和lambda的比例不合适。ADMM 的收敛性对rho敏感rho太小会导致对偶变量更新步长过大产生震荡。解决办法是增大rho比如从 1 调到 5或者实现自适应rho每 20 次迭代检查原始残差和对偶残差的比值如果比值大于 10 就增大rho小于 0.1 就减小rho。代码里可以加一个简单的自适应逻辑if primal_res 10 * dual_res rho rho * 2; Y Y / 2; % 乘子需要同步缩放 elseif dual_res 10 * primal_res rho rho / 2; Y Y * 2; end5.3 现象tsvd报错「矩阵维度不一致」这通常是因为输入张量的第三维大小n3和tsvd内部假设的不一致。有些版本的tsvd.m假设n3是奇数因为 FFT 后频率切片是对称的只处理一半。如果你的n3是偶数需要检查代码里是否有floor(n3/2)之类的硬编码。解决办法是打开tsvd.m看for i 1:n3这个循环的范围。如果是1:ceil(n3/2)说明它只处理了前半部分频率后半部分需要共轭对称补全。这种情况下要么把数据第三维补零到奇数要么改代码支持偶数。5.4 现象去噪后的张量出现块状伪影这是张量分解类方法的典型问题根源在于 t-SVD 沿第三维做 FFT 时默认数据在第三维是周期性的。如果你的数据第三维不是周期信号比如时间序列有趋势项FFT 会在边界产生频谱泄漏导致重构后出现块状伪影。解决办法有两个一是去趋势在跑 ADMM 之前先对第三维做差分或减去均值二是加窗在 FFT 之前乘一个 Hann 窗但这样会损失一部分能量需要做重叠保留处理。我一般用第一种简单有效。5.5 现象换用自己的数据后 PSNR 计算报错psnr_ssim.m通常假设输入是uint8类型且范围在 0 到 255 之间。如果你的数据是double类型且范围在 0 到 1直接传进去会得到错误的 PSNR 值。解决办法是在调用前做类型转换和范围归一化if isa(X_clean, double) max(X_clean(:)) 1 X_clean uint8(X_clean * 255); L_denoised uint8(L_denoised * 255); end psnr_val psnr_ssim(X_clean, L_denoised);6. 二次开发把 ADMM 换成近端梯度法以及一个我常用的验证习惯这套代码的架构留了很好的扩展点。admm_solver.m是一个独立的求解器输入输出接口清晰你完全可以写一个pgd_solver.m用近端梯度法替代它然后在main_denoise.m里换一行调用就行。近端梯度法的迭代格式是$$\mathcal{L}{k1} \text{prox}{\eta |\cdot|_*}\left(\mathcal{L}_k - \eta \nabla g(\mathcal{L}_k)\right)$$其中 $g(\mathcal{L}) \lambda |\mathcal{X} - \mathcal{L}|_1$梯度 $\nabla g$ 在 L1 范数下是次梯度实际实现时用sign函数近似。步长 $\eta$ 需要满足 $\eta 2 / L_g$其中 $L_g$ 是 $g$ 的 Lipschitz 常数。对于 L1 项$L_g \lambda$所以 $\eta 2/\lambda$。function [L, obj_history] pgd_solver(X, lambda, eta, max_iter, tol) % 近端梯度法求解: min ||L||_* lambda*||X - L||_1 L X; % 初始化 obj_history zeros(max_iter, 1); for k 1:max_iter % 梯度步: g(L) lambda * ||X - L||_1 grad lambda * sign(L - X); L_half L - eta * grad; % 近端步: 核范数近端算子 L_new svt_tensor(L_half, eta); % 收敛判断 obj_history(k) tensor_nuclear_norm(L_new) lambda * sum(abs(X(:) - L_new(:))); if norm(L_new(:) - L(:)) / norm(L(:)) tol obj_history obj_history(1:k); break; end L L_new; end end近端梯度法的优势是实现简单不需要维护对偶变量劣势是收敛速度通常比 ADMM 慢尤其是当问题条件数不好时。我一般会先用 ADMM 跑出结果作为基准再用 PGD 验证如果两者 PSNR 差距在 0.5 dB 以内说明 ADMM 的结果可信。另一个常用的改造方向是换用不同的张量分解。t-SVD 只是其中一种还有 CP 分解、Tucker 分解。CP 分解的秩定义更直观但计算 NP 难通常用 ALS 交替最小二乘近似。Tucker 分解更灵活但核范数的定义更复杂。如果你想试 CP 分解可以把svt_tensor换成 CP-ALS 的软阈值版本但要注意 CP 分解的秩选择比 t-SVD 更敏感。验证改造是否成功我有一个固定习惯不管换什么求解器或分解方式都先在gen_test_data.m生成的合成数据上跑一遍。合成数据的好处是干净张量已知PSNR 可以直接算而且噪声类型可控高斯、拉普拉斯、椒盐。如果合成数据上 PSNR 达标再换真实数据。真实数据没有 ground truth只能靠视觉检查和下游任务指标比如分类准确率来判断。从那以后我每次改求解器都强制走一遍「合成数据验证 → 真实数据试跑 → 参数微调」这个流程省了很多来回折腾的时间。希望帮到你。本文还有配套的精品资源点击获取