简介：面向图像处理与凸优化入门读者的MATLAB实践包，围绕“图像上采样”这一典型问题，示范如何借助CVX把插值权重的选择建模为凸优化问题，并用PSNR衡量上采样图像与原始图像的接近程度。相比最近邻、双线性等常规插值，方案更强调以优化方式抑制伪影、提升重建质量，适合具备MATLAB基础、希望理解凸函数与优化建模如何落到图像任务的中级学习者。压缩包整体约22KB，文件总数显示为0、主要文件类型明细暂未提供，下载前可留意包内实际内容与运行环境。已有797人学习下载，说明该主题在课程实验与小规模复现中具有一定参考度，对课程设计、算法对比实验中的小规模验证与思路复现也较为友好。读者可据此梳理从读取图像、定义优化问题、调用CVX求解权重，到插入像素、计算PSNR并迭代调参的完整流程，也可将其作为高分辨率显示、图像分析与计算机视觉相关实验的入门样例。

1. 上采样不是插值：为什么 bicubic 之后再套一层凸优化能多出几个 dB

把一张 128×128 的图放大四倍，直觉上这是插值问题：bicubic 一个函数调用就完事。但插值隐含了一个很强的先验——图像局部可以用光滑曲面逼近。这个假设在平坦区成立，在边缘处必然过冲，放大后的图在轮廓两侧出现一圈振铃。换一条路：把上采样写成 y = DHx 的逆问题，在 matlab 里用凸优化（cvx 建模）求解，把先验从「光滑」换成「分段光滑」，同一张低分辨率输入，PSNR 常常能多出 1～3 dB，边缘也更干净。

关键不在调包，而在三个问题：D 和 H 从哪来、凸函数在这件事里到底保证了什么、cvx 的 DCP 规则为什么逼着你换一种写法。搞清楚这三件事，代码只是副产品。

2. 观测模型：y = DHx 里的三个算子分别从哪来

整条链路的起点不是优化，是建模。绝大多数上采样代码翻车的第一个地方，不是求解器不收敛，而是 A 这个矩阵压根就不对——模型错了，后面 PSNR 再高也是自欺欺人。

线性观测模型写成：

y = D H x + n

x 是高分辨率图拉直后的列向量（H·W 维），y 是低分辨率观测（m·n 维），H 是模糊矩阵，D 是下采样矩阵，n 是噪声。A = DH 就是「上采样算子」的逆——它描述的是高分辨率图怎么变成低分辨率图的，而我们要做的是解这个方程把 x 反推回来。

2.1 为什么模型里必须先有一个模糊核 H

如果只写 y = Dx，D 取隔行隔列抽样，那你等于在说：低分辨率图的每个像素恰好等于高分辨率图某个像素的值，其余全被扔掉。这在物理上不成立。真实成像里，传感器的像元有面积，光学系统有衍射极限，降采样之前必然先做了一次低通。忽略这个低通，模型就欠了最关键的一环。

后果是可见的：重建时，优化器会把混叠（aliasing）产生的虚假高频当成真实细节去拟合。低分辨率图里一条斜边上的锯齿，会被放大成摩尔纹或者一组平行条纹。这不是算法不行，是模型在骗它。

模糊核的取法有讲究。常用的实现思路是各向同性高斯，标准差 σ 与缩放因子 s 挂钩，一个常见的经验是 σ 取 0.5s 到 0.6s 之间；也有人直接用 s×s 的盒式滤波，等价于理想像元积分。选哪个会影响边缘的锐度：高斯衰减慢，重建结果偏软但稳定；盒式滤波锐利，但频谱有旁瓣，容易在强边缘附近引入振铃。我在小图上做对比时，一般先用盒式把链路跑通，再换高斯微调，因为盒式的物理含义更直白——像元积分，不会有人质疑参数怎么来的。

