企业官网建设流程全解析

简介：面向高维数据去噪与张量分解的学习者和工程人员，这份Matlab开源代码包以凸优化为求解思路，通过张量低秩近似完成含噪数据恢复，可用于图像处理、视频分析、多模态数据挖掘等场景，适合具备Matlab基础、希望理解算法实现与实验复现的中高级读者。压缩包共48个文件，约55KB，其中45个为.m函数与脚本，另有1个Markdown说明文档及git配置类文件，体量轻便、目录集中，便于快速检索核心算法与辅助工具。已有739人学习。代码围绕张量/矩阵的ADM求解、软阈值处理、张量展开与重构、去噪及补全实验、结果绘图等环节组织，读者可据此梳理凸优化张量分解的实现脉络，调整迭代次数、阈值等参数开展自己的去噪实验，并复用其中的函数与实验脚本搭建验证流程。

1. 项目定位：一张被噪声糊住的图像，和一个低秩假设

去年帮做遥感方向的 A 同学看一批图，原图干净得能数清纹理，加完高斯噪声后细节全糊成一片，中值滤波一上去边缘也跟着软掉，调参数调得有点像是在赌运气。这份 MATLAB 工程处理的就是这类问题：它不靠卷积核滑动，而是把图像块堆成高阶张量，假设干净数据在某个展开维度上是低秩的，噪声则被当成稀疏扰动，再用凸优化的方式把两者拆开，低秩那部分留作去噪结果，稀疏那部分直接丢掉。工程是系统开源的，张量算子、求解器、测试图、评估脚本分开摆放，改算法不用动主干。适合两类人：要复现张量去噪基线做对比实验的研究生，和想在 MATLAB 里快速验证低秩分解想法的开发者。

2. 环境与运行：把这份 MATLAB 工程跑起来

2.1 目录结构与依赖清单

拿到压缩包先别着急双击 .m 文件，我一般是先把目录结构过一遍再决定从哪个脚本进。这类低秩张量去噪工程的目录划分有比较固定的习惯，核心是「算子层 / 求解器层 / 数据层 / 评估层」四块分开，主脚本只负责串流程：

tensor-denoise/ ├── main.m % 总入口，跑官方示例 ├── setup.m % 路径挂载，每个新会话先跑它 ├── core/ % 张量基础算子 │ ├── unfold.m fold.m │ ├── mode_n_product.m │ └── svt.m soft_threshold.m ├── solver/ % 凸优化求解器 │ ├── admm_rpca_tensor.m │ └── proximal_gradient.m ├── utils/ % 切块、归一化、评估 │ ├── im2col_patch.m col2im_patch.m │ └── psnr_local.m ├── data/ % 测试图 + 固定的噪声种子 └── results/ % 输出目录，首次运行自动创建

依赖上不复杂，MATLAB R2018b 之后的版本基本都能跑，主要用到 Image Processing Toolbox 里的imread/imwrite/imresize/rgb2gray这几个。如果你的 MATLAB 装的是精简版，没有图像工具箱，那读图这一步换成imread的替代实现也能凑合，但归一化和缩放会麻烦一点，建议还是补上。另外工程里如果引用了第三方的张量工具箱函数（常见的是提供ttm、hosvd、tucker_als这类名字的一组 .m 文件），把整个文件夹塞进core/或者单独建一个thirdparty/挂路径就行，不用装成 MATLAB 的正式 toolbox。

提示：MATLAB 的路径是全局的，同名函数冲突是这类工程最常见的问题。如果之前挂过别的张量工具箱，先用which -all ttm查一遍，出现两三个同名文件就先清路径。

2.2 三条命令跑通第一个去噪示例

我习惯把初始化单独写成一个脚本，而不是每个文件里都addpath一遍。原因很实在：一旦你把相对路径写死在子函数里，换个工作目录跑就全崩，这种问题排查起来非常费时间。

