企业官网建设流程全解析

1. 为什么要把PDHG单独拎出来讲

这一讲是“最优化方法”系列里比较特别的一篇，主题叫PDHG，全称是Primal-Dual Hybrid Gradient，国内一般译成“原始-对偶混合梯度法”，但更多人干脆叫它Chambolle-Pock算法，因为2011年Chambolle和Pock那篇《A first-order primal-dual algorithm for convex problems with applications to imaging》把这套框架完整地定理化了。

在山东大学软件学院的最优化方法课上，我们一路从凸集、凸函数、次梯度、共轭函数，讲到近端算子、对偶理论、ADMM，到了第25讲终于轮到PDHG登场。这课安排到这时候，其实很讲究：PDHG几乎把所有前置概念串起来了，没有凸优化和对偶理论打底，直接看迭代格式会觉得“这到底在干嘛”。反过来，如果前面的知识都跟上了，PDHG的推导就像一块拼图刚好归位，你会发现它把“非光滑”“复合函数”“线性算子”这几个最麻烦的东西拆得干干净净。

PDHG能解的问题长这样：

min_x f(x) + g(Ax)

其中f和g都是凸函数，允许不光滑，A是线性算子。这个形式覆盖了大量工程和机器学习问题，Lasso回归、全变分图像去噪、压缩感知重构、最优传输、多任务学习，等等。所以这篇文章不只是讲一个算法，更想帮大家把“怎么把实际问题改写成这个形式”和“改完之后怎么确定每一步的更新公式”这两件事彻底搞明白。代码部分我用图像去噪做例子，因为每一步都直观，调试也方便，适合对照学习。

2. 原始-对偶形式：为什么非得绕这么一圈

2.1 直接算prox_{g∘A}基本是死路

如果你已经学过近端梯度法，看到min f(x) + g(Ax)这种形式，第一反应很可能是：用近端梯度做啊，光滑项走梯度，非光滑项走prox。问题在于近端梯度法的prox作用对象是g(Ax)这个整体，也就是要计算：

prox_{g∘A}(x) = argmin_z [ g(Az) + 1/(2t)||z - x||² ]

这个式子绝大多数情况下是算不出来的。为什么？因为A在函数内部跟g复合在一起，除非A是单位阵，否则近端算子里那个二次项和线性算子纠缠在一起，没有闭式解，硬要用内层迭代做的话，整个算法就变成双层循环，代价直接翻几倍。

我当初第一次碰这个问题，天真地选了Lasso的一个带测量矩阵的版本去试。Lasso的目标是min_x 1/2||y - Bx||² + λ||x||₁，其中B是一个随机高斯矩阵而不是单位阵。FISTA看着很美好，但迭代里需要prox_{λ||·||₁}(·)，这倒好算，可问题在于FISTA的梯度步是相对于整体目标函数的。整体目标里含1/2||y - Bx||²，这个函数对x的梯度是B^T(Bx - y)，需要计算矩阵乘法没问题，但如果你想在这个光滑部分上加速，需要估计B^T B的条件数，这个估计在大规模场景下非常费劲。

好，那换一种思路，把测量矩阵B吸收到线性算子A里，写成min_x f(x) + g(Ax)，其中f(x) = λ||x||₁，g(v) = 1/2||v - y||²，A = B。可惜g(Ax)里的A还是消不掉，prox_{g∘A}照样算不出来。

PDHG的关键洞察是：不跟这个g∘A死磕，转而把问题改写成鞍点形式，让A从prox里“退”出来，变成两个独立的prox，每个都有闭式解。这是整个算法的命门。

2.2 用Fenchel共轭把问题拉成min-max

对g做Fenchel共轭：

g(Ax) = sup_y [ ⟨Ax, y⟩ - g*(y) ]

其中g*是g的共轭函数。这一步看起来只是换了个写法，但它改变了问题的结构。原来的min_x f(x) + g(Ax)变成了：

min_x sup_y [ f(x) + ⟨Ax, y⟩ - g*(y) ]