模糊矩阵本身是块 Toeplitz 结构，但没必要真的写成 Toeplitz。可分离核（高斯、盒式都可分离）用 Kronecker 积拼更快：

% 一维高斯核，长度 2r+1 t = -r:r; g = exp(-(t.^2)/(2*sigma^2)); g = g / sum(g); % 行方向模糊：作用于每一行的 W 个像素 Gh = kron(speye(H), toeplitz([g(ceil(end/2):end), zeros(1, W-numel(g))], ... [g(ceil(end/2):-1:1), zeros(1, W-numel(g))])); % 列方向同理，作用在 H 维 Gv = kron(toeplitz([g(ceil(end/2):end), zeros(1, H-numel(g))], ... [g(ceil(end/2):-1:1), zeros(1, H-numel(g))]), speye(W)); Hb = Gv * Gh; % H*W × H*W 的稀疏模糊矩阵

kron(speye(H), Gh_row)的意思是「把同一个一维模糊算子独立地作用到每一行」，kron(Gv_col, speye(W))是「作用到每一列」。两次可分离操作合起来等价于二维卷积，但非零元数量从 O(N·k²) 降到 O(N·k)，这是能不能在普通机器上跑起来的分水岭。

提示：模糊矩阵的边界处理必须和后面差分算子的边界处理保持一致。H 用了 zero padding（卷积核超出边界补零），差分算子却用 Neumann 边界，两者在图像最外一圈会打架，重建结果边缘出现一圈暗边或者亮边，PSNR 掉得莫名其妙。

2.2 用 Kronecker 积把二维差分算子拼出来

TV 正则需要图像的梯度，梯度是线性算子，也得写成矩阵。二维前向差分拆成水平和垂直两个方向：

function [Dh, Dv] = grad_ops(H, W) % 前向差分，保持方阵形式：末行/末列的梯度自然为 0（Neumann 边界） e = ones(W,1); D1 = spdiags([-e e], [0 1], W, W); % W×W，第 i 行是 x(i+1) - x(i) Dh = kron(speye(H), D1); % 沿行方向，逐行独立 e2 = ones(H,1); D2 = spdiags([-e2 e2], [0 1], H, H); Dv = kron(D2, speye(W)); % 沿列方向，逐列独立 end

这里有个容易被忽略的细节：spdiags([-e e], [0 1], W, W)造出来的是方阵，最后一行只有一个-1，因为x(W+1)不存在，超出的部分被丢掉了——等价于把边界梯度定义为 0。这正是 Neumann 边界条件。如果改成W-1 × W的瘦长矩阵，梯度向量维度就变成 (W-1)H，和后面的变量维度对不上，norms那一步会直接报维度不匹配。

方阵写法的代价是多了 H+W 个恒为零的约束，收益是Dh*x和Dv*x与x同维，可以直接拼成 2×N 的矩阵交给 CVX 的norms原子。这个取舍值得做。

Dh*x和Dv*x分别返回 H·W 维列向量，第 (i,j) 个位置的元素就是该像素的水平/垂直差分。把它们按行拼起来：

G = [Dh*x, Dv*x].'; % 2 × (H*W)，每一列是一个像素的梯度向量

这个 2×N 的矩阵就是 TV 项的原料。每一列是一个二维梯度向量 (∂x/∂h, ∂x/∂v)，对它取 2-范数得到梯度模长，再对所有像素求和，就是各向同性全变分。为什么是「各向同性」，原因就在这里——先求模长再求和，旋转不变；如果分开对两个方向求绝对值再相加，那是各向异性 TV，会把斜边缘惩罚得比水平/垂直边缘更重，结果斜边更容易被抹平。

2.3 下采样算子的两种写法与半像素错位

D 最直接的写法是挑点：

function D = down_op(H, W, s) % 每 s 行、每 s 列取一个点，等价于理想 decimation [r, c] = ndgrid(1:s:H, 1:s:W); lin = sub2ind([H W], r(:), c(:)); D = sparse(1:numel(lin), lin, 1, numel(lin), H*W); end

