学习目标：
- 深度理解几何直觉：掌握 Bregman 距离如何测量凸函数的“非线性程度”，以及它如何将约束问题转化为无约束序列。
- 掌握算法推导：从原始的 Bregman 迭代推导至 Split Bregman，理解“把误差加回去”（Add-back-the-noise）的深刻含义。
- 精通求解细节：学会处理全变分（TV）模型中的各向同性与各向异性收缩算子，以及利用 FFT 快速求解线性子问题。
- 理清算法谱系：彻底搞懂 Split Bregman 与增广拉格朗日（ALM）、ADMM 之间的等价性与符号对应。
在图像处理的变分模型中，我们经常遇到如下形式的约束优化问题： \(\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\)
这里存在一个著名的权衡困境：
2008-2009 年间，Goldstein 和 Osher 提出的 Split Bregman 算法改变了游戏规则。它的核心哲学是：不要试图一次性用巨大的 $\lambda$ 强行满足约束，而是用一个适中的 $\lambda$ 多次求解，每次把上一次的“违反量”加回到数据中去。
这种方法不仅速度极快（通常 10-50 次迭代即可收敛），而且由于每一步的 $\lambda$ 保持常数且适中，数值稳定性极佳。
要理解这个算法，不能只看公式，必须理解其背后的几何量——Bregman 距离。
设 $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>
| |
|----------|-------------------------> x
v u
核心解读：
考虑去噪问题 $\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$ 的递推关系），上述迭代等价于：
直觉 (Adding back the noise)： 假设 $f$ 是含噪图像。第一次迭代后，$u^1$ 变得平滑，丢失了纹理。 残差 $r = f - u^1$ 中包含了“噪声 + 丢失的纹理”。 下一步我们让算法去逼近 $f + r = f + (f-u^1)$。这就好比告诉算法：“既然你上次去噪去狠了，这次我把丢失的部分加倍补给你，请你务必保留下来。”
原始 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$ 耦合在二次项里，我们可以固定一个解另一个：
这就是著名的 Split Bregman 迭代三步曲。
让我们把上述抽象框架应用到具体的 TV 模型中。 模型：$\min_u |\nabla u|_1 + \frac{\mu}{2}|Ku - f|^2$ 其中 $K$ 是模糊核（去噪时 $K=I$）。
引入 $d_x \approx \nabla_x u, d_y \approx \nabla_y u$。 约束为：$d_x = \nabla_x u, \quad d_y = \nabla_y u$。
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，速度极快。
2. 更新 $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)
3. 更新 $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})\)
许多初学者会对 Split Bregman 和 ADMM (Alternating Direction Method of Multipliers) 感到困惑。它们看起来非常相似，但名字不同。
结论：对于凸问题，Split Bregman 算法 完全等价于 ADMM 应用在特定的增广拉格朗日形式上。
对照表 (Rosetta Stone)
| 概念 | 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)$ | 实质是梯度上升更新乘子 |
区别在哪？
分母为零 (Division by Zero)：
在各向同性 TV 的 shrinkage 公式中，出现 $\frac{s_x}{S}$。当 $S=0$（即平坦区域）时，程序会崩溃或产生 NaN。
Fix：实现时使用 max(S, 1e-10) 或者 np.where 掩码处理。
FFT 的中心化 (Shift)：
MATLAB/Python 的 fft2 默认零频在角落，而高斯模糊核或差分算子通常定义在中心。
Fix：必须使用 ifftshift 将算子调整到与 fft2 兼容的格式（零频在 (0,0)），或者显式构造 psf2otf。算子相位不对会导致重建图像发生平移或产生振铃。
$\lambda$ 选得太激进： 虽然 Split Bregman 对 $\lambda$ 鲁棒，但如果 $\lambda$ 设得极大（例如 $10^5$），子问题 $u$ 的更新几乎完全被 $\lambda \Delta$ 主导，忽略了数据项，导致收敛极其缓慢（Stiffness）。 Rule-of-thumb：对于归一化到 [0,1] 的图像，$\lambda$ 通常取在 $0.1 \mu$ 到 $10 \mu$ 之间。
各向同性 vs 各向异性混用： 如果代码里 $u$ 子问题用了拉普拉斯算子（隐含各项同性假设），而 $d$ 子问题用了简单的分开软阈值（各向异性），这会导致求解的模型方向偏差，出现细微的棋盘格伪影。
忘记更新 $b$： 这是新手最容易犯的错。如果忘了 $b$，算法就退化成了普通的交替最小化（Coordinate Descent），无法精确满足约束，导致图像看起来像被“过度磨皮”了。