第8章 算子分裂 II：原始-对偶 (PDHG/Chambolle–Pock)

本章摘要

如果说 ADMM 是变分图像处理的“瑞士军刀”（通用、稳定但有时笨重），那么 原始-对偶混合梯度法（Primal-Dual Hybrid Gradient, PDHG），特别是 Chambolle-Pock 算法，就是一把精密的“手术刀”。

当你面对以下情况时，PDHG 是不二之选：

  1. 算子 $A$ 极其巨大且无法快速求逆（如 3D CT、非均匀 MRI、大型模糊核）。
  2. 模型包含多个非光滑项（如 TV-L1 去模糊，数据项和正则项都不可导）。
  3. 需要极高的并行度（PDHG 的每一步都是简单的矩阵向量乘法，极其适合 GPU）。

本章将带你深入 PDHG 的数学核心、工程实现技巧（如预条件）以及如何通过“算子堆叠”技巧优雅地解决复杂组合问题。

8.1 为什么要超越 ADMM？

在第7章中，我们看到 ADMM 通过引入辅助变量将问题解耦。其核心代价在于 $x$-更新步通常需要求解如下线性方程： $$ (A^TA + \rho D^TD)x = \text{rhs} $$ 其中 $A$ 是成像算子，$D$ 是正则化算子（如梯度）。

8.1.1 矩阵求逆的诅咒

  • 好情况：如果 $A$ 是恒等算子（去噪）或循环卷积（去模糊），且边界条件是周期的，我们可以用 FFT 在 $O(N \log N)$ 内完成求逆。
  • 坏情况：
    • 非周期边界：FFT 失效。
    • 空间变化模糊：模糊核随像素位置变化，$A$ 没有循环结构。
    • 断层扫描 (CT/PET)：$A$ 是 Radon 变换，没有快速的伪逆。
    • 丢失像素 (Inpainting)：$A$ 是对角阵但含 0，与 $D^TD$ 结合后结构复杂，无法对角化。

8.1.2 PDHG 的“显式”优势

PDHG 算法属于 显式方法 (Explicit Methods)。它不需要求解线性方程组，只需要计算：

  1. 前向投影：$v = Kx$
  2. 反向投影（伴随）：$u = K^* y$

这使得 PDHG 的每一次迭代复杂度严格为 $O(N)$（对于稀疏算子），并且不仅避免了矩阵求逆，甚至避免了显式存储矩阵 $A$——你只需要写出实现 $Ax$ 和 $A^Ty$ 的函数即可。

8.2 鞍点问题与算法推导

8.2.1 标准形式

我们关注通用的线性约束凸优化问题： $$ \min_{x \in X} F(Kx) + G(x) $$

  • $K$: 线性算子（将原始域映射到对偶域）。
  • $G(x)$: 原始空间上的凸函数（通常是数据项或非负约束）。
  • $F(z)$: 对偶空间上的凸函数（通常是正则项，如 $L_1$ 范数）。

利用 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) $$ 这是一个零和博弈：

  • 原始变量 $x$ 试图最小化能量。
  • 对偶变量 $y$ 试图最大化能量（寻找最“严苛”的约束违反方向）。

8.2.2 Chambolle-Pock 算法 (First-Order Primal-Dual Algorithm)

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$)：

  1. 对偶上升 (Dual Ascent)： $$ y_{n+1} = \text{prox}_{\sigma F^*} (y_n + \sigma K \bar{x}_n) $$

    解读：先沿着 $K \bar{x}_n$ 的方向走一步（试图最大化 $\langle Kx, y \rangle$），然后用 $F^*$ 的近端算子把它拉回来。注意这里用到的是 $\bar{x}$ (预测值)。

  2. 原始下降 (Primal Descent)： $$ x_{n+1} = \text{prox}_{\tau G} (x_n - \tau K^* y_{n+1}) $$

    解读：利用最新的 $y_{n+1}$ 计算梯度 $K^* y_{n+1}$，沿负方向走一步（最小化交互项），然后用 $G$ 的近端算子处理。

  3. 外推/动量 (Extrapolation)： $$ \bar{x}_{n+1} = x_{n+1} + \theta (x_{n+1} - x_n) $$

    解读：这是算法的灵魂。简单的交替梯度下降（Arrow-Hurwicz 方法）通常不收敛。这一步相当于告诉对偶变量：“根据 $x$ 的运动趋势，我猜下一次 $x$ 会在这里，请基于这个预测来更新你的 $y$。”