ndgrid生成所有采样点的行列坐标，sub2ind转成线性索引，sparse用一个 1 的对角把这些位置挑出来。这个矩阵每行只有一个 1，极度稀疏，构造开销可以忽略。

另一种写法是先模糊再抽样，把 D 和 H 合成一个矩阵 A = D*Hb。合成的好处是 CVX 里只出现一个线性映射，表达式更短；坏处是 A 的非零元是 Hb 每行的非零元数乘以采样率，核稍大一点内存就上去了。用分离形式D*(Hb*x)在 CVX 里同样合法（两个常数矩阵夹一个变量，是线性表达式），建模阶段不会真的把乘积展开成一个稠密矩阵。

真正会翻车的是栅格对齐。imresize(img, 1/s, 'bicubic')的像素中心映射是(i-0.5)/s + 0.5，也就是低分辨率第 i 个像素的中心对到高分辨率坐标(i-0.5)s + 0.5；而down_op取的是1:s:H，即高分辨率第 1、1+s、1+2s… 个像素。这两套栅格差半个像素。半个像素的平移对 PSNR 的杀伤力有多大？在纹理丰富的图上，足以让优化后的结果比 bicubic 还低 2～3 dB，而且肉眼看不出问题，只觉得「怎么糊了」。

自洽的做法是：要么低分辨率图也用down_op生成（做仿真实验时最稳），要么把 D 的采样起点改成ceil(s/2):s:H去对齐imresize的约定。两条路选一条走到底，别混用。

3. 凸函数在这里救了什么：从无穷多解到全局最优

模型搭好，接下来是解方程。A是 (m·n) × (H·W) 的矩阵，行数远少于列数——四倍上采样时，方程数量只有未知数的十六分之一。这是一个严重欠定的系统，解集是一整个仿射子空间，不是单点。

3.1 零空间与病态性：没有正则项时解为什么满天飞

只最小化数据项：

min ||A x - y||²

这个问题的解集是 x₀ + null(A)，null(A) 的维数至少是 H·W - m·n。也就是说，在满足「降采样之后和观测一致」的前提下，你可以往解里加任意多的高频成分而目标函数值不变。最小二乘会给一个最小范数解，但那个解在视觉上很糟——它把所有不在 A 的值域里的高频都设成 0，结果就是一张糊图，比 bicubic 还糊。

更麻烦的是数值层面。A 的条件数极大，因为模糊核把高频压到了接近 0，再做 decimation，很多奇异值小到 1e-10 量级。直接求伪逆，噪声会被放大成高频雪花。这就是典型的病态逆问题：不是没解，是解太多且都不稳定。

正则化就是在这时候进场的：

min ||A x - y||² + λ R(x)

R(x) 编码的是「什么样的图看起来像自然图像」这个先验。λ 控制先验和数据之间的相对权重。R 的选择决定了重建结果长什么样。

3.2 L2 与 TV 的凸性：几何直觉和一条能自己验证的思路

R(x) = ||x||²（Tikhonov）是最省事的，它让整个目标变成严格凸的二次型，有闭式解(AᵀA + λI)⁻¹Aᵀy，一次线性求解搞定。代价是它惩罚大梯度，等于在说「图像应该处处平坦」，边缘被均匀地糊掉。放大后的图看着像隔着一层毛玻璃。

R(x) = TV(x) = Σᵢ √((D_h x)ᵢ² + (D_v x)ᵢ²) 换了个说法：只惩罚梯度的模长，不惩罚梯度的大小。平坦区梯度接近 0，惩罚也接近 0；边缘处梯度很大，但因为是 L1 型惩罚（模长的一次方，不是平方），代价线性增长而不是平方增长。结果就是优化器宁可保留少量大梯度（边缘），也要把大量小梯度（噪声、纹理起伏）压到 0。这就是 TV 保边的机制来源。