把这一步看成一个博弈：x在最小化，y在最大化。在凸性条件满足时，min和sup可以交换顺序，解就是拉格朗日函数的鞍点。这就是为什么算法叫“原始-对偶”，因为迭代里同时更新原始变量x和对偶变量y。

拉格朗日函数写成：

L(x, y) = f(x) + ⟨Ax, y⟩ - g*(y)

这里定义域的问题先不细抠，你可以粗略认为，如果某个y让⟨Ax, y⟩ - g*(y)趋向正无穷，那么g(Ax)就是无穷，这个y不会成为最大值点。严谨的凸分析里需要检查相对内部条件（比如存在某个x让g(Ax)有限且连续），但在绝大多数实际问题上条件都满足。

到了这一步，PDHG的思路就很清晰了：对x做一次近端步来最小化f(x)（并同时带上一个线性项），对y做一次近端步来最大化⟨Ax, y⟩ - g*(y)（本质上是求g*的近端）。这样A就只在两个prox之间传递信息，不再藏在任何prox里面。

2.3 鞍点视角下的直观理解

鞍点问题相当于在玩一个双人零和博弈：x选手想让L尽量小，y选手想让L尽量大。两个选手轮流行动，各自朝自己的方向走一小步。如果每一步都走得合适，最终会收敛到“谁都不愿意先动”的点，也就是鞍点。

这就是PDHG的骨架。至于每一步“走”的步长怎么定、什么时候收敛、加速怎么做，都需要专门的收敛分析来回答。这也是为什么PDHG比单纯背迭代公式复杂——它不是一个固定的闭式算法，而是一族可以调节的方法，参数变了，行为和收敛速度都不同。

3. PDHG迭代格式：三步更新逐项拆解

3.1 算法完整形式

Chambolle-Pock算法（θ=1版本，也是最常用的）迭代过程如下：

给定初始x⁰, y⁰，取θ∈[0,1]，τ>0，σ>0，满足τσ‖A‖² < 1（后续细说），对k=0,1,2,...重复：

x^(k+1) = prox_{τf}(x^k - τ A^T y^k)

x̄^(k+1) = x^(k+1) + θ(x^(k+1) - x^k)

y^(k+1) = prox_{σg*}(y^k + σ A x̄^(k+1))

三行公式，第一行更新原始变量x，第三行更新对偶变量y，中间一行是外推步。很多人第一次看会问：为什么中间要插个x̄而不是直接用x^(k+1)？

这个外推步正是PDHG能加速的关键。直观上说，对偶更新时如果只用刚更新的x^(k+1)，相当于少看了一步“动量信息”；用x^(k+1)和x^k的差做外推，相当于给对偶更新一点“冲量”，类似加速梯度法里的外推。数学上，外推步能把收敛率从O(1/k)提到O(1/k²)吗？严格说，PDHG在一般凸问题上的收敛率是O(1/N)，不管有没有外推；但外推让常数因子大幅改善，实际收敛速度快得多。最重要的一点是，收敛性证明并不依赖θ的精确取值，只要0≤θ≤1就行，工程上直接取θ=1几乎总是最优选择。

3.2 原始更新：prox_{τf}里发生了什么

先看第一行：

x^(k+1) = prox_{τf}(x^k - τ A^T y^k)

如果不看prox，只看括号里的x^k - τ A^T y^k，这像不像梯度下降？对L关于x求梯度，如果不考虑f的不可微性，梯度的一部分就是A^T y。所以“x^k - τ A^T y^k”就是在沿着让拉格朗日函数减小的方向走一步。走了这一步后，还需要把f的不可微性吸收进来——prox_{τf}做的就是这件事。

prox_{τf}(v)的定义是：

prox_{τf}(v) = argmin_z [ f(z) + 1/(2τ)||z - v||² ]

它找一个点z，让f(z)别太大，同时z也别离v太远。1/(2τ)是平衡系数。τ越小，z越靠近v；τ越大，f的权重越大，z可以离v更远。从这个角度理解，prox可以看成是“考虑完f的惩罚之后，把点拉到一个新位置”。

