第9章 Split Bregman 迭代：Bregman 距离与“罚项不够硬”的修复

学习目标：

  1. 深度理解几何直觉：掌握 Bregman 距离如何测量凸函数的“非线性程度”，以及它如何将约束问题转化为无约束序列。
  2. 掌握算法推导：从原始的 Bregman 迭代推导至 Split Bregman，理解“把误差加回去”（Add-back-the-noise）的深刻含义。
  3. 精通求解细节：学会处理全变分（TV）模型中的各向同性与各向异性收缩算子，以及利用 FFT 快速求解线性子问题。
  4. 理清算法谱系：彻底搞懂 Split Bregman 与增广拉格朗日（ALM）、ADMM 之间的等价性与符号对应。

9.1 引言：罚函数法的困境与突围

在图像处理的变分模型中，我们经常遇到如下形式的约束优化问题： $$ \min_u E(u) \quad \text{s.t.} \quad Au = f $$ 或者含有很难直接优化的正则项的问题（如 $\ell_1$ 或 TV）。为了求解，最直观的方法是罚函数法（Penalty Method），将其转化为无约束问题： $$ \min_u E(u) + \frac{\lambda}{2} |Au - f|_2^2 $$

9.1.1 进退维谷的 $\lambda$

这里存在一个著名的权衡困境：

  • $\lambda \to \infty$：约束 $Au=f$ 被精确满足，但目标函数变得极度病态（Condition Number 爆炸），导致常规梯度下降法收敛极慢，甚至数值溢出。
  • $\lambda$ 较小：问题良态，求解快，但约束满足得差。在去噪任务中，这意味着图像过度平滑，对比度（Contrast）丢失；在压缩感知中，这意味着无法精确重建信号。

9.1.2 破局思路

2008-2009 年间，Goldstein 和 Osher 提出的 Split Bregman 算法改变了游戏规则。它的核心哲学是：不要试图一次性用巨大的 $\lambda$ 强行满足约束，而是用一个适中的 $\lambda$ 多次求解，每次把上一次的“违反量”加回到数据中去。

这种方法不仅速度极快（通常 10-50 次迭代即可收敛），而且由于每一步的 $\lambda$ 保持常数且适中，数值稳定性极佳。

9.2 数学基石：Bregman 距离

要理解这个算法，不能只看公式，必须理解其背后的几何量——Bregman 距离。

9.2.1 定义与几何意义

设 $J: \Omega \to \mathbb{R}$ 是一个凸泛函（Convex Functional），且 $u, v \in \Omega$。$J$ 在 $u$ 点关于 $v$ 的 Bregman 距离 定义为： $$ D_J^p(u, v) = J(u) - J(v) - \langle p, u-v \rangle $$ 其中 $p \in \partial J(v)$ 是 $J$ 在 $v$ 点的次梯度（Subgradient）。如果 $J$ 是光滑的，则 $p = \nabla J(v)$。

ASCII 几何直觉

想象 $J(u)$ 是一个碗状的凸函数曲线。我们在 $v$ 点做一条切线（或切平面）。

        y
        ^
        |                Func J(x)
        |               /
        |             /   
        |           /     
   J(u) | - - - - /|         <--- 实际函数值
        |       /  |
        |     /    |  gap = D_J(u,v)
        |   /      |
  Linear| /        |         <--- 线性化近似 (一阶泰勒展开)
 Approx |/         |              L(u) = J(v) + <p, u-v>

 Approx |/         |              L(u) = J(v) + <p, u-v>
        |          |
        |----------|-------------------------> x

        v          u

核心解读：

  1. $D_J(u,v)$ 是什么？ 它是函数 $J(u)$ 的实际值与它在 $v$ 点的一阶泰勒展开值之间的差。
  2. 凸性测量：它衡量了函数“有多弯”。如果 $J$ 是线性的，Bregman 距离恒为 0；$J$ 越“凸”，距离越大。
  3. 特例： - 若 $J(u) = |u|^2$，则 $D_J(u,v) = |u-v|^2$（欧氏距离平方）。 - 若 $J(u) = u^T A u$，则 $D_J(u,v) = (u-v)^T A (u-v)$（马氏距离）。 - 若 $J(u) = \sum u_i \log u_i$（负熵），则 $D_J(u,v)$ 是 KL 散度。

9.2.2 为什么不是“真”距离？

  • 不对称性：通常 $D_J(u,v) \neq D_J(v,u)$。
  • 不满足三角不等式。 但在优化中，我们利用的是它的非负性（$D_J(u,v) \ge 0$）和它可以消除正则项中复杂非线性结构的特性。

9.3 核心推导：从 Bregman 迭代到 Split Bregman

9.3.1 原始 Bregman 迭代 (Bregman Iteration)

