简介：这套源代码基于凸优化理论实现张量分解与去噪，面向从事图像处理、视频分析、社交网络和多模态数据挖掘的研究者与工程师，解决高维数据中噪声去除与信号恢复问题。代码包共48个文件，其中45个.m脚本构成主体，涵盖张量1范数优化、交替方向乘子法求解、随机张量生成、去噪与补全实验、误差对比绘图等功能模块，另含README说明及Git工程配置，压缩包仅55KB，轻量且易于修改复用。目前已有738人学习下载，适合具备Matlab基础和凸优化背景的读者进阶使用。借助示例与测试脚本，可深入理解CP、Tucker等低秩近似在去噪中的具体实现，掌握迭代次数、阈值参数对恢复效果的影响，并迁移至自身科研或工程场景，是探索张量恢复与优化方法的实用开源工具。

1. 去噪与补全：为什么我把这个张量分解项目留在工具箱里

如果你处理过视频去噪或图像补全，大概率遇到过这样的场景：一段视频里混入了随机椒盐噪声，用矩阵化的方法一帧一帧去噪，结果把时间维上的连续性完全丢掉了，噪声是少了，画面却也变得闪烁。这个基于凸优化的张量分解Matlab代码，正是为解决这类问题设计的，它把带噪声、带缺失的数据当作一个三维张量整体处理，一次性完成去噪和补全两项任务，而不是逐帧做矩阵分解。对正在做图像恢复、视频修复、多维信号处理相关研究的人来说，这份资源最大的价值是给了你一套完整的、可复现的工程实现：核心算法基于张量核范数与交替方向乘子法（ADMM），不依赖某个闭源库，拿到手就能在Matlab里跑通并继续改造。

开篇说清楚定位：这是一个偏学术风格、但代码组织非常工整的Matlab实现，适合用来做算法对比实验、论文复现和二次开发。接下来我按实际拆代码的顺序，从目录结构、核心模块、运行方式到改造建议，一步步把它讲透。

2. 项目定位与技术选型：凸优化、张量分解和工程落地的平衡

2.1 为什么用张量分解而不是矩阵分解

传统去噪方法把一段视频、一组多光谱图像或一批传感器数据看作多个独立矩阵，分别做低秩近似。这样做的问题在于忽略了各切片之间的关联性。比如监控视频的前后帧之间、高光谱图像的相邻波段之间，数据本身存在强相关，矩阵化处理相当于把这些关联切断，导致去噪结果在时间或光谱维度上出现跳变。这个项目采用张量分解，本质上是把数据整体建模为一个三阶张量，利用张量低秩性去捕获跨维度的全局结构，这使得它在处理强相关多维数据时，效果明显好于逐帧矩阵低秩逼近。

张量分解的具体形式有CP分解、Tucker分解和张量奇异值分解（t-SVD）等，本项目走的是凸优化路线：构造一个关于张量核范数的目标函数，通过最小化该目标函数来寻找低秩近似。相比非凸的张量分解方法，凸优化方法的优势是有全局最优解作为理论支撑，迭代过程稳定，不容易陷入糟糕的局部极小值。当然，代价是迭代次数更多，计算量更大，这一点在大型三维数据上需要通过参数调节来平衡。

2.2 技术栈全貌：Matlab生态、ADMM求解器与t-SVD核心

梳理一下这个项目涉及的核心技术栈，我在拆解时大致归纳为以下四层，每一层都有清晰的替换和扩展空间：

ADMM在这里承担的是“调度者”角色。原始问题是一个带线性约束的最小化问题，直接求解析解几乎不可能。ADMM把它拆成几个容易求解的子问题：一个是带张量核范数的低秩逼近子问题，一个是带线性约束的保真子问题，再加上对偶变量的更新。每个子问题都有闭式解或成熟求解路径，所以迭代体面简洁，调参也直观。

从工程角度看，Matlab选择得非常务实。Matlab的矩阵运算与FFT函数在张量处理上写得非常高效，代码表达与论文公式高度一致，能直接把算法公式映射成m文件，这在原型验证阶段比Python加NumPy的组合更顺手，尤其适合需要快速调整算法细节并立即看到数值结果的场景。