两者都是凸函数，这一点可以自己验证，不用背结论：

• 2-范数||·||₂是凸的，这是三角不等式直接给的；
• 仿射映射复合保凸：D_h x是 x 的线性函数，所以||D_h x||₂关于 x 凸；
• 非负加权和保凸：Σᵢ 凸函数 = 凸函数。

所以 TV(x) 是「凸函数 ∘ 线性映射」的和，凸。数据项||Ax - y||²是凸二次函数，加一个 λ>0 的凸正则，整体还是凸。

凸性带来的不是「解更好看」，而是「解可预测」。非凸目标用梯度下降，换一个初值、换一个学习率、换一个 λ，落到的局部极小可能完全不同，调试过程变成玄学：这次跑出来边缘很锐，下次同样的参数跑出来糊成一团，你根本不知道是代码改了还是随机种子变了。凸问题不存在这个问题——最优性由 KKT 条件完全刻画，只要解满足 KKT，它就是全局最优，没有第二个。λ 连续变化时，解沿着一条连续路径移动，你可以放心地扫 λ 看趋势。

TV 有个小毛病：它凸但不严格凸。在整块平坦区域里，把这块区域整体加一个常数，TV 不变，数据项因为经过模糊和下采样，可能也几乎不变。所以解不唯一。实践中不影响，因为目标函数值唯一，PSNR 也就唯一，只是求解器可能在不同次运行给出略有差异的图。

3.3 λ 不是玄学：用偏差原则定量纲

λ 的取值经常被当成调参玄学，其实它有明确的下界。如果噪声是方差为 σ² 的高斯白噪声，那么在任何「合理」的解上，残差应该和噪声量级相当：

||A x - y||² ≈ (m·n)·σ²

这就是偏差原则（Morozov discrepancy principle）。λ 太小，残差远小于噪声水平，说明优化器在拟合噪声，图上有雪花；λ 太大，残差远大于噪声水平，说明数据项被忽略，解趋近于「纯 TV 最小化」的结果，图上出现阶梯状的色块。

实操上，把 y 归一化到 [0,1] 之后，λ 的合理区间通常在 1e-3 到 1e-1。先用 1e-2 跑一次，看残差：如果||Ax-y||²比(m·n)σ²小一个数量级，λ 加一个数量级；反之减。两三轮就能落到合理区间。比网格搜索快得多，而且每一步都有物理依据，不是碰运气。

注意：λ 的量纲依赖 A 的尺度。如果 A 的列没有归一化，或者 y 还是 0-255 的 uint8，λ 的绝对数值会差好几个数量级，抄别人的参数必然翻车。先把 y 归一化、把 A 的谱范数估出来，再谈 λ。

4. CVX 建模：DCP 规则下 TV 项为什么不能直接写 sqrt

模型清楚了，落到 cvx 上。第一次写 TV 的人几乎都会撞同一堵墙——把数学公式原样抄进代码，然后被 DCP 报错拦下来。

4.1 数据项与 sum_square 的锥转换

数据项||Ax - y||²在 CVX 里对应sum_square(A*x - y)。这是一个合法的凸原子，CVX 内部会把它转成二阶锥约束：

||Ax - y||² ≤ t ⟺ ||[2(Ax-y); t-1]||₂ ≤ t+1

内点法求解器（SDPT3、SeDuMi）处理的是标准锥形式，所以每一个平方项最终都会变成锥。这也是为什么 CVX 的规模上不去——锥的数量和变量维度同阶，而内点法每步迭代要做一次大规模线性方程组求解，复杂度随锥数量超线性增长。

sum_square和square_pos(sum(...))不是一回事。前者是「先逐元素平方再求和」，后者是「先求和再平方」，梯度结构完全不同，别混。数据项要的是前者。

4.2 用 norms 原子写逐像素梯度模长

数学上 TV 是sum(sqrt(dh.^2 + dv.^2))。直接写：