8.3 核心工具：Moreau Identity 与近端计算

在步骤 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) $$

常见函数的对偶 Prox 实例

| 原始函数 $F(z)$ | 物理含义 | 共轭函数的 Prox: $\text{prox}_{\sigma F^*}(v)$ | 实现代码思路 (Python) |

原始函数 $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$		 (不做任何改变，常数项)

8.4 步长选择与预条件 (Preconditioning)

PDHG 的收敛性严格依赖于步长 $\tau$ (原始) 和 $\sigma$ (对偶) 的选择。

8.4.1 基本收敛条件

$$ \tau \sigma L^2 < 1 $$ 其中 $L = |K|_2$ 是算子 $K$ 的算子范数（最大奇异值）。

Rule-of-Thumb：

  1. 对于 TV 正则化（$K=\nabla$），$L^2 \approx 8$（在 2D 图像中）。
  2. 对于一般算子，使用 幂迭代法 (Power Method) 估算 $L$：
    • 随机初始化 $x$。
    • 迭代 $k=1..10$: $x \leftarrow K^*(Kx)$; $x \leftarrow x / |x|$。
    • 最后的瑞利商 $\sqrt{|K^*(Kx)|/|x|}$ 即为 $L$ 的估计。
  3. 安全设置：$\tau \sigma = 1 / (1.5 L^2)$。

8.4.2 预条件：让算法提速 10 倍

当 $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$）

  • 对偶步长矩阵 $\Sigma$：由 $K$ 的行和（Row sum）决定。
  • 原始步长矩阵 $T$：由 $K$ 的列和（Column sum）决定。

直觉：如果某个像素被很多测量值覆盖（$K$ 的列和很大），它的更新步长 $\tau$ 就应该小一点，以免震荡；反之亦然。

8.5 实战模型构建：算子堆叠技巧

如何用 PDHG 求解复杂模型？关键在于构造 $K, F, G$。

案例 1：TV-L1 去模糊 (Robust Deblurring)

模型： $$ \min_x |Ax - b|_1 + \lambda |\nabla x|_1 + \iota_{x \ge 0}(x) $$ 此模型对椒盐噪声和运动模糊具有鲁棒性。

拆解：

  • $G(x) = \iota_{x \ge 0}(x)$ （非负约束）。
  • $\text{prox}_{\tau G}(x) = \max(x, 0)$。
  • 难点在于有两个 $L_1$ 项涉及线性算子。我们将它们堆叠起来。

令 $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 $$ 迭代步骤：

  1. 对偶更新（分块并行）：

    • $u_1 = y_1 + \sigma A \bar{x}$

    • $y_1^{new} = \text{prox}_{\sigma F_1^*}(u_1)$ $\to$ 投影到 $L_\infty$ 球并平移（处理 $-b$）。 (具体：$y_1 = \text{clip}(u_1 - \sigma b, -1, 1)$)

    • $u_2 = y_2 + \sigma \nabla \bar{x}$

    • $y_2^{new} = \text{prox}_{\sigma F_2^*}(u_2)$ $\to$ 投影到 $\lambda$-球: $\text{clip}(u_2, -\lambda, \lambda)$。

  2. 原始更新：

    • 计算总梯度：$v = x - \tau (A^* y_1^{new} + \nabla^* y_2^{new})$
    • $x^{new} = \text{prox}_{\tau G}(v) = \max(v, 0)$

  3. 外推：

    • $\bar{x} = 2x^{new} - x$

通过这种堆叠算子 (Stacking Operators) 的方式，PDHG 可以轻松处理任意数量的正则项和数据项组合。

8.6 动态加速：如果 $G$ 或 $F$ 是强凸的

如果目标函数中包含强凸项（例如 $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 $$ 这种策略不需要预知最优步长，参数会自动适应，收敛速度显著提升。

8.7 本章小结

  1. PDHG vs ADMM：PDHG 避免了线性方程求解，只依赖矩阵向量乘法，适合算子复杂或大规模问题；ADMM 在子问题易解（如 FFT）时收敛步数通常更少。
  2. 算法三部曲：对偶上升 (Prox on Dual) $\to$ 原始下降 (Prox on Primal) $\to$ 外推 (Extrapolation)。
  3. Moreau Identity：是实现对偶更新的桥梁，切记公式 $\text{prox}_{f^*}(v) = v - \sigma \text{prox}_{f/\sigma}(v/\sigma)$。
  4. 工程秘籍：
    • 使用 预条件 处理非均匀算子。
    • 使用 算子堆叠 处理多项混合模型。
    • 必须通过 点积测试 验证 $K$ 和 $K^*$。