% setup.m —— 把工程目录挂到搜索路径，每个新会话跑一次 root = fileparts(mfilename('fullpath')); % 以当前文件为基准，规避 pwd 变化 addpath(genpath(fullfile(root, 'core'))); % 张量算子、近端算子 addpath(genpath(fullfile(root, 'solver'))); % ADMM / 近端梯度 addpath(genpath(fullfile(root, 'utils'))); % 切块、评估 addpath(genpath(fullfile(root, 'data'))); % 测试数据 if ~exist(fullfile(root, 'results'), 'dir') mkdir(fullfile(root, 'results')); % 先建目录，避免 imwrite 直接报错 end fprintf('工程根目录: %s\n', root);

fileparts(mfilename('fullpath'))这一句是关键，它拿到的是「setup.m 自己所在的目录」，不管你 cd 到哪都是对的。genpath会把子目录一起加进来，省得一层层写。如果你只想跑一次不污染全局路径，把addpath换成addpath(..., '-begin')，保证优先级压过那些同名函数。

接下来跑主示例。这类工程的主入口一般长这样，参数集中在opts结构体里，改参数不用动函数体：

% main.m —— 最小可运行示例：灰度图 + 高斯噪声 + 低秩稀疏去噪 setup(); % 先挂路径 img = im2double(imread(fullfile('data', 'test_gray.png'))); % 读图并归一化到 [0,1] rng(2024); % 固定种子，结果可复现 sigma = 25 / 255; % 噪声强度，按 0-255 尺度折算 noisy = img + sigma * randn(size(img)); % 加性高斯白噪声 noisy = min(max(noisy, 0), 1); % 截断，防止越界 opts.lambda = 0.06; % 稀疏项权重，越大去噪越狠，细节也越容易丢 opts.maxIter = 200; % 外层最大迭代次数 opts.tol = 1e-4; % 相邻两次迭代的相对残差阈值 opts.patch = 8; % 分块边长，8x8 是低秩分解的甜点区 opts.stride = 4; % 滑窗步长，越小越稳但耗时线性上升 [denoised, info] = tensor_denoise(noisy, opts); % 主入口，返回去噪图和迭代信息 fprintf('PSNR: %.2f dB -> %.2f dB\n', psnr_local(noisy, img), psnr_local(denoised, img)); imwrite(denoised, fullfile('results', 'denoised_gray.png'));

几个参数值得单独说清楚。sigma是噪声标准差，用 0-255 尺度除以 255 是因为图像已经im2double归一化到 [0,1]，不折算的话后面lambda的量级会完全对不上，去噪结果要么原样返回、要么整张糊成灰块。opts.lambda是低秩项和稀疏项的权衡系数，它不是「越大越好」的旋钮，本质上是你在「去噪强度」和「纹理保真」之间挑一个点，后面第 4 章会给一个扫参脚本。patch和stride决定重叠块的数量，8×8 配步长 4 意味着每个像素被覆盖 4 次，重叠平均能压掉块效应，代价是计算量翻四倍。

跑通之后results/里应该有一张图，命令行能看到两个 PSNR 数值，前者是含噪图的（大概 20 dB 出头），后者是去噪后的，正常情况下能涨 3～6 dB。如果第二个数反而更低，别急着怀疑算法，先看第 5 章的维度检查。

2.3 换自己的图：尺寸、通道与归一化

换成自己的图，有三件事必须先对齐，否则跑出来的结果看着像随机噪声，会浪费很多时间去查算法。

第一是数据类型。工程内部统一按double且落在 [0,1] 处理，如果你的图是uint8，一定要先im2double。直接double(img)会得到 0-255 的浮点，MSE 会大六个数量级，PSNR 直接变负数。

第二是尺寸。滑窗切块要求图像边长能被patch整除，或者工程里自己做了边界填充。常见做法是在im2col_patch里用padarray补边，补完再切、还原时裁掉。如果你的图是 1000×1000 这种不规整的尺寸，先用imresize缩到整数倍，或者在 opts 里加一个opts.pad = true的开关。

第三是通道。彩色图不要直接展开成三阶张量丢进求解器，除非你的求解器确实按 t-SVD 的管状纤维（tube fiber）在处理第三模。稳妥的做法是按通道循环，每个通道独立做一次低秩稀疏分解，最后合并：

% 彩色图逐通道处理，避免三个通道的尺度差异互相污染 img = im2double(imread(fullfile('data', 'test_color.png'))); [H, W, C] = size(img); out = zeros(H, W, C); for c = 1:C out(:,:,c) = tensor_denoise(img(:,:,c), opts); % 单通道复用同一个入口 end