3. 核心模块拆解：从数据构造到ADMM迭代求解的完整链路

3.1 目录结构与模块职责速览

这个项目的目录结构不是那种把几十个m文件扔在同一层的仓库，而是做了清晰的功能分区，我拆解后简化为以下框架：

tensor-completion-denoising/ ├── main_completion.m # 主入口脚本，负责数据读取、参数设置、调用求解 ├── optimize/ │ ├── admm_solver.m # ADMM主循环，各子问题的交替更新 │ ├── update_X.m # 求解带张量核范数的低秩逼近子问题 │ ├── update_E.m # 更新稀疏噪声项 │ └── update_M.m # 更新观测约束项 ├── tensor_ops/ │ ├── t_svd.m # 张量奇异值分解核心实现 │ ├── t_product.m # 张量t-product运算 │ └── tensor_norm.m # 张量核范数与Frobenius范数计算 ├── data/ │ ├── gen_noisy_data.m # 生成带噪声与缺失的合成数据 │ └── load_real_data.m # 读取真实图像或视频数据 └── metrics/ └── cal_psnr_rse.m # 计算PSNR与相对误差RSE

主入口脚本main_completion.m相当于项目的“总控台”，它把数据准备、参数初始化、算法调用和结果评估串成一条流水线。我建议你从它开始阅读代码，顺着函数的调用关系就能理清整个算法流程。optimize目录下的几个函数是核心求解逻辑，对应ADMM迭代中的各个更新步骤；tensor_ops目录提供的是支撑运算，尤其t_svd.m是整个算法的基础组件；data与metrics目录则分别解决实验数据从哪来以及去噪效果如何量化的问题。

3.2 合成数据的生成逻辑：可控噪声与缺失模式

要验证去噪算法，第一步是先有一个“真值张量”和叠加噪声后的观测张量。gen_noisy_data.m完成的是这样一个过程：它先构造一个低秩的三阶张量作为干净数据，然后按指定比例随机添加高斯噪声或椒盐噪声，再随机挖去一定比例的元素来模拟缺失场景。

function [Y_obs, M_mask, X_true] = gen_noisy_data(sz, rank, noise_ratio, miss_ratio) % 生成合成张量数据：低秩真值 + 噪声 + 缺失掩码 % sz: 张量维度 [n1, n2, n3] % rank: 真值张量的CP秩 % noise_ratio: 噪声强度系数 % miss_ratio: 缺失比例 [0,1] rng(2025); % 固定随机种子保证可复现 X_true = randn(sz(1), sz(2), rank) * randn(rank, sz(3)) * 0.1; X_true = reshape(X_true, sz); % 构造低秩真值张量 Y_noise = X_true + noise_ratio * randn(sz); % 叠加高斯噪声 M_mask = double(rand(sz) > miss_ratio); % 生成缺失掩码，1表示观测到 Y_obs = Y_noise .* M_mask; % 观测张量：保留非缺失位置 end

这段代码展示了数据构造中的几个关键设计：真值张量先通过两个随机矩阵的乘积生成低秩结构，这样我们确切知道它的真实秩是多少，为后续算法的恢复精度提供了对比基准；噪声强度用noise_ratio系数控制，实验时逐个调节这个参数就能观察算法在不同信噪比条件下的表现；缺失掩码M_mask是独立的0-1矩阵，1代表该位置有观测值，0代表缺失。观测张量Y_obs通过逐元素乘法获得，缺失位置直接置零，这与后续ADMM优化中的约束条件设计刚好对应。

参数说明很直接，但这里有一个工程细节值得留意：如果sz(3)表示视频帧数，那么rank参数一般取很小的值比如10到30。取值太大，真值张量的低秩特性不明显，算法几乎无从区分信号和噪声；取值太小，又可能让真值过于简单，导致算法“胜之不武”。实践上我会先取rank=20做基准实验，再调整到10和30看结果的敏感度。

3.3 ADMM主循环：一个迭代框架如何协调三个子问题