8.8 练习题

基础题

  1. 对偶 Prox 推导： 设 $F(z) = |z|_1$。证明其共轭 $F^*(y)$ 是集合 $C = \{y : |y|_\infty \le 1\}$ 的指示函数。并由此推导出 $y_{k+1}$ 的更新公式是一个简单的截断操作（Clipping）。

    Hint 由定义 $F^*(y) = \sup_z (\langle z, y \rangle - |z|_1)$。利用 Holder 不等式，当 $|y|_\infty > 1$ 时，上确界为无穷大。

  2. 点积测试 (Dot Product Test)： 编写伪代码或文字描述，如何验证你实现的 divergence 算子确实是 gradient 算子的伴随（负转置）。

    Hint 生成随机 $u$ (图像) 和随机 $v$ (向量场)。计算 $s1 = \langle \nabla u, v \rangle$ 和 $s2 = \langle u, -\text{div } v \rangle$。检查 $|s1 - s2|$ 是否小于 $10^{-10}$。

  3. 参数计算： 若 $K = \nabla$ (图像梯度)，$G(x) = \frac{\alpha}{2}|x-f|^2$。这是 Tikhonov 正则化模型。请写出 $\text{prox}_{\tau G}(v)$ 的解析解。

挑战题

  1. 去雨模型 (Deraining)： 考虑图像分解模型 $I = B + R$（背景 + 雨线）。 $$ \min_{B,R} |\nabla B|_1 + \lambda |R|_1 \quad \text{s.t.} \quad B+R=I, B \ge 0, R \ge 0 $$ 请构建算子 $K$，定义 $x, F, G$，并写出 PDHG 的迭代步骤。

    Hint 令原始变量 $x = [B; R]$。约束 $B+R=I$ 可以放入 $G$（作为指示函数），或者作为硬约束放入 $F$（通过引入拉格朗日乘子）。建议将约束 $Ax=b$ 的形式放入 $K$，对应 $F$ 为指示函数 $\delta_{\{I\}}$.

  2. 内存占用分析： 对于 $N \times N$ 的图像，TV 去噪的 PDHG 需要存储哪些变量？总共需要多少个 $N^2$ 的存储空间？与 Split Bregman 相比如何？

8.9 常见陷阱与错误 (Gotchas)

1. 伴随算子的边界陷阱

这是新手 99% 会遇到的 Bug。

  • 错误：正向梯度算子用了 np.diff (导致尺寸减小)，反向散度算子简单地用 padding 补回去，但没仔细推导边界项。
  • 现象：算法收敛到一坨噪声，或者边缘出现剧烈振铃。
  • 对策：必须通过 8.8 题中的点积测试。如果测试失败，不要进行任何算法迭代，先修好算子。

2. 步长 $\sigma$ 在 Moreau Identity 中的位置

  • 错误代码：y = clip(y_tilde, -lam, lam)。
  • 正确代码：y = clip(y_tilde, -lam * sigma, lam * sigma) 吗？
  • 分析：也不完全对。对于 $F(z) = \lambda |z|_1$，其共轭 $F^*$ 是 $\lambda$-盒约束的指示函数。$F^*(y)$ 要求 $|y|_\infty \le \lambda$。Prox 算子是投影到该集合。
  • 最终正确：$y_{k+1} = \text{clip}(y_k + \sigma K \bar{x}_k, -\lambda, \lambda)$。注意这里 $\sigma$ 不乘在 $\lambda$ 上，因为它是指示函数的投影半径，与步长无关！但如果是 $L_2$ 数据项，$\sigma$ 会出现在分母中。请务必根据公式 $\text{prox}_{\sigma F^*}(v) = v - \sigma \text{prox}_{F/\sigma}(v/\sigma)$ 仔细推导。

3. 忘记外推

  • 现象：目标函数值下降极其缓慢，或者在某个值附近震荡。
  • 原因：没有使用 $\bar{x}$ 更新 $y$，而是用了 $x$。
  • Check：代码中必须有一行 x_bar = x_new + theta * (x_new - x_old)。

4. 预条件矩阵没更新

  • 如果在处理不同图像尺寸时，硬编码了预条件矩阵 $T, \Sigma$，会导致完全错误的步长，引发发散。