逐通道的问题是会丢掉通道间的相关性，理论上不如直接做三阶分解，但胜在稳、好调、内存占用小。如果你要做对比实验，两种都跑一遍，通常彩色图上三阶分解能高出 0.3～0.8 dB，差距不算大。

3. 核心模块拆解：张量展开、低秩稀疏模型与两个近端算子

3.1 张量基础算子：unfold / fold / mode-n 乘积

整套算法建立在「张量可以沿任意一个模展开成矩阵」这个操作上。低秩假设在很多场景下不是对张量本身成立的，而是对某个展开矩阵成立，所以展开和还原的顺序一旦搞错，后面所有结果都是错的，而且是那种「不报错但结果很差」的错，非常难查。

function X = unfold(T, k) % 把 N 阶张量沿第 k 模展开成矩阵：第 k 模做行，其余模按原顺序排在列方向 sz = size(T); n = ndims(T); order = [k, setdiff(1:n, k)]; % 第 k 模提到最前 X = reshape(permute(T, order), sz(k), []); % 先换轴再拉平，两步不能合并 end function T = fold(X, k, sz) % unfold 的逆操作，sz 是原张量的尺寸向量 n = numel(sz); order = [k, setdiff(1:n, k)]; T = permute(reshape(X, sz(order)), [2:n, 1]); % 先按换轴后尺寸还原，再把轴换回去 end

permute(T, order)把第 k 模挪到第一维，reshape随后把这个维度保留为行数、其余所有元素压成列。写成一行reshape(T, sz(k), [])也能跑，但那是对内存里的原始顺序拉平，等价于沿第 1 模展开，和第 k 模完全是两回事，这是最典型的翻车点。fold里的[2:n, 1]是把第一维挪回第 k 位，逆置换就是原来的order循环右移一位。

有了展开，mode-n 乘积就是「把张量展开、乘矩阵、再折回去」：

function Y = mode_n_product(X, A, n) % 张量与矩阵沿第 n 模相乘：Y_(n) = A * X_(n) sz = size(X); Xn = unfold(X, n); Y = fold(A * Xn, n, [size(A,1), sz(setdiff(1:ndims(X), n))]); end

注意fold传进去的尺寸向量要把第 n 维换成size(A,1)，因为乘完之后那一维的长度变了。这一步如果还是传原尺寸，reshape会直接报元素数不匹配，属于运气好的情况；如果元素数恰好相等，就会静默地把数据错位，属于运气差的情况。

3.2 低秩加稀疏模型与 ADMM 求解器

去噪的核心模型可以写成一句话：观测张量 Y 等于低秩分量 L 加稀疏分量 S。干净图像在块展开后是强相关的，落在低秩子空间里；噪声在像素层面是独立扰动，落在稀疏项里。写成优化问题就是

min ||L||_* + λ||S||_1 s.t. Y = L + S

核范数（奇异值之和）是秩的凸松弛，L1 范数是稀疏度的凸松弛，两个都是凸的，加起来还是凸的，所以能用 ADMM 这种「分而治之」的框架拆成两个有闭式解的子问题轮流求解：

function [L, S, hist] = admm_rpca_tensor(Y, lambda, maxIter, tol) % 低秩 + 稀疏分解：min ||L||_* + lambda*||S||_1 s.t. Y = L + S [L, S, Z] = deal(zeros(size(Y))); % 原始变量 L、S，对偶变量 Z mu = 1e-3; % 增广拉格朗日步长，太大不收敛，太小极慢 hist = zeros(maxIter, 1); for it = 1:maxIter % L 步：对 (Y - S + Z/mu) 做奇异值阈值，天然压秩 L = svt(Y - S + Z / mu, 1 / mu); % S 步：对 (Y - L + Z/mu) 做软阈值，逐元素收缩 S = soft_threshold(Y - L + Z / mu, lambda / mu); % 对偶变量更新 + 相对残差监控 R = Y - L - S; Z = Z + mu * R; hist(it) = norm(R, 'fro') / (norm(Y, 'fro') + eps); if hist(it) < tol, break; end % 提前退出，省一半时间 end end