tv = sum(sqrt((Dh*x).^2 + (Dv*x).^2)); % 会被 CVX 拒绝

报错大致是Illegal operation: sqrt( {convex} )。原因在 DCP 规则：sqrt在 CVX 里被登记为凹函数（concave），规则是「凹函数只能出现在 ≥ 约束的左边，或者 maximize 的目标里」。这里 sqrt 的自变量是凸的，凸函数套凹函数，CVX 无法判定整体凸性，直接拒绝。

绕过的方法是用 CVX 专门为「范数之和」准备的原子norms：

cvx_begin quiet cvx_solver sdpt3 variable xh(H, W) x = xh(:); G = [Dh*x, Dv*x].'; % 2 × (H*W)，每列一个梯度向量 tv = norms(G, 2, 1); % 逐列 2-范数，再求和 minimize( 0.5*sum_square(A*x - y) + lambda*tv ) subject to xh >= 0; xh <= 1; cvx_end

norms(G, 2, 1)的语义是：沿第一维把 G 切成一个个列向量，分别求 2-范数，再把所有范数加起来。G 是 2×N，沿第一维切，切出来的就是每个像素的二维梯度向量，正好对上 TV 的定义。CVX 内部把它一次性转成 N 个三维二阶锥，建模阶段的符号处理是向量化的，速度可以接受。

注意：norms的维度参数在不同 CVX 版本上的行为不完全一致，有的版本dim=1表示按列分组，有的表示按行。第一次跑如果报维度错误，或者解出来的图完全没有 TV 效果（梯度没被惩罚），把1换成2试一次，两分钟能确认。

4.3 手写二阶锥约束：循环写法为什么在 CVX 里会卡死

想绕开norms的话，可以引入辅助变量 t 手动写锥约束：

cvx_begin quiet variable xh(H, W) variable t(H, W) x = xh(:); minimize( 0.5*sum_square(A*x - y) + lambda*sum(t(:)) ) subject to for i = 1:H*W norm([Dh(i,:)*x; Dv(i,:)*x]) <= t(i); end cvx_end

数学上没错，norm(·) ≤ t就是标准的二阶锥约束。但这段代码在 128×128 的图上就能让你等到怀疑人生：H·W = 16384 次循环，每次循环 CVX 都要做一遍表达式树的 DCP 检查、注册一个新的锥约束、维护约束列表。CVX 是解释执行的建模层，不是编译型 DSL，单次约束注册的开销在毫秒量级，一万六千次就是几十秒到几分钟，而且内存占用随约束数量线性上涨。256×256 时是 65536 次，基本等于挂死。

向量化写法和循环写法在数学上等价，在 CVX 里的建模时间差两三个数量级。这是 CVX 使用中最值得记住的一条经验：能用原子就别用循环，原子内部的锥转换是 C 层实现的。

顺带说一句，上面这段循环写法还有一个隐藏 bug：Dh(i,:)*x里Dh是稀疏矩阵，逐行取出来做矩阵向量乘，在 MATLAB 里每次都要重新索引稀疏结构，比稠密矩阵慢得多。真要写循环，先把Dh、Dv转成 dense，但那样内存又爆了。两头堵。

4.4 求解器与精度：SDPT3、SeDuMi 和 cvx_precision

CVX 默认调 SDPT3。对纯 SOCP（没有半定锥）的问题，SeDuMi 往往内存更省，尤其是在锥数量很大的时候。切换只需要一行：

cvx_begin quiet cvx_solver sedumi cvx_precision medium ... cvx_end

cvx_precision medium把内点法的终止容差从默认的 ~1e-8 放宽到 ~1e-4。对图像重建，1e-4 的精度早就够了，肉眼分辨不出来，但迭代次数能少三到四成。这个开关在调参阶段特别值——你要跑十几个 λ 看趋势，每个都高精度是浪费。