看完数据构造，接下来是项目的算法心脏——admm_solver.m。ADMM求解的凸优化问题形式如下：目标函数是张量核范数加上稀疏项L1范数，同时要求恢复出的张量在观测位置上与输入数据一致。求解思路是把这个问题拆成带等式约束的形式，然后对每个子问题轮换求解并更新对偶变量。

function [X_out, E_out, history] = admm_solver(Y_obs, M_mask, params) % 基于ADMM的张量去噪与补全主循环 % params: 结构体，包含lambda, rho, max_iter, tol等超参数 lambda = params.lambda; % 稀疏噪声项权重 rho = params.rho; % ADMM惩罚参数 max_iter = params.max_iter; % 最大迭代次数 tol = params.tol; % 相对变化收敛阈值 [n1, n2, n3] = size(Y_obs); X = Y_obs; % 待恢复张量，用观测数据初始化 M = Y_obs; % 引入的辅助变量 E = zeros(n1, n2, n3); % 稀疏噪声项 U = zeros(n1, n2, n3); % 对偶变量 for iter = 1:max_iter % 子问题1：求解带张量核范数的低秩逼近 X = update_X(M - E + U / rho, lambda / rho); % 子问题2：更新稀疏噪声项，带L1范数约束 E = update_E(M - X + U / rho, lambda / rho); % 子问题3：更新辅助变量M，强制满足观测约束 M = (Y_obs + rho * (X + E - U / rho)) / (1 + rho); M = M .* M_mask + X .* (1 - M_mask); % 对偶变量更新 U = U + rho * (X + E - M); % 检查收敛条件 history.rel_change(iter) = norm(X - X_prev, 'fro') / (norm(X_prev, 'fro') + eps); if history.rel_change(iter) < tol break; end X_prev = X; end X_out = X; E_out = E; end

这个主循环的巧妙之处在于它把一个复杂的约束优化问题分解成三个相对容易的子问题。第一步update_X是核心，它执行的是张量核范数最小化，本质上是张量域的软阈值操作，在后续小节详讲。第二步update_E处理稀疏噪声项，使用软阈值收缩，直观理解就是把绝对值小于阈值的位置置零，相当于不断“剔除”那些幅度较小的噪声成分。第三步update_M则保证最终结果在观测位置和输入数据一致，在缺失位置则完全由模型推断，这正是补全功能的来源。

收敛判定用的是张量间的Frobenius范数相对变化率，当前后两次迭代结果几乎不再变化时，就认为找到了稳定解。这里需要注意rho的选取：rho太小，对偶变量更新缓慢，收敛很慢；rho太大，前期震荡明显甚至发散。我一般在L曲线搜索时从1e-3起步，如果发现相对变化曲线出现锯齿状，就增大rho到1e-2重新跑一遍。

3.4 张量核范数与t-SVD：追踪核心运算的实现细节

张量核范数是在凸优化框架下衡量张量“秩”的关键工具，本项目的t_svd.m实现方式直接决定了算法能否在可接受时间内完成迭代。t-SVD的核心思想是对张量的第三维做傅里叶变换，然后在频域逐片执行矩阵SVD，最后通过逆傅里叶变换回到时域。

function [U, S, V] = t_svd(A) % 张量奇异值分解，基于傅里叶变换与逐片矩阵SVD % 输入A: n1 x n2 x n3 实张量 % 输出U, S, V: t-SVD分解的三个因子张量 [n1, n2, n3] = size(A); Af = fft(A, [], 3); % 沿第三维做FFT U = zeros(n1, min(n1,n2), n3, 'like', A); S = zeros(min(n1,n2), min(n1,n2), n3, 'like', A); V = zeros(n2, min(n1,n2), n3, 'like', A); % 对每个频域切片的矩阵执行标准SVD for k = 1:n3 [U(:, :, k), S(:, :, k), V(:, :, k)] = svd(Af(:, :, k), 'econ'); end % 逆FFT回到时域，并做后处理 U = ifft(U, [], 3); S = ifft(S, [], 3); V = ifft(V, [], 3); end