mu的选择直接决定收敛速度，工程里常见做法是从 1e-3 起步，每 50 次迭代乘 1.1 做一次自适应提升，或者干脆按mu = 1.25 / norm(Y(:), 1)估一个初值。Z/mu是对偶变量的缩放形式，写成Z除以步长是标准写法，不用改。hist存相对残差，出问题的时候把它画出来看一眼，正常情况下应该单调下降并在 50～100 次内落到 1e-4 以下，如果震荡或者干脆上升，基本就是mu或者lambda量级不对。

3.3 奇异值阈值与软阈值：两个近端算子

ADMM 的两个子问题都有闭式解，对应两个近端算子，代码短但每一步都有讲究：

function X = svt(A, tau) % 奇异值阈值：所有奇异值整体向 0 收缩 tau，小于 tau 的直接归零 [U, S, V] = svd(A, 'econ'); % econ 模式，只要前 min(m,n) 个奇异值，省内存 s = diag(S); s = max(s - tau, 0); % 截断在 0，得到的就是降秩后的矩阵 X = U * diag(s) * V'; end function X = soft_threshold(A, tau) % 软阈值：逐元素收缩，正负两侧都要处理 X = sign(A) .* max(abs(A) - tau, 0); end

svd一定要加'econ'。以 8×8 块拉平成 64×N 的矩阵为例，N 可能是几万个，全量 SVD 会算出大量零奇异值对应的左右奇异向量，白占内存还拖慢速度。max(s - tau, 0)这一步既做了收缩又做了降秩，tau越大留下的奇异值越少，低秩约束越强，去噪越狠。软阈值里的sign(A) .* max(abs(A) - tau, 0)是不能拆成max(A - tau, 0)的，后者会把所有负值直接抹平，图像会整体偏亮，这也是个不报错但结果明显不对的坑。

这两个算子单独看都很简单，但它们构成了整个算法的全部计算开销。真正跑大图的时候，SVD 占了 90% 以上的时间，所以任何优化都优先从这里下手：块尺寸调小、用随机化 SVD 近似前 k 个奇异值、或者干脆把svd换成svds(A, k)只算前 k 个。

4. 去噪主流程：从加载图像到 PSNR 评估的完整链路

4.1 主入口的执行链路

把前面几块拼起来，主函数负责的是「切块 → 逐块求解 → 重叠平均 → 输出」这条流水线：

function [Xout, info] = tensor_denoise(Y, opts) % 输入 Y: HxW 或 HxWxC 含噪图像；opts 见 main.m if ndims(Y) == 2, Y = reshape(Y, size(Y,1), size(Y,2), 1); end [H, W, C] = size(Y); p = opts.patch; s = opts.stride; Xout = zeros(H, W, C); info.iter = zeros(1, C); for c = 1:C % 1) 滑窗切块：每块拉成一列，堆成 (p*p) x Npatch [P, idx] = im2col_patch(Y(:,:,c), p, s); % 2) 低秩稀疏分解，只保留低秩分量 L，稀疏分量 S 丢掉 [L, ~, hist] = admm_rpca_tensor(P, opts.lambda, opts.maxIter, opts.tol); % 3) 重叠块加权平均还原，权重用双明窗压块效应 Xout(:,:,c) = col2im_patch(L, [H, W], p, s, idx); info.iter(c) = numel(hist); end end

im2col_patch返回的P是「每列一个块」的矩阵，idx记录每块在原图中的位置，col2im_patch靠它把块贴回去。重叠区域取平均是最朴素的做法，讲究一点就用双明窗（hann(p) * hann(p)'）做加权，边缘权重低、中心权重高，块效应能压掉一大半。info.iter记录每个通道的实际迭代次数，这个字段在扫参时很有用，能看出哪些参数组合压根没收敛就被maxIter截断了。

链路里唯一容易忽略的是col2im_patch的归一化。如果重叠区域直接累加不除权重和，输出图像整体会偏亮好几倍；如果除以固定的重叠次数而不是逐像素的权重和，边缘区域会偏暗。这两种错误都不会报错，只会在结果图上体现为一圈明显的亮边或暗边。

4.2 噪声注入与 PSNR / SSIM 评估