对于很多常见f，prox有闭式解。比如：

• f(x) = 0（没有惩罚），prox_{τf}(v) = v，退化成普通梯度步骤
• f(x) = δ_C(x)（集合指示函数），prox_{τf}(v) = P_C(v)，即投影到集合C
• f(x) = λ||x||₁，prox_{τf}(v) = soft-thresholding算子：S_{τλ}(v_i) = sign(v_i)·max(|v_i| - τλ, 0)
• f(x) = (1/2)||x - b||²，prox_{τf}(v) = (v + τb)/(1+τ)

这些闭式解是PDHG实用性的基础。如果没有它们，每次迭代都要解一个内层优化，算法就失去意义了。

3.3 外推步：x̄的作用是什么

中间那一行x̄^(k+1) = x^(k+1) + θ(x^(k+1) - x^k)，其实就是把刚算出来的x^(k+1)沿着它的更新方向再推一段。

当θ=1时：

x̄^(k+1) = x^(k+1) + (x^(k+1) - x^k) = 2x^(k+1) - x^k

这是当前点相对于前一步的“镜像点”，在加速梯度法里很常见。当θ=0时，x̄=x，退化成无外推版本，收敛慢但在某些带约束问题里更稳。

为什么对偶更新要用x̄而不是x？核心在于算法的收敛证明。Chambolle-Pock在证明中构造了一个关于迭代点的势能函数，外推步能让势能函数沿着迭代单调下降得更快。如果去掉外推，算法仍然收敛，但界差一截；保持θ=1，可以以几乎零成本获得更快的实际收敛速度和更好的常数因子。

3.4 对偶更新：prox_{σg*}与Moreau分解

第三行：

y^(k+1) = prox_{σg*}(y^k + σ A x̄^(k+1))

括号里的y^k + σA x̄是对拉格朗日函数关于y做梯度上升步。因为L对y的梯度贡献是Ax，所以这一步沿着Ax方向走。然后prox_{σg*}处理g*的部分。

问题来了：g的近端算子好算吗？如果你已经会算g的近端，那么g的近端可以用Moreau分解得到，这是一个非常实用的恒等式：

prox_{σg*}(y) = y - σ·prox_{g/σ}(y/σ)

这个式子是什么意思？它告诉你：不用去求g的显式表达式，只要会算g的近端，就能算出g的近端。这里g/σ表示函数h(z) = g(z)/σ，也就是除以σ。这个恒等式的推导基于Fenchel共轭和Moreau包络的性质，我不展开证明了，但实际中极其好用。

举个例子。在图像去噪里经常g(v) = λ||v||₁，一次范数。它的共轭是：

g*(w) = δ_{[-λ, λ]}(w)

也就是在[-λ, λ]区间上的指示函数（如果是向量，就是无穷范数球的指示函数）。所以prox_{σg*}(y)就是把y投影到[-λ, λ]区间上，即clamp操作。

用Moreau分解也能得到同样的结论，但直接记住这些常见对偶近端算子的结果更方便。下表是常用函数及其对偶近端：

有了这些闭式公式，PDHG的每一步更新都只需要做加法和逐元素运算，没有内层迭代，非常快。

3.5 和ADMM的关系：别搞混了

很多同学学了ADMM之后再看PDHG，觉得两者很像，但其实关系要分几个层面说。

ADMM针对的问题形式是min f(x) + g(z), s.t. Ax + Bz = c，核心手段是引入辅助变量并构造增广拉格朗日，然后用交替方向依次更新x、z和对偶变量。ADMM的一大优点是稳定，但代价是如果A或B对应的子问题没有闭式解，需要内层迭代。PDHG针对的问题形式是min f(x) + g(Ax)，不需要显式约束，需要的是A的伴随算子A^T以及f和g的近端算子，这要求更低。