如果求解器报Failed或者Infeasible，先别怀疑模型。最常见的原因是约束写冲突，比如同时写了xh >= 0和某个隐含要求负值的项；其次是数据里有 NaN（imresize在极端缩放下可能产生），CVX 遇到 NaN 会静默失败，报一个和 NaN 毫无关系的错。养成习惯：进 CVX 之前先any(isnan(y))检查一次。

5. 避坑与排查：PSNR 算错、DCP 报错和内存爆掉时先看哪里

写这类代码，时间基本不花在「想出方法」上，花在「为什么跑出来不对」上。下面几条是踩得最深的。

5.1 PSNR 的两个隐形错误：数据范围和通道

PSNR 的定义简单：

function p = psnr_metric(ref, est, peak) % ref, est: 同尺寸矩阵，double 类型，取值范围 [0, peak] if nargin < 3, peak = 1; end d = double(ref) - double(est); mse = mean(d(:).^2); p = 10*log10(peak^2 / mse); end

两个参数能让你算出的数值偏离真值几十 dB，而且看不出错。

第一个是peak。数据归一化到 [0,1] 时 peak=1；数据是 0-255 时 peak=255。搞混了会差20*log10(255) ≈ 48 dB。更隐蔽的是：重建结果归一化到 [0,1]，参考图还是 uint8 的 0-255，两者相减得到的 mse 巨大，PSNR 变成负数——你会以为算法彻底失败，其实是单位没统一。

第二个是通道。彩色图直接对 R、G、B 三个通道求平均 PSNR，数值会虚高 1～2 dB，因为色度通道经过下采样后本来就平滑，重建误差小，拉高了平均值。学界惯例是转 YCbCr 只算 Y 通道：

ycc_ref = rgb2ycbcr(ref); ycc_est = rgb2ycbcr(est); p = psnr_metric(ycc_ref(:,:,1), ycc_est(:,:,1), 255);

还有两个细节：上采样后的图外围有 s-1 圈像素是靠边界外推得到的，误差偏大，算 PSNR 前裁掉；两张图如果有整像素平移，先做位移搜索对齐再算，否则纹理区错位会让 PSNR 平白掉 2～3 dB。

5.2 五条踩坑记录

一、Cannot assign values to a CVX variable

现象：想用 bicubic 的结果当初值，xh = imresize(y, s, 'bicubic')写在variable xh(H,W)后面，直接报错。

原因：CVX 的变量是符号对象，不是数值容器，声明之后不允许赋值。CVX 也不支持热启动，内点法本来也不需要初值。

解决：删掉赋值语句。如果确实需要初值加速（比如跑 FISTA），那就不该用 CVX，换一阶方法。

二、Illegal operation: sqrt( {convex} )

现象：把 TV 的数学定义原样抄进minimize。

原因：sqrt是凹原子，套在凸表达式外面违反 DCP。

解决：用norms(G,2,1)；或者引入辅助变量 t，写norm([dh;dv]) <= t，目标里用sum(t(:))。后者在变量维度小的时候更直观。

三、解出来整张图偏灰或者几乎全黑

现象：PSNR 极低，图看着像被压平了。

原因：三种可能。y 是 uint8 没转 double，CVX 隐式转换时丢精度甚至报错；λ 大了两三个数量级，数据项被彻底压制；变量没加 [0,1] 约束，解跑到范围外，显示时被 clamp 成一片。

解决：进 CVX 前统一y = double(y)/255；用 3.3 节的残差判据定 λ；把xh >= 0; xh <= 1;加上——这两条不只是为了显示，它们是真的在缩小可行域，能显著加快收敛。

四、重建 PSNR 比 bicubic 还低

现象：优化跑完了，一算比插值还差。

原因：栅格错位（2.3 节）或者模糊核 σ 与缩放因子不匹配。σ 取太小，模型以为低分辨率图里保留了更多高频，重建时过度锐化；σ 取太大，模型以为高频全丢了，重建结果偏软。