考虑去噪问题 $\min_u J(u) + \frac{\lambda}{2}|f-u|^2$。如果不加修改，结果会丢失对比度。Osher 提出求解如下序列： $$ u^{k+1} = \arg\min_u D_J^p(u, u^k) + \frac{\lambda}{2}|f - u|^2 $$ 这看起来比原问题更难，因为引入了 $D_J$。但经过惊人的代数化简（利用次梯度条件 $p$ 的递推关系），上述迭代等价于：

  1. 求解原问题（但输入数据变了）： $$ u^{k+1} = \arg\min_u J(u) + \frac{\lambda}{2} |f^k - u|^2 $$

  2. 更新数据（把残差加回去）： $$ f^{k+1} = f^k + (f - u^{k+1}) $$

    直觉 (Adding back the noise)： 假设 $f$ 是含噪图像。第一次迭代后，$u^1$ 变得平滑，丢失了纹理。 残差 $r = f - u^1$ 中包含了“噪声 + 丢失的纹理”。 下一步我们让算法去逼近 $f + r = f + (f-u^1)$。这就好比告诉算法：“既然你上次去噪去狠了，这次我把丢失的部分加倍补给你，请你务必保留下来。”

9.3.2 变量分裂 (The "Split")

原始 Bregman 迭代虽然修正了偏差，但第一步 $\min J(u) + \dots$ 仍然很难解（特别是当 $J(u)$ 是不可微的 $\ell_1$ 或 TV 时）。

Split Bregman 的神来之笔在于：结合 变量分裂 和 Bregman 迭代。

考虑广义问题（如压缩感知、去模糊）： $$ \min_u |u|_1 + \frac{\mu}{2} |Au - f|^2 $$ Step 1: 引入辅助变量 $d$ $$ \min_{u,d} |d|_1 + \frac{\mu}{2}|Au-f|^2 \quad \text{s.t.} \quad d = u $$ Step 2: 转化为无约束优化（应用 Bregman 迭代） 将约束 $d=u$ 视为罚项 $\frac{\lambda}{2}|d-u|^2$，并对其应用 Bregman 迭代策略（引入变量 $b$）： $$ (u^{k+1}, d^{k+1}) = \arg\min_{u,d} |d|_1 + \frac{\mu}{2}|Au-f|^2 + \frac{\lambda}{2} |d - u - b^k|^2 $$ $$ b^{k+1} = b^k + (u^{k+1} - d^{k+1}) $$ Step 3: 交替最小化 (Alternating Minimization) 由于 $u$ 和 $d$ 耦合在二次项里，我们可以固定一个解另一个：

  1. 子问题 $u$（可微，二次）： $$ u^{k+1} = \arg\min_u \frac{\mu}{2}|Au-f|^2 + \frac{\lambda}{2} |d^k - u - b^k|^2 $$

  2. 子问题 $d$（$\ell_1$ 范数，Prox）： $$ d^{k+1} = \arg\min_d |d|_1 + \frac{\lambda}{2} |d - u^{k+1} - b^k|^2 $$

  3. 更新 $b$： $$ b^{k+1} = b^k + (u^{k+1} - d^{k+1}) $$ 这就是著名的 Split Bregman 迭代三步曲。

9.4 经典实战：全变分（TV）去噪与去模糊

让我们把上述抽象框架应用到具体的 TV 模型中。 模型：$\min_u |\nabla u|_1 + \frac{\mu}{2}|Ku - f|^2$ 其中 $K$ 是模糊核（去噪时 $K=I$）。

9.4.1 分裂设置

引入 $d_x \approx \nabla_x u, d_y \approx \nabla_y u$。 约束为：$d_x = \nabla_x u, \quad d_y = \nabla_y u$。

9.4.2 详细迭代步骤

  1. 更新 $u$（求解线性方程组） 目标函数关于 $u$ 求导并令为 0，得到正规方程（Euler-Lagrange方程）： $$ (\mu K^T K + \lambda (\nabla_x^T \nabla_x + \nabla_y^T \nabla_y)) u = \mu K^T f + \lambda [\nabla_x^T (d_x^k - b_x^k) + \nabla_y^T (d_y^k - b_y^k)] $$ 注意到 $\nabla_x^T \nabla_x + \nabla_y^T \nabla_y$ 在离散下就是拉普拉斯算子 $-\Delta$。 $$ (\mu K^T K - \lambda \Delta) u = \text{RHS} $$