另一个实际区别：ADMM更新里有二次罚项，相当于每步解一个带强凸扰动的近端问题；PDHG没有增广项，而是靠步长τ、σ来控制发散风险。所以PDHG的实现更简洁，对内存的占用也更少，特别适合变量维度极高的成像问题。但也要强调：ADMM在很多带等式约束的问题上仍然有不可替代的优势，两者不是简单的谁替代谁。

4. 收敛条件与参数选择：τσ‖A‖² < 1背后是什么

4.1 收敛定理和证明思路

Chambolle-Pock的收敛定理大致是这样：如果f和g都是闭凸函数，且问题存在鞍点，取0≤θ≤1，步长τ和σ满足τσ‖A‖² < 1，那么序列(x^k, y^k)收敛到某个鞍点。更进一步，对任意下标N，有：

|f( x̄_N ) + g(A x̄_N) - (f(x*) + g(Ax*))| ≤ C / N

其中C是与初值和步长有关的常数，x̄_N是前N步迭代的加权平均。也就是说，在目标函数值上，收敛速度是O(1/N)。这个界对一般凸问题是最优的。

证明思路用到了“部分对偶”和“势能函数”的技巧。构造一个关于(x^k, y^k)的二次型势能函数：

E_k = (1/(2τ))||x^k - x*||² + (1/(2σ))||y^k - y*||²

然后证明每一步迭代都会让这个势能函数下降一定量，下降量正比于当前点与鞍点之间的原始对偶间隔。这里τσ‖A‖² < 1的作用就是保证某个中间矩阵是正定的，从而让“下降量”总是正的。如果取τσ‖A‖² ≥ 1，这个矩阵可能出现零特征值甚至负特征值，算法很可能震荡甚至发散。

4.2 谱范数‖A‖怎么算：从理论到工程

“‖A‖”在这里指A从L²空间到L²空间的算子范数，即‖A‖ = sup_{‖x‖≠0} ‖Ax‖ / ‖x‖。理论分析上，很多算子需要你自己推导这个界；工程实践中，至少有三种方法：

第一种，从算子定义出发推导。比如图像处理里的前向差分梯度算子∇，在二维图像上有经典结论‖∇‖² ≤ 8，也就是‖∇‖ ≤ 2√2。这个界来自离散梯度算子和Laplacian特征值的估计，用起来很方便。

第二种，用幂迭代法数值估计。把一个随机向量反复乘A^T A并归一化，收敛到的最大奇异值就是‖A‖²的估计。这个方法对任何线性算子都通用，代价是要多做几次矩阵向量乘，但PDHG每次迭代本来就要做A和A^T的乘法，所以在预处理阶段跑几十步幂迭代，成本完全可接受。

第三种，如果A是结构化的（如傅里叶变换、小波变换），可以用Parseval定理等工具算精确值，比如傅里叶变换的算子范数就是1。

以图像梯度算子为例说明怎么用更有意思。假设图像大小m×n，则梯度算子输出两个m×n的矩阵，其实可以把梯度算子看成一个2mn×(m+n)左右的矩阵，但没人真的存储它，只是在算法里用“差分+边界处理”来实现乘法。计算谱范数的严格上界是2√2，工程上更安全的做法是开头跑一步幂迭代，然后乘以0.99或0.95留个余量。

4.3 步长τ和σ的选法

理论上只需要τσ‖A‖² < 1，但τ和σ分别取多少直接影响收敛速度。最常见的选择有两种：

• 对称取法：τ = σ = 1/‖A‖，此时τσ‖A‖² = 1，刚好在边界上。为了安全，通常取τ = σ = 0.99/‖A‖。这个取法的好处是代码简单，收敛曲线比较平滑。
• 非对称取法：当f或g的光滑性层级不同时，可以一个取大一个取小。比如f是二次光滑函数，τ可以适当放大；g包含指示函数（约束），σ保持适中。

经验法则：如果迭代中出现振荡，特别是对偶变量y反复跳动，先不要动步长，检查是不是‖A‖估小了；步长直接乘以0.5看是否稳定。如果稳定了，说明问题出在谱范数估计上，而不是目标函数本身。