评估这件事值得单独写一段，因为很多人跑完只看一眼图，觉得「看着挺干净」就结束了，实际上 PSNR 可能只涨了 0.5 dB，属于参数没调对。评估函数自己写一个更可控：

function v = psnr_local(A, B, peak) % 逐像素 MSE 换算 PSNR，A、B 均为 [0,1] 的 double if nargin < 3, peak = 1; end mse = mean((A(:) - B(:)).^2); v = 10 * log10(peak^2 / max(mse, eps)); % eps 防除零，全黑图也能算 end

peak取 1 是因为数据归一化到了 [0,1]，如果你的数据是 0-255 的uint8，这里必须是 255，否则算出来的 PSNR 会偏高 48 dB 左右，看着漂亮但完全没意义。噪声注入那一步同理，sigma跟数据尺度必须匹配，25/255 对应的是 0-255 尺度下的 σ=25，这是图像去噪论文里最常用的设定，拿来做横向对比时直接用这个值就行。

想要更贴近主观观感，可以再算一个 SSIM，MATLAB 自带ssim函数（需要 Image Processing Toolbox），窗口用默认的 11×11 高斯窗即可。一般情况下 PSNR 涨 3 dB 对应 SSIM 涨 0.02～0.05，如果 PSNR 涨了但 SSIM 掉，说明去噪把结构抹平了，得把lambda往回压。

4.3 结果文件保存与批量扫参

单张图跑通之后，下一步基本都是扫参。手改opts.lambda跑十遍太慢，写个循环：

lambdas = [0.02 0.04 0.06 0.08 0.10 0.15]; best = struct('psnr', -inf, 'lambda', NaN, 'img', []); for lam = lambdas opts = struct('lambda', lam, 'maxIter', 200, 'tol', 1e-4, ... 'patch', 8, 'stride', 4); out = tensor_denoise(noisy, opts); p = psnr_local(out, img); fprintf('lambda=%.2f PSNR=%.2f dB\n', lam, p); if p > best.psnr best = struct('psnr', p, 'lambda', lam, 'img', out); % 记住当前最优 end end imwrite(best.img, fullfile('results', sprintf('best_lam_%.2f.png', best.lambda)));

扫参范围不用太宽，lambda在 0.02 到 0.15 之间基本能覆盖全部有意义的区间，再大就整张图糊掉了。真正需要联合扫的是patch和stride，块越大低秩性越强但局部细节丢得越多，常见组合是 6/3、8/4、16/8。我一般固定stride = patch/2，只扫patch，把搜索空间压到一维，跑起来快很多。

注意：扫参时一定要固定随机种子。噪声重新生成一遍，PSNR 会有 ±0.2 dB 的抖动，这个抖动足以让你在两个接近的参数之间做出错误判断。

5. 常见问题排查：几个容易出错的典型情况

5.1 维度与通道相关的问题

现象一：程序正常跑完，没有报错，但输出图像是错的，纹理错位成棋盘格或者斜条纹。原因基本锁定在unfold/fold的轴顺序不一致。很多人只改了一边，比如unfold里用了[k, setdiff(1:n,k)]，fold里却按[1:k-1, k+1:n, k]还原，两种排列在 k=1 或 k=n 时恰好等价，中间模就全错了，而且元素总数一致，reshape不会报错。 解决：写一个自检函数，随机生成 3×4×5 的张量，对每个 k 跑一遍fold(unfold(T,k), k, size(T))，用isequal对比是否还原成功。这个自检只要三行，放在测试脚本里跑一次，能省下几个小时的排查时间。

现象二：灰度图结果正常，彩色图整体偏色或者某个通道明显更糊。原因是三个通道的数值分布差异大，比如蓝色通道整体偏暗，同样的lambda对它来说就过强了。 解决：按通道做尺度归一化，每个通道单独减均值、除标准差，去噪完再反变换回去；或者干脆把lambda按通道的噪声估计值缩放，opts.lambda * sigma_c / sigma_mean。

5.2 收敛与参数相关的问题