工程实现 (Fast Solver)： 如果是周期性边界条件，卷积算子 $K$ 和拉普拉斯算子 $\Delta$ 均可被 FFT 对角化。 $$ u = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(\text{RHS})}{\mu |\hat{K}|^2 + \lambda (|\hat{\nabla}_x|^2 + |\hat{\nabla}_y|^2)} \right) $$ 这使得 $u$ 的更新仅需几次 FFT，速度极快。

  1. 更新 $d$（广义软阈值 / Shrinkage） 我们需要解决： $$ \min_{d_x, d_y} |(d_x, d_y)|_1 + \frac{\lambda}{2} (|d_x - \nabla_x u - b_x|^2 + |d_y - \nabla_y u - b_y|^2) $$ 这里需要区分两种 TV 定义：

  • Case A: 各向异性 TV (Anisotropic TV) 定义为 $|\nabla_x u|_1 + |\nabla_y u|_1$。此时 $d_x, d_y$ 完全解耦。 $$ d_x = \text{shrink}(\nabla_x u + b_x, 1/\lambda) $$ $$ d_y = \text{shrink}(\nabla_y u + b_y, 1/\lambda) $$ 其中 $\text{shrink}(x, \gamma) = \frac{x}{|x|} \max(|x|-\gamma, 0)$。

  • Case B: 各向同性 TV (Isotropic TV) 定义为 $\sum \sqrt{(\nabla_x u)^2 + (\nabla_y u)^2}$。此时 $d_x, d_y$ 在空间上耦合。 令中间变量 $s_x = \nabla_x u + b_x, \quad s_y = \nabla_y u + b_y$。 幅度 $S = \sqrt{s_x^2 + s_y^2}$。 $$ d_x = \max(S - 1/\lambda, 0) \frac{s_x}{S} $$ $$ d_y = \max(S - 1/\lambda, 0) \frac{s_y}{S} $$ (注：程序中需处理 $S=0$ 的除零情况，通常设结果为 0)

  1. 更新 $b$ $$ b_x^{k+1} = b_x^k + (\nabla_x u^{k+1} - d_x^{k+1}) $$ $$ b_y^{k+1} = b_y^k + (\nabla_y u^{k+1} - d_y^{k+1}) $$

9.5 Split Bregman 与 ADMM 的等价性

许多初学者会对 Split Bregman 和 ADMM (Alternating Direction Method of Multipliers) 感到困惑。它们看起来非常相似，但名字不同。

结论：对于凸问题，Split Bregman 算法 完全等价于 ADMM 应用在特定的增广拉格朗日形式上。

对照表 (Rosetta Stone)

| 概念 | Split Bregman 表示 | ADMM 表示 | 备注 |

概念 Split Bregman 表示 ADMM 表示 备注
主变量 $u$ $x$ 待恢复图像
辅助变量 $d$ $z$ 梯度/稀疏系数
惩罚参数 $\lambda$ $\rho$ 控制收敛速度
对偶变量 $b$ $u_{dual}$ 或 $\eta/\rho$ 累积误差 / 乘子
更新逻辑 $b \leftarrow b + (u-d)$ $u_{dual} \leftarrow u_{dual} + (x-z)$ 实质是梯度上升更新乘子

区别在哪？

  • 历史渊源：ADMM 源于 70 年代的优化理论；Split Bregman 源于 2000 年代末的图像处理社区，从 Bregman 距离的角度独立推导出来。
  • 直觉侧重：ADMM 侧重于“对偶上升”和乘子法；Split Bregman 侧重于“误差回填”和“迭代正则化”。Split Bregman 的视角对于理解为什么该算法能恢复由于正则化导致的纹理损失（Contrast recovery）更加直观。

9.6 本章小结

  1. 核心思想：Split Bregman = 变量分裂 (把难解的 $\ell_1$/TV 拆出来) + Bregman 迭代 (通过累积误差 $b$ 修正线性惩罚带来的偏差)。
  2. 三步范式： - 线性求解：利用 FFT 或 Gauss-Seidel 快速求逆。 - Shrinkage：利用解析解处理非光滑正则项。 - Bregman 更新：简单的加法，保证约束最终被满足。
  3. 效率：它是解决 $\ell_1$ 和 TV 正则化问题最高效的一阶算法之一，特别适合大规模图像问题。
  4. 哲学：它解决了优化中的“软硬之争”——用“软”的惩罚（便于计算）实现了“硬”的约束效果（高精度）。

9.7 练习题