还有一些自适应的变体，比如预条件PDHG，它把固定步长换成随坐标变化的预条件矩阵，收敛速度能明显提升，但对初学来说，先把固定步长吃透更重要。

4.4 原始-对偶间隔：怎么判断收敛了

在实际代码里，不能去算真实误差，因为鞍点解是未知的。一个靠谱的替代指标是原始-对偶间隔（duality gap）：

Gap(x, y) = [f(x) + g(Ax)] - [ -f* (-A^T y) - g*(y) ]

注意到由弱对偶，f(x) + g(Ax) ≥ -f*(-A^T y) - g*(y)，所以Gap非负。当Gap趋近0时，说明当前(x,y)离鞍点非常近。这个指标好在它严格、可以计算，不需要知道真实解。但坏处是如果数值误差控制不好，Gap可能不会严格降到0。工程上常用做法是记录Gap的对数值，看到它下降缓慢了，就判定收敛。

另一个简单粗暴的判据是看x相邻迭代的变化量||x^(k+1) - x^k||是否小于某个绝对阈值。但在慢收敛场景下，这个判据可能误报早停，所以还是用Gap更严谨。

5. 实战：用PDHG做图像全变分去噪

5.1 问题建模

图像去噪的全变分（TV）模型是PDHG最经典的例子之一。观测图像u₀，含噪声，想恢复原图u，目标函数：

min_u 1/2||u - u₀||² + λ·TV(u)

TV(u)可以取各向异性形式TV(u) = ∑_i |(∇u)_i|，也可以取各向同性形式TV(u) = ∑_i √[ ((∇u)_i^x)² + ((∇u)_i^y)² ]。为方便推导，下面用各向同性TV，但它对应的对偶投影跟各向异性稍有不同，代码里我会写清楚。

把它映射到min f(x) + g(Ax)的框架，需要把定义里的分离变量设成：

• u作为原始变量，对应x
• f(u) = (1/2)||u - u₀||²
• A = ∇，即梯度算子
• g(v) = λ∑_i ||v_i||，其中v_i是每个像素处的梯度向量

这样g(Au) = λ∑_i ||(∇u)_i|| = λ·TV(u)，完美匹配。

5.2 对偶近端：各向同性TV的投影

看对偶变量的更新。g(v) = λ∑_i ||v_i||，其中每个v_i是二维向量。它的共轭是：

g*(p) = δ_E(p)，其中E = { p : ||p_i|| ≤ λ, ∀i }

这就是逐像素约束在每个位置上的对偶变量长度不超过λ。所以prox_{σg*}(p)就是对每个像素，把p_i投影到半径为λ的圆盘内：

p_i ← p_i · min(1, λ / ||p_i||)

如果选各向异性TV，则约束变成每个分量|p_i^x|≤λ且|p_i^y|≤λ，投影变成分别对每个分量做clamp，这就是两种TV在对偶形式层面的区别。代码实现时注意这一点，避免各向异性的模型写了各向同性的投影。

5.3 Python实现

我写一个完整的Python实现，用numpy做矩阵运算，用梯度算子的前向差分和它的负伴随（散度）算子。注意边界条件用Neumann边界，也就是图像边界外补零复制，这在成像问题里是标准操作。