这段实现充分利用了傅里叶域里卷积定理带来的对角化特性：三维张量的t-product在频域变成逐片矩阵乘法，因此张量SVD可以分解成n3个独立的矩阵SVD，这正是并行化的突破口。svd(..., 'econ')使用经济型SVD，截断到矩阵的最小维度上，因为张量核范数求和只需要奇异值向量而不需要补零，能降低存储开销。

在实现中有两个常被忽略的细节。一是在k=1与对应共轭对称的切片上，FFT结果具有特殊结构，如果不做对称处理，逆FFT后会出现轻微的数值虚部噪声；严谨的做法是先对U、S、V做一次共轭对称修正，或者使用real()截断。另一个是执行SVD前对数据做归一化以避免数值溢出。我在最初运行这个函数时曾因为输入数据量级在1e4~1e5之间没有归一化，导致部分频点奇异值计算准确度下降，后续迭代出现NaN。加上数据归一化处理后，整个流程才稳定下来。

3.5 指标评估模块：PSNR与RSE不是只看一个数字

跑通算法只是第一步，判断恢复效果好坏还得靠定量的指标评估。metrics/cal_psnr_rse.m实现了两个最常用的图像恢复指标：峰值信噪比PSNR与相对误差RSE。PSNR基于最大像素值与均方误差之比，单位是dB，值越大说明恢复越接近真值；RSE则直接衡量恢复张量与真值之间的相对差异，值越接近0说明效果越好。

function [PSNR_val, RSE_val] = cal_psnr_rse(X_est, X_true) % 计算峰值信噪比与相对误差 % X_est: 算法恢复的完整张量 % X_true: 真实无噪声张量 diff = X_est - X_true; mse_val = mean(diff(:).^2); % 均方误差 max_val = max(X_true(:)); PIXEL_MAX = max_val; % 像素最大值，通常取255或1 PSNR_val = 10 * log10(PIXEL_MAX^2 / (mse_val + eps)); RSE_val = norm(diff(:)) / norm(X_true(:)); % 相对误差 end

使用这个模块时我会同时记录PSNR和RSE两个数，因为它们反映的效果维度不同：PSNR对噪声的随机误差敏感，适合横向对比不同算法在同一噪声水平下的性能；RSE则能看出恢复结果的整体结构保真度，比如纹理细节是否被过度平滑。一个常见做法是设置一个固定噪声梯度序列noise_ratio = [0.05, 0.1, 0.2, 0.4]，分别记录指标后绘制曲线，这样可以看到算法性能随噪声强度增加是如何衰减的，帮助判断它的鲁棒边界。

4. 运行复现与参数调试：从第一个demo到批量实验

4.1 环境准备与启动顺序

运行这个项目不需要额外的Matlab工具箱，它的核心计算基于内置的fft、svd和基本矩阵运算，这一点对复现非常友好。环境要求其实很低，R2018之后的Matlab版本都能直接运行，因为用到了结构体数组的扩展语法和字符串处理函数，早年版本可能需要小幅适配。推荐的启动顺序是先在Matlab命令行进入项目根目录，然后运行main_completion.m，它会自动调用data、optimize、metrics三个目录下的函数。

启动前建议手动检查两个环境细节。第一是确认当前目录下optimize和tensor_ops子目录是否在Matlab搜索路径中，如果直接运行出现“函数或变量无法识别”的错误，用addpath命令把子目录加进搜索路径即可。第二是运行gen_noisy_data.m时，如果Matlab提示随机数种子相关的警告，把rng(2025)换成rng('default')就能消除。

4.2 逐步解释：params超参数的定位与调优

admm_solver.m接收的params结构体是调参的核心入口，它在main_completion.m中设置。以下是完整的参数表及其调节建议，这部分我整理成了可直接对照的清单：

调参的顺序有讲究。我一般先固定lambda=0.01，调整rho找收敛稳定的区间；之后固定rho，用lambda在0.001到0.1之间做一次网格搜索，观察PSNR和RSE的曲线变化。lambda过小时稀疏噪声项约束形同虚设，去噪不彻底；lambda过大又会过度收缩信号本身，导致恢复结果偏模糊，这在图像上表现为边缘被磨平。建议每次改完参数后跑通一个demo，计算时间通常在一到两分钟，如果你发现一次跑通超过了五分钟，检查是否在FFT循环里加了不必要的矩阵复制操作。

