本章摘要
如果说 ADMM 是变分图像处理的“瑞士军刀”（通用、稳定但有时笨重），那么 原始-对偶混合梯度法（Primal-Dual Hybrid Gradient, PDHG），特别是 Chambolle-Pock 算法，就是一把精密的“手术刀”。
当你面对以下情况时，PDHG 是不二之选：
- 算子 $A$ 极其巨大且无法快速求逆（如 3D CT、非均匀 MRI、大型模糊核）。
- 模型包含多个非光滑项（如 TV-L1 去模糊，数据项和正则项都不可导）。
- 需要极高的并行度（PDHG 的每一步都是简单的矩阵向量乘法，极其适合 GPU）。
本章将带你深入 PDHG 的数学核心、工程实现技巧（如预条件）以及如何通过“算子堆叠”技巧优雅地解决复杂组合问题。
在第7章中，我们看到 ADMM 通过引入辅助变量将问题解耦。其核心代价在于 $x$-更新步通常需要求解如下线性方程： \((A^TA + \rho D^TD)x = \text{rhs}\) 其中 $A$ 是成像算子，$D$ 是正则化算子（如梯度）。
PDHG 算法属于 显式方法 (Explicit Methods)。它不需要求解线性方程组，只需要计算：
这使得 PDHG 的每一次迭代复杂度严格为 $O(N)$（对于稀疏算子），并且不仅避免了矩阵求逆，甚至避免了显式存储矩阵 $A$——你只需要写出实现 $Ax$ 和 $A^Ty$ 的函数即可。
我们关注通用的线性约束凸优化问题： \(\min_{x \in X} F(Kx) + G(x)\)
利用 Fenchel 共轭 $F(z) = \max_y \langle z, y \rangle - F^*(y)$，我们将原问题转化为 原始-对偶鞍点问题 (Min-Max Problem)：
\[\min_{x} \max_{y} \mathcal{L}(x, y) = \langle Kx, y \rangle + G(x) - F^*(y)\]这是一个零和博弈：
Chambolle 和 Pock (2011) 提出的算法通过交替更新 $x$ 和 $y$，并加入关键的 外推 (Extrapolation) 步骤来打破震荡，确保收敛。
初始化：$x_0, y_0, \bar{x}_0 = x_0$。步长 $\tau, \sigma > 0$。松弛因子 $\theta \in [0, 1]$（通常取 1）。
迭代循环 ($n \to n+1$)：
解读：先沿着 $K \bar{x}_n$ 的方向走一步（试图最大化 $\langle Kx, y \rangle$），然后用 $F^*$ 的近端算子把它拉回来。注意这里用到的是 $\bar{x}$ (预测值)。
解读：利用最新的 $y_{n+1}$ 计算梯度 $K^* y_{n+1}$，沿负方向走一步（最小化交互项），然后用 $G$ 的近端算子处理。
解读：这是算法的灵魂。简单的交替梯度下降（Arrow-Hurwicz 方法）通常不收敛。这一步相当于告诉对偶变量：“根据 $x$ 的运动趋势，我猜下一次 $x$ 会在这里，请基于这个预测来更新你的 $y$。”
在步骤 1 中，我们需要计算 $\text{prox}_{\sigma F^}$。然而，通常我们只知道 $F$ 的性质（例如 $F(z)=|z|_1$），而不熟悉其共轭 $F^$。
Moreau 分解 (Moreau’s Identity) 解决了这个问题： \(x = \text{prox}_{\gamma f}(x) + \gamma \text{prox}_{f^*/\gamma}(x/\gamma)\)
将此公式变形，我们可以用 $F$ 的 Prox 来实现 $F^*$ 的 Prox：
重要公式： \(\text{prox}_{\sigma F^*}(v) = v - \sigma \cdot \text{prox}_{F/\sigma}(v/\sigma)\)
| 原始函数 $F(z)$ | 物理含义 | 共轭函数的 Prox: $\text{prox}_{\sigma F^*}(v)$ | 实现代码思路 (Python) |
|---|---|---|---|
| $L_1$ 范数 $\lambda |z|_1$ |
TV 正则, 稀疏编码 | $L_\infty$ 球投影 (Projection) $v \to \text{clip}(v, -\lambda, \lambda)$ |
np.clip(v, -lam, lam) |
| $L_2$ 范数 $\lambda |z|_2$ |
Group Lasso | $L_2$ 球投影 $v \to v \cdot \min(1, \lambda/|v|_2)$ |
v * min(1, lam / norm(v)) |
| 指示函数 $\delta_{C}(z)$ |
硬约束 (如 $Ax=b$) | 无约束 (恒等变换) $v \to v$ |
(不做任何改变，常数项) |
PDHG 的收敛性严格依赖于步长 $\tau$ (原始) 和 $\sigma$ (对偶) 的选择。
\(\tau \sigma L^2 < 1\) 其中 $L = |K|_2$ 是算子 $K$ 的算子范数（最大奇异值）。
Rule-of-Thumb：
- 对于 TV 正则化（$K=\nabla$），$L^2 \approx 8$（在 2D 图像中）。
- 对于一般算子，使用 幂迭代法 (Power Method) 估算 $L$：
- 随机初始化 $x$。
- 迭代 $k=1..10$: $x \leftarrow K^*(Kx)$; $x \leftarrow x / |x|$。
- 最后的瑞利商 $\sqrt{|K^*(Kx)|/|x|}$ 即为 $L$ 的估计。
- 安全设置：$\tau \sigma = 1 / (1.5 L^2)$。
当 $K$ 的各行或各列范数差异很大时（例如 CT 成像中穿过中心的射线与边缘射线的积分长度不同），全局标量步长 $\tau, \sigma$ 会受限于“最坏”的那个维度，导致整体收敛极慢。
对角预条件 (Diagonal Preconditioning) [Pock & Chambolle 2011]： 我们将标量步长替换为对角矩阵 $\Sigma$ 和 $T$。
\(\sigma_j = \frac{1}{\sum_i |K_{i,j}|^{2-\alpha}}, \quad \tau_i = \frac{1}{\sum_j |K_{i,j}|^{2-\alpha}}\) （通常取 $\alpha=1$）
直觉：如果某个像素被很多测量值覆盖（$K$ 的列和很大），它的更新步长 $\tau$ 就应该小一点，以免震荡；反之亦然。
如何用 PDHG 求解复杂模型？关键在于构造 $K, F, G$。
模型： \(\min_x \|Ax - b\|_1 + \lambda \|\nabla x\|_1 + \iota_{x \ge 0}(x)\) 此模型对椒盐噪声和运动模糊具有鲁棒性。
拆解：
令 $K = \begin{bmatrix} A \ \nabla \end{bmatrix}$，则 $Kx = \begin{bmatrix} Ax \ \nabla x \end{bmatrix}$。
对应地，对偶变量 $y$ 也分为两部分 $y = \begin{bmatrix} y_1 \ y_2 \end{bmatrix}$，其中 $y_1$ 对应数据项，$y_2$ 对应梯度项。
函数 $F$ 变为可分离求和形式： \(F(Kx) = F_1(Ax) + F_2(\nabla x) = \|Ax-b\|_1 + \lambda \|\nabla x\|_1\)
迭代步骤：
通过这种堆叠算子 (Stacking Operators) 的方式，PDHG 可以轻松处理任意数量的正则项和数据项组合。
如果目标函数中包含强凸项（例如 $G(x) = \frac{\mu}{2}|x|^2$），标准 PDHG 只能达到 $O(1/N)$ 收敛率。通过动态调整参数，我们可以达到 $O(1/N^2)$ 的加速收敛。
Alg 2 (加速版 PDHG): 若 $G$ 是 $\mu$-强凸的。 在每一步迭代中更新参数： \(\theta_n = \frac{1}{\sqrt{1 + 2\mu \tau_n}}, \quad \tau_{n+1} = \theta_n \tau_n, \quad \sigma_{n+1} = \sigma_n / \theta_n\) 这种策略不需要预知最优步长，参数会自动适应，收敛速度显著提升。
divergence 算子确实是 gradient 算子的伴随（负转置）。
这是新手 99% 会遇到的 Bug。
np.diff (导致尺寸减小)，反向散度算子简单地用 padding 补回去，但没仔细推导边界项。y = clip(y_tilde, -lam, lam)。y = clip(y_tilde, -lam * sigma, lam * sigma) 吗？x_bar = x_new + theta * (x_new - x_old)。