import numpy as np import matplotlib.pyplot as plt def gradient(u): """前向差分梯度，Neumann边界条件 输入u: (H, W)灰度图 返回(gx, gy)，每个都是(H, W)形状 """ gx = np.zeros_like(u) gy = np.zeros_like(u) gx[:-1, :] = u[1:, :] - u[:-1, :] gy[:, :-1] = u[:, 1:] - u[:, :-1] return gx, gy def divergence(px, py): """梯度算子的负伴随，即散度算子 对应gradient的转置操作 """ dx = np.zeros_like(px) dy = np.zeros_like(py) dx[1:, :] = px[1:, :] - px[:-1, :] dx[0, :] = px[0, :] dy[:, 1:] = py[:, 1:] - py[:, :-1] dy[:, 0] = py[:, 0] return dx + dy def prox_f(u, u0, tau): """f(u) = 0.5*||u - u0||^2 的近端算子 闭式解: (u + tau*u0)/(1+tau) """ return (u + tau * u0) / (1.0 + tau) def project_ball(px, py, lam): """逐像素投影到半径为lam的L2球内 用于各向同性TV的prox_{sigma*g*} """ norm = np.sqrt(px**2 + py**2) # 避免除零：norm为0的位置不需要投影 mask = norm > lam scale = np.ones_like(norm) scale[mask] = lam / norm[mask] return px * scale, py * scale def pdhg_tv(u0, lam=0.1, tau=0.25, sigma=0.25, theta=1.0, max_iter=300): """PDHG求解TV去噪问题 参数: u0: 带噪声的观测图像 lam: TV正则系数 tau, sigma: 原始和对偶步长，需满足 tau*sigma*||gradient||^2 < 1 theta: 外推参数，通常取1 max_iter: 最大迭代次数 """ u = u0.copy() ubar = u.copy() px = np.zeros_like(u0) py = np.zeros_like(u0) gap_history = [] obj_history = [] for k in range(max_iter): # 1. 原始更新 gx, gy = gradient(ubar) # u = u - tau * (A^T p)，注意梯度算子的伴随是负散度 v = u - tau * divergence(px, py) u_new = prox_f(v, u0, tau) # 2. 外推步 ubar = u_new + theta * (u_new - u) # 3. 对偶更新 gx, gy = gradient(ubar) px = px + sigma * gx py = py + sigma * gy px, py = project_ball(px, py, lam) # 记录目标函数值，每10轮打印一次 if k % 10 == 0 or k == max_iter - 1: gx, gy = gradient(u) tv_val = lam * np.sum(np.sqrt(gx**2 + gy**2)) data_term = 0.5 * np.sum((u - u0)**2) obj = data_term + tv_val obj_history.append(obj) # 原始-对偶间隔（简化版） gap = np.abs(obj_history[-1] - (-np.sum(u0 * u) + 0.5 * np.sum(u0**2)) ) \ if False else None u = u_new return u, np.array(obj_history) # 生成测试数据：一个简单方块图加高斯噪声 np.random.seed(42) H, W = 128, 128 u_true = np.zeros((H, W)) u_true[40:90, 45:85] = 0.8 u_true[60:75, 65:75] = 0.3 noise = 0.15 * np.random.randn(H, W) u0 = np.clip(u_true + noise, 0, 1) u_denoised, hist = pdhg_tv(u0, lam=0.1, tau=0.25, sigma=0.25, max_iter=300) print("迭代完成，最终目标函数值:", hist[-1])

这段代码有几个点要解释。第一，divergence函数是gradient的负伴随，这一步是从“A^T y”项里来的。在梯度算子用前向差分定义时，它的伴随作用在东西上返回的正好是负的散度。如果这里符号搞反了，迭代直接发散，这是初学者最常踩的坑。第二，各向同性TV的对偶投影是逐像素的L2球投影，不是两个分量分开做clamp。第三，步长我取τ=σ=0.25，因为梯度算子谱范数上界是2√2≈2.83，所以τσ‖A‖² = 0.25² × 8 = 0.5 < 1，满足收敛条件。

5.4 运行效果和参数调优

在实际运行这个代码时，你会发现目标函数值在前几十轮快速下降，后面就慢慢趋平。如果λ取太大，图像会被磨得很平，细节全没了；λ取太小，噪声根本去不干净。我试了几组λ值的心得是：λ在0.05到0.15之间对高斯噪声比较合适；σ越大，去噪越激进，但可能过平滑；σ越小，保留细节越多但噪声残留严重。

一个更省心的做法是先用大λ跑一遍看整体轮廓，再逐步降λ微调细节。PDHG的好处是每次迭代代价很低，跑几百轮也就一两秒的事，所以可以多试几组参数再选合适的。想要更好的效果，还可以在目标函数里加一个数据保真项的权重，让不同区域自适应，但那是另一篇文章的内容了。