4.3 可视化验证：查看恢复结果的两种方式

算法跑完，除了看PSNR数字，还要用眼睛检查恢复效果。main_completion.m末尾会输出一个简单的可视化结果：分别展示原始真值张量的中间切片、带噪声的观测切片和算法恢复的切片，这三张图并排放置可以直观看到去噪和补全效果。

% 在main_completion.m中查看某一切片的恢复效果 k = 10; % 选择第10帧（或第10个切片） figure('Name', '去噪/补全结果对比'); subplot(1,3,1); imshow(X_true(:,:,k), []); title('真值'); subplot(1,3,2); imshow(Y_obs(:,:,k), []); title('观测（含噪声与缺失）'); subplot(1,3,3); imshow(X_out(:,:,k), []); title('恢复结果');

关于可视化有两个值得留意的点。一是imshow的第二个参数空矩阵[]表示自动缩放灰度范围，让每个切片独立显示其动态范围，这样对于对比算法更友好，也方便观察谱数据中的细节结构。二是当数据不是图像而是传感器阵列数据时，把三维张量按第三维展开成矩阵热力图来观察，heatmap函数比imshow更直观。我在跑高光谱图像去噪实验时，习惯查看第1、第50、第100个波段的恢复对比，观察不同频带上算法表现的差异，这比单看一个切片的PSNR更能判断算法的适应性。

5. 避坑与常见问题排查：张量代码里那些防不胜防的坑

5.1 报错“矩阵维度必须一致”：问题多半出在掩码与观测张量的维度匹配上

现象：运行gen_noisy_data.m或主流程时，Matlab直接在逐元素乘法或加法步骤报错，提示矩阵维度必须一致。

原因：这类错误最常见的情况是用矩阵掩码去乘三维张量，或是在构建掩码时尺寸少了第三维。例如M_mask写成二维的rand(n1, n2)，却对sz=[n1, n2, n3]的三维数据做逐元素乘法。

解决：统一所有张量的维度为三维。构造掩码时显式使用rand(sz)而不是rand(n1, n2)，或者在乘法前用repmat(M_mask, [1, 1, n3])扩充到第三维，维度对齐后再进入后续计算。

5.2 迭代中出现NaN：多与未归一化和SVD数值溢出有关

现象：ADMM迭代进行到中途，history.rel_change开始出现NaN，恢复结果变成全零或乱码。

原因：这类问题大多数时候不是算法本身发散，而是子问题中SVD计算的输入量级过大或过小，导致在计算奇异值时出现数值越界。另一个隐藏原因是t_svd.m的逆FFT结果产生轻微复数分量，后续取模运算没有做实部截断，最终污染了更新过程。

解决：第一步先做数据归一化，把观测张量的元素值缩放到[0, 1]区间，方法是Y_obs = (Y_obs - min(Y_obs(:))) / (max(Y_obs(:)) - min(Y_obs(:)))。第二步在t_svd返回前使用real()去除虚部残留。这两步处理完，迭代过程的数值稳定性会明显改善。

5.3 收敛缓慢或震荡：rho参数的调整没有章法

现象：迭代次数已经到500的上限，相对变化曲线却还在缓慢抖动，PSNR始终达不到合理水平。

原因：ADMM的惩罚参数rho直接影响对偶变量更新的步长。rho过小，对偶变量走得慢，收敛自然拖沓；rho过大，前期震荡明显，甚至直接越过最优点。这两个方向在曲线上的表现需要仔细区分。

解决：观察history.rel_change的曲线形态。如果曲线单调下降但是很平缓，增大rho到原来的5到10倍；如果曲线在开始阶段明显锯齿状振荡，减小rho到原来的1/5左右。调一次跑一次，一两次就能锁定合适的区间。经验数据是当rho在0.005到0.05之间时，大部分中等规模实验的收敛表现稳定。

5.4 恢复结果边缘模糊：稀疏噪声项权重lambda压得太狠