基础题

  1. Shrinkage 推导： 给定优化问题 $\min_x |x| + \frac{\lambda}{2}(x - c)^2$，请通过对 $x>0, x<0, x=0$ 分类讨论导数，证明其解为 $x^* = \text{sign}(c)\max(|c| - 1/\lambda, 0)$。

    Hint 当 $x \neq 0$ 时，求导得 $\text{sign}(x) + \lambda(x-c)=0$。观察 $x$ 和 $c$ 同号。

  2. 参数角色： 在去噪模型中，$\mu$ 是模型参数（决定去噪强度），$\lambda$ 是算法参数（分裂惩罚）。

  • 如果 $\mu$ 变大，最终图像会更清晰还是更模糊？
  • 如果 $\lambda$ 变大，算法的收敛速度通常会变快还是变慢？（假设在合理范围内）
    Hint $\mu$ 控制保真项，越大越接近原图（噪点多）。$\lambda$ 控制 $u$ 和 $d$ 的耦合紧密度，过大导致 $u$ 步长太小（Stiff），过小导致约束满足慢。
  1. Bregman 距离计算： 设 $J(x) = x^T Q x$，其中 $Q$ 是正定矩阵。计算 $D_J(u, v)$ 并化简。它是什么几何形状？
    Hint 利用 $\nabla J(v) = 2Qv$ 代入公式。结果应为 $(u-v)^T Q (u-v)$，即马氏距离。

挑战题

  1. 去卷积的边界效应： 在使用 FFT 求解 $u$ 子问题时，我们隐含了图像是周期延拓的。对于非周期性的真实图片，这会在边界产生明显的十字伪影。请提出两种减少这种伪影的工程方法。

    Hint 1. `edgetaper`：模糊算子对边缘做平滑衰减。2. 边界填充：将图像向外扩充（padding）后再做 FFT，解完后裁掉。

  2. CS-MRI 重建： 写出压缩感知 MRI 重建（Partial Fourier + TV）的 Split Bregman 迭代公式。 模型：$\min_u |\nabla u|_1 \quad \text{s.t.} \quad R \mathcal{F} u = y$。其中 $R$ 是采样掩码。

    Hint 这是约束型问题。除了 $d=\nabla u$，还需要处理 $R \mathcal{F} u = y$。可以直接把数据项看作硬约束，或者用罚函数 $\frac{\mu}{2}\|R\mathcal{F}u - y\|^2$。如果是后者，更新 $u$ 时矩阵求逆会变成 $( \mu \mathcal{F}^{-1} R^T R \mathcal{F} - \lambda \Delta )$，利用 $R^T R$ 是对角阵（掩码）的特性求解。

  3. 收敛性陷阱： 证明如果去掉 $b$ 的更新步骤（即 $b^{k+1}=0$ 恒成立），该算法将无法收敛到原问题 $\min |\nabla u|_1 + \frac{\mu}{2}|u-f|^2$ 的解，而是收敛到某个近似解。这个近似解对应什么优化问题？

    Hint 退化为标准的二次罚函数法，解的是 $\min \|\nabla u\|\_1 + \frac{\lambda}{2}\|\nabla u - \nabla u\|^2 + \text{data term}$（此处逻辑需修正为：实际上退化为交替最小化求解 $\min \|d\|\_1 + \frac{\lambda}{2}\|d-\nabla u\|^2 + \text{data}$，这实际上是在解一个平滑过的 TV 模型，且 $\lambda$ 有限时必有偏差）。

9.8 常见陷阱与错误 (Gotchas)

  1. 分母为零 (Division by Zero)： 在各向同性 TV 的 shrinkage 公式中，出现 $\frac{s_x}{S}$。当 $S=0$（即平坦区域）时，程序会崩溃或产生 NaN。 Fix：实现时使用 max(S, 1e-10) 或者 np.where 掩码处理。

  2. FFT 的中心化 (Shift)： MATLAB/Python 的 fft2 默认零频在角落，而高斯模糊核或差分算子通常定义在中心。 Fix：必须使用 ifftshift 将算子调整到与 fft2 兼容的格式（零频在 (0,0)），或者显式构造 psf2otf。算子相位不对会导致重建图像发生平移或产生振铃。

  3. $\lambda$ 选得太激进： 虽然 Split Bregman 对 $\lambda$ 鲁棒，但如果 $\lambda$ 设得极大（例如 $10^5$），子问题 $u$ 的更新几乎完全被 $\lambda \Delta$ 主导，忽略了数据项，导致收敛极其缓慢（Stiffness）。 Rule-of-thumb：对于归一化到 [0,1] 的图像，$\lambda$ 通常取在 $0.1 \mu$ 到 $10 \mu$ 之间。

  4. 各向同性 vs 各向异性混用： 如果代码里 $u$ 子问题用了拉普拉斯算子（隐含各项同性假设），而 $d$ 子问题用了简单的分开软阈值（各向异性），这会导致求解的模型方向偏差，出现细微的棋盘格伪影。

  5. 忘记更新 $b$： 这是新手最容易犯的错。如果忘了 $b$，算法就退化成了普通的交替最小化（Coordinate Descent），无法精确满足约束，导致图像看起来像被“过度磨皮”了。