现象三：迭代跑到maxIter还没停，hist曲线在前面几十次快速下降之后就贴着某个值不动了。这是典型的不收敛，mu太小导致对偶变量更新跟不上。也可能是lambda设得太小，稀疏项几乎不起作用，整个问题退化成一个近似恒等映射，残差自然降不下去。 解决：先把tol放宽到 1e-3 看能不能停，能停说明是精度要求过严；还是不停就把mu乘 10 再试。如果hist是震荡而不是缓慢下降，那是mu太大，反过来除以 3。工程里比较稳的做法是让mu自适应，每 20 次迭代比较一次残差，上升就减、下降就增。

现象四：去噪结果看着很干净，但 PSNR 只涨了不到 1 dB，边缘全部发虚。过平滑，lambda太大或者块尺寸太大。低秩约束本质上是「用少数几个基表示这个块」，块越大、秩越低，能表示的纹理就越少。 解决：先把patch从 16 降到 8 甚至 6，再把lambda减半，两件事一起做通常能找回 2 dB 左右。如果还是发虚，说明这个场景的低秩假设本身就站不住（比如强纹理的织物、密集的草地），换非局部相似块分组会更合适。

5.3 内存与路径相关的问题

现象五：小图能跑，换成 512×512 以上的图就报 Out of memory，或者 MATLAB 直接卡死。不是内存真不够，是中间变量炸了。滑窗切块后如果stride=1，一张 512×512 的图会产生 25 万块，每块拉成 64×1，堆起来是 64×250000 的矩阵，svd在这个尺寸上做全量分解会直接申请几个 G。 解决：把stride提到patch/2，块数直接降到四分之一；svd换成svds(A, k)只算前 k 个奇异值（k 取 5～15 就够）；再不行就按行分带处理，一次只处理几十行，结果写回大矩阵。这三招叠起来，普通笔记本跑 1024×1024 的灰度图没压力。

现象六：换台机器就报「未定义函数或变量 ttm」之类的错误。路径没挂全，genpath只挂了core/，第三方工具箱在别的目录。也可能是另一个同名函数在路径里排前面，MATLAB 找到了错的版本。 解决：which -all ttm查所有同名文件，restoredefaultpath之后重新跑setup.m，确认addpath的顺序里工程目录在最前。血泪经验是别把工程目录和别的项目放在同一个 MATLAB 工作区里来回切，路径残留是排查起来最没有线索的一类问题。

6. 二次开发：把求解器抽出来做模块替换

这份代码最值得改的地方不是主流程，是solver/那一层。主流程只是「切块—求解—还原」，跟具体用哪种分解没关系，真正的算法差异全在求解器里。我一般的改法是把求解器做成函数句柄注入，而不是在主函数里写if strcmp(method, 'tucker')这种开关：

% 用函数句柄把分解算法换掉，主流程一行不动 opts.rank = [20 20 3]; % Tucker 各模秩，第三模对应通道 opts.solver = @tucker_hooi; % 换成 t-SVD 就 @tsvd_solver opts.lambda = 0.06; [out, info] = tensor_denoise(noisy, opts); % 入口签名保持不变

改造时有一个参数表值得先背下来，它决定了你换算法之后要重新调哪些东西：

验证改造有没有效果，别只看最终 PSNR，那是黑匣子。我会在求解器里加一个返回值，把每次迭代的hist和低秩分量的秩一起吐出来，画成两条曲线。改造正确的标志是：残差单调下降且收敛迭代次数没有明显变多，低秩分量的有效秩从初始的高值稳定降到目标区间。如果秩一直降不下来，说明1/mu这个阈值太小，核范数约束根本没起作用；如果秩一步就掉到 1，说明阈值太大，整个块被压成了常数，去噪图会像被磨皮磨过。

另一个性价比很高的改法是把svd换成随机化 SVD。块矩阵动辄几千列，只取前十几个奇异值的场景下，随机投影 + 两次幂迭代就能逼近到 1e-3 的精度，速度能快三到五倍，PSNR 掉不到 0.1 dB。改的时候注意随机投影矩阵的列数取k + 10左右，取得太小（比如k + 2）精度会崩，这个参数没有理论最优值，得在你的数据上试。

从那以后我每次动求解器，都强制先跑一遍三项自检：unfold/fold的可逆性、阈值算子的符号对称性、以及固定种子下前后两次运行的 PSNR 是否完全一致。这三项里任何一项挂了，后面的调参全是白费功夫。希望帮到你。

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