现象：PSNR数值不低，但肉眼观察恢复图像时，边缘和细节纹理被明显抹平，画面有一种过度平滑的塑料感。

原因：lambda设置的权重偏高，L1范数对稀疏项的压力过大，算法为了减少稀疏项的能量愿意牺牲掉一部分结构细节，结果就是噪声被去除的同时细节也一并被压缩了。

解决：把lambda缩小到原来的1/3左右，重新跑一遍，对比生成的PSNR和图像观感。需要理解的是，去噪任务始终存在一个“保结构”与“去噪声”之间的张力，lambda在中间某个区间能同时兼顾两者，而偏大偏小都会让某一方受损。严格评估时，我会用一组不同lambda值跑出曲线，找PSNR曲线的峰值点作为理论最优值。

6. 二次开发路线：把张量分解拓展到视频修复与多光谱领域

如果你已经跑通了上面的demo，二次开发的方向就完全打开了。从工程改造的角度，这里给出三条不同深度的扩展路线，你可以根据自己手头任务的类型选择。

第一条是更换数据输入，把合成数据换成真实场景数据。我在做视频去噪时把gen_noisy_data.m的随机张量替换成一段真实监控视频的灰度帧堆叠，具体方法是读取视频每一帧转成灰度图，再用cat(3, frames{:})堆叠成三维张量。这里需要注意三点：首先是灰度化之后像素值域通常还是0-255，需要先归一化到[0,1]再输入算法；其次是帧数不要一次性全部堆进去，我一般每30帧作为一个批处理窗口，这是因为张量核范数在时间维度上假设切片间存在周期或平稳特性，短窗口的恢复效果更稳定；最后是真实视频里的噪声通常不是简单的高斯白噪声，可能需要把噪声类型修改为脉冲噪声，对应调整E项的更新函数。

第二条是改进子问题的求解方式。update_X.m目前是用软阈值操作近似求解张量核范数最小化，如果你处理的数据有特殊的物理约束，可以替换为其他低秩逼近策略，比如奇异值硬阈值法或加权核范数最小化。硬阈值法只保留最大若干奇异值，其余置零，缺点是非凸；加权核范数则对不同奇异值赋予不同权重，能更好保持图像纹理。这条路线对正在做算法创新研究的开发者比较合适，因为修改的地方高度集中，实验设计和效果对比都很方便。

第三条是部署到更大规模的数据上。当前t_svd.m对n3个频域切片逐个执行SVD，当n3是一个很大的数字时，可以改成parfor并行循环，前提是Matlab安装了Parallel Computing Toolbox。我在一组256x256x200的张量上测试过，parfor大概能带来3到5倍的加速，同时要确保每个worker都能访问到所需的数据。需要注意的是t_svd的输出U、S、V在时域中需要保证各频点的对应顺序，parfor写索引时不要乱序。

更深一层的改造方向是把张量分解与机器学习结合。一种做法是调用张量分解得到的U因子矩阵作为低维特征，接上分类器做活动识别任务；另一种是把它嵌入到深度学习pipeline的前置模块，先通过ADMM做数据预清洗，再输入到神经网络训练。这两个方向目前在相关领域都比较活跃，大部分尝试都集中在如何更自然地衔接张量分解的输出和下游模块的输入维度。

作为一次完整的二次开发走查，我建议你完成以下动作：先把第3.2节的合成数据换成自己的真实数据；然后调整参数跑出一个baseline效果；再对比修改update_X和t_svd两个文件后PSNR和RSE的变化；最后写一个简单的系数实验结果表来论证你的改动是否真的有效。这套流程走下来，你对这个项目的理解就完整了。

在做张量分解项目初期，我为了省事常常直接套用默认参数跑实验，结果在某个高光谱数据上PSNR凭空低了两三个dB，排查才发现是数据量级没有归一化导致的。从那以后我每次拿到新的输入数据都会强制走一遍“归一化、掩码检查、参数网格搜索”三步流程，这套习惯到现在都还在用。希望这篇拆解能帮你绕过一些我踩过的坑，更顺畅地跑通并改进这个张量分解项目。

本文还有配套的精品资源，点击获取