6. 常见问题与排查实录

6.1 迭代发散了，图像变成彩色噪点

这是PDHG最典型的失败模式。原因几乎都出在步长上：τσ‖A‖² ≥ 1。排查顺序如下：

第一步，确认收敛条件。如果梯度算子谱范数估计不准，去跑一个幂迭代，或者把当前τ和σ各乘以0.5试试。第二步，检查伴随算子实现。A和A^T必须满足<Au, v> = <u, A^T v>，用一个小随机向量验证，否则算法在数学上就是错的。第三步，检查prox是不是写错了，尤其是投影函数里的半径λ，如果投影半径写错，对偶更新方向就错了。

我在教学过程中发现，90%的“参数不好调”问题其实不是参数问题，而是代码里A^T实现错了。用一个小尺寸的随机数组，手算一次梯度+散度的数值，对比代码输出，五分钟就能定位。

6.2 收敛很慢，目标函数几乎不动

这种通常是参数比例失衡。τ和σ一个太大一个太小。当你发现原始变量x迭代像“蜗牛爬”但对偶变量y跳得飞快，大概率σ太大τ太小，反之亦然。尝试让τ和σ保持在同一数量级，或者用自适应预条件版本，比如对角线缩放预条件器，能显著加速。

另一个原因是目标函数本身病态。如果f很“平坦”，梯度信息弱，算法要绕很多圈才能走到鞍点。这时可以考虑把f改成强凸版本，比如加一个小的二次正则项，这在某些场景下反而更快。

6.3 原始-对偶间隔不下降

理论上Gap一定收敛到0，但如果你的问题没有鞍点，或者数值误差把Gap卡在某一个值附近，就需要检查问题是否严格可行。比如约束集合是空集，鞍点不存在，那算法会在某个边界来回震荡。更常见的是投影操作里数值裁剪太激进（比如把很小的数直接裁剪成0），导致Gap长期停滞在上界。解决办法是用更高精度的浮点类型，或者调整收敛判据，不再追求Gap到0，而是看相对变化率。

6.4 内存不足，因为A矩阵被显式存储了

PDHG在成像问题上最大的优势之一就是你永远不需要存储A矩阵，只需要实现A和A^T的乘法算法。如果初学者把A写成显式稀疏矩阵甚至稠密矩阵，到图像尺寸大一点就会爆内存。正确的做法永远是写成函数，在内部做差分、变换、采样等操作。这条建议适用所有算子分裂类算法。

6.5 常见问题速查表

7. 收尾：一些实际操作中的体会

PDHG用到现在，给我的感觉是“易上手，难精”。易上手是因为迭代格式只有三行，轮子搭起来非常快；难精是因为你对问题结构理解得越深，越能把步长、预条件、外推参数调得更漂亮。我自己踩过的坑里，印象最深的是伴随算子的一致性：花了整整一个晚上调试发散，最后发现是梯度算子边界条件的符号错了。所以强烈建议，写任何一个新问题前，先花五分钟用数值方法验证<Au, v>和<u, A^T v>相等，这一步能省下后面无数排查时间。

另外一个想特别提醒的点是，PDHG不是银弹。遇到强凸+高精度需求的问题，或者带复杂耦合约束的问题，ADMM和半光滑牛顿法仍然值得优先考虑。PDHG最好的应用场景是目标函数有结构、维度极高、精度要求中等的大规模问题，尤其是成像和信号处理里的经典正则化问题。

最后分享一个小技巧：如果你要在实际项目里反复跑PDHG，建议把谱范数估计、gap计算和参数调度都封装成通用模块，下次换问题只需要换prox函数和A算子。我在自己的工具库里就是按这个模式组织的，从图像去噪换到CT重建，只需要复用那套外骨架，改两个算子就完事。这也是PDHG这类算法在工程里最舒服的打开方式。