解决：仿真实验里用自己构造的 D、H 生成低分辨率图，保证模型自洽，先确认代码链路对，再换真实降采样图；换真实图时扫一遍 σ ∈ {0.4s, 0.5s, 0.6s, 0.7s}，看 PSNR 曲线的峰值在哪。

五、求解几十分钟不收敛，或者报内存不足

现象：128×128 还能跑，256×256 卡死。

原因：变量维度 H·W 和锥约束数量同阶，内点法的每次迭代要做一次大规模稀疏分解，内存随问题规模超线性增长。256×256 是 65536 维变量、65536 个三维锥，SDPT3 吃掉几 GB 很正常。

解决：按优先级——先开cvx_precision medium；换cvx_solver sedumi；把图分成 64×64 的有重叠块分别解再融合（重叠区加权平均，权重用余弦窗）；都不行就换 ADMM。CVX 本来就不是给大图准备的。

6. 设计取舍：CVX 该停在多大规模，什么时候必须换 ADMM

把 CVX 用在这个问题上，最需要想清楚的是它的定位。它不是生产工具，是一把标尺。

内点法的迭代次数对精度是对数级的，看起来很美，但每次迭代的代价是求解一个规模随锥数量增长的线性系统。CVX 在建模阶段把问题转成标准锥形式，这个转换本身是 O(N) 的符号操作，真正吃资源的是求解器内部的稀疏 Cholesky 分解。经验边界很清楚：128×128 以内，几十秒到几分钟能出结果；256×256 开始进入「泡杯茶回来看看」的区间；再往上，内存和时间都不划算。

替代方案是 ADMM。它把问题拆成两块：

min f(x) + g(z) s.t. Dx = z

x 子问题是最小二乘，在循环边界下AᵀA + ρDᵀD可以被 FFT 对角化，一次正反变换加逐点除法就解完，复杂度 O(N log N)；z 子问题是 TV 的近端算子，逐像素做二维收缩（把梯度向量整体往原点缩，模长小于阈值就归零），闭式解。每次迭代都是 O(N log N)，没有矩阵分解，内存只有几个 N 维向量。

表里的时间只是量级，具体数字取决于机器和核大小，但相对关系是稳的：同一张 256×256 的图，ADMM 跑几百次迭代的总时间，可能还不到 CVX 建模阶段的时间。

那为什么还要用 CVX？因为它是后悔药。手写 ADMM 的时候，x 子问题的边界条件、z 子问题的收缩阈值、ρ 的取值、残差判据，任何一处写错，结果都不会报错，只会悄悄变糊或者出条纹。你没法判断是「算法收敛到最优了但 λ 不对」还是「算法压根没收敛」。这时候在小图上用 CVX 跑一遍，拿到目标函数的最优值和重建图，就有了一个绝对参照：ADMM 跑出来的目标值如果比 CVX 高 0.5%，那就是没收敛；如果一样但图不同，那就是解不唯一（TV 的老毛病），不用管。

这个做法我在好几个项目里都用过。有一次写一个带方向梯度的各向异性 TV，ADMM 跑出来边缘总有锯齿，调了两天参数没结果。后来在 64×64 的小图上用 CVX 跑同样的模型，发现 CVX 的结果也带锯齿——问题不在算法，在模型本身，各向异性 TV 对斜边就是不友好。如果没有那个基准，可能还要再调两天。

从那以后我每加一个新正则项，都先用 CVX 在 64×64 的小图上跑一遍当基准，确认目标函数和解的形状都对，再换 ADMM 上大图。多花的这十分钟，省下的是后面几小时的自我怀疑。资源里的脚本按「构造算子 → CVX 求解 → PSNR 评估」拆成了几段，建议先单独跑算子构造那一段，把 A 的维度和稀疏度打出来看一眼，再往下走——很多问题在这一步就能暴露。希望帮到你。

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