opt_vision_tutorial

第7章算子分裂 I：ADMM、Douglas–Rachford 与变量分裂

学习目标：

核心逻辑：理解如何通过“引入变量”将一个不可微、耦合严重的优化难题，拆解为两个可微或有闭式解的简单子问题。

数学推导：掌握从增广拉格朗日（ALM）到 ADMM 缩放形式（Scaled Form）的完整推导过程。

求解技术：熟练运用 FFT 和矩阵求逆引理（Sherman-Morrison）快速求解 ADMM 中的线性子问题。

理论联系：理解 ADMM 与 Douglas-Rachford 分裂法的几何联系及其收敛性直觉。

7.1 引言：不可微与线性算子的“那扇窄门”

在图像复原与重建中，我们经常面临如下结构的复合优化问题： $\min_x \underbrace{f(x)}_{\text{数据保真}} + \underbrace{g(Ax)}_{\text{正则化}}$

$f(x)$ 通常是光滑的凸函数（如 $\frac{1}{2}|x-y|^2$）。
$g(\cdot)$ 通常是非光滑的凸函数（如 $\ell_1$ 范数 $|\cdot|_1$）。
$A$ 是一个线性算子（如梯度算子 $\nabla$、小波变换 $W$ 或卷积 $K$）。

困境：如果 $A=I$（单位阵），问题就是近端梯度法（Proximal Gradient）的主场。但当 $A$ 存在且不可逆时（例如全变分 TV 中的梯度算子），$g(Ax)$ 的近端算子 $\text{prox}_{g(A\cdot)}$ 通常没有闭式解。$A$ 就像扭曲了空间的透镜，让简单的 $\ell_1$ 球变得极其复杂。

策略：分而治之（Divide and Conquer） 算子分裂的核心思想极其朴素：既然 $g$ 和 $A$ 纠缠在一起很难解，那就引入一个新变量 $z$ 把它们分开，然后强迫新旧变量相等。

这就像两个人要同时侧身穿过一扇窄门会卡住（难解），我们引入一个“分身”，让真身处理 $f(x)$ 和线性变换，分身处理复杂的非光滑项 $g(z)$，最后通过“弹簧”（二次惩罚）和“胶水”（拉格朗日乘子）让真身和分身不断同步。

7.2 变量分裂与增广拉格朗日 (ALM)

7.2.1 变量分裂的标准形式

我们将原问题重写为约束优化问题： $\min_{x,z} f(x) + g(z) \quad \text{s.t.} \quad Ax - z = 0$ 此时，目标函数完全解耦（separable），但约束条件将它们重新绑在了一起。

7.2.2 为什么要用“增广”拉格朗日？

传统的拉格朗日乘子法 $\mathcal{L}(x,z,y) = f(x) + g(z) + y^T(Ax-z)$ 虽然理论上可行，但在数值上极其不稳定，且对对偶变量 $y$ 的更新步长非常敏感。

增广拉格朗日（Augmented Lagrangian） 引入了一个二次惩罚项 $(\rho/2)|Ax-z|^2$： $\mathcal{L}_{\rho}(x, z, y) = f(x) + g(z) + y^T(Ax - z) + \frac{\rho}{2}\|Ax - z\|_2^2$ 直觉：

拉格朗日项 ($y^T\dots$)：负责将解推向可行域（满足约束），保证收敛到精确解。
二次惩罚项 ($\rho \dots$)：像一根强力的弹簧，将 $Ax$ 和 $z$ 拉近。它为目标函数增加了强凸性，极大地提升了数值稳定性和收敛速度。

7.2.3 缩放形式（Scaled Form）：实现的基石

在实际工程代码中，带着 $y^T(Ax-z)$ 这一项进行配方非常繁琐。我们通常使用 缩放对偶变量 $u$。

推导：观察 $\mathcal{L}_{\rho}$ 的后两项： $\begin{aligned} y^T(Ax - z) + \frac{\rho}{2}\|Ax - z\|^2 &= \frac{\rho}{2} \left( \frac{2}{\rho}y^T(Ax - z) + \|Ax - z\|^2 \right) \\ &= \frac{\rho}{2} \left( \|Ax - z + \frac{1}{\rho}y\|^2 - \|\frac{1}{\rho}y\|^2 \right) \end{aligned}$ 令 $u = \frac{1}{\rho}y$，这就是缩放对偶变量。忽略常数项 $-\frac{\rho}{2}|u|^2$，我们的优化目标变为更简洁的形式： $\mathcal{L}_{\rho}(x, z, u) = f(x) + g(z) + \frac{\rho}{2} \|Ax - z + u\|_2^2$

Rule of Thumb: 永远在代码中使用缩放形式。它将所有约束相关项折叠进了一个简单的 $\ell_2$ 范数中，使得子问题看起来就像是标准的去噪问题（Proximal problem）。

7.3 ADMM：交替方向乘子法

ADMM (Alternating Direction Method of Multipliers) 的精髓在于“交替”。如果同时最小化 $x$ 和 $z$，问题依然很难。但如果我们固定一个，更新另一个，事情就变得简单了。

7.3.1 算法流程

给定 $k$ 时刻的 $z^k, u^k$，迭代如下：

x-最小化步骤 (x-update)： $x^{k+1} = \arg\min_x \left( f(x) + \frac{\rho}{2} \|Ax - z^k + u^k\|_2^2 \right)$
- 这是一个关于 $x$ 的平滑凸优化问题（如果 $f$ 是平滑的），通常对应求解线性方程组。
z-最小化步骤 (z-update)： $z^{k+1} = \arg\min_z \left( g(z) + \frac{\rho}{2} \|Ax^{k+1} - z + u^k\|_2^2 \right)$
- 注意这里使用了最新的 $x^{k+1}$。这正是 Proximal Operator 的定义形式： $z^{k+1} = \text{prox}_{g, 1/\rho} (Ax^{k+1} + u^k)$
对偶变量更新 (u-update)： $u^{k+1} = u^k + (Ax^{k+1} - z^{k+1})$
- 这就是梯度上升法更新对偶变量（步长为 $\rho$，但在缩放形式中被隐藏了）。本质是累积残差。

7.3.2 流程与数据流图 (ASCII Flow)

      [Start: x0, z0, u0, rho]
                 |
      +---> [ x-update ] <--------------------------------+
      |     Solve: min f(x) + (rho/2)||Ax - (z-u)||^2     |
      |     (Often a Linear Solver / FFT)                 |
      |          |                                        |
      |     (Pass x_new)                                  |
      |          v                                        |
      |     [ z-update ]                                  |
      |     Solve: min g(z) + (rho/2)||z - (Ax+u)||^2     |
      |     (Often a Prox Operator / Thresholding)        |
      |          |                                        |
      |     (Pass z_new)                                  |
      |          v                                        |
      |     [ u-update ]                                  |
      |     u_new = u_old + (Ax_new - z_new)              |
      |     (Accumulate Residual)                         |
      |          |                                        |
      +----------+                                        |
                 |                                        |
         [Check Convergence] -----------------------------+
         If primal_res < eps AND dual_res < eps
                 |
               [End]

7.4 典型成像应用与子问题求解技巧

ADMM 的威力取决于能否高效求解 $x$ 和 $z$ 的子问题。在图像处理中，我们有两大法宝：FFT 和 解析解。

7.4.1 案例一：TV 去噪 (FFT 加速)

模型：$\min_x \frac{1}{2}|x-y|_2^2 + \lambda |\nabla x|_1$ 令 $A=\nabla$（梯度算子），分裂变量 $z=\nabla x$。

x-update: $\min_x \frac{1}{2}\|x-y\|_2^2 + \frac{\rho}{2}\|\nabla x - z^k + u^k\|_2^2$ 求导令其为 0，得到线性方程组： $(I + \rho \nabla^T \nabla) x = y + \rho \nabla^T(z^k - u^k)$ 这看起来涉及大矩阵求逆。但在周期边界条件假设下，$\nabla$ 是循环矩阵卷积，$\nabla^T \nabla$ 也是卷积。根据卷积定理，频域中的卷积等于点乘。

FFT 解法： $x^{k+1} = \mathcal{F}^{-1} \left( \frac{\mathcal{F}(y) + \rho \mathcal{F}(\nabla^T)\odot \mathcal{F}(z^k - u^k)}{ \mathbf{1} + \rho (|\mathcal{F}(\nabla_h)|^2 + |\mathcal{F}(\nabla_v)|^2) } \right)$
- 分母可以在循环外预计算（Precomputation）。
- 每次迭代只需做正反 FFT，复杂度 $O(N \log N)$。
z-update: $\min_z \lambda \|z\|_1 + \frac{\rho}{2} \|\nabla x^{k+1} - z + u^k\|_2^2$ 这是一个标准的 $\ell_1$ 极小化，解为 软阈值（Soft Thresholding）： $z^{k+1} = \text{soft}(\nabla x^{k+1} + u^k, \frac{\lambda}{\rho})$ 其中 $\text{soft}(v, \tau) = \text{sign}(v) \max(|v| - \tau, 0)$，逐元素操作。

7.4.2 案例二：去模糊与压缩感知 (矩阵求逆引理)

模型：$\min_x \frac{1}{2}|Kx-y|_2^2 + \lambda |x|_1$ 这里 $K$ 是模糊核或降采样矩阵。这次我们令 $z=x$（将 $x$ 从 $\ell_1$ 中解救出来）。

x-update: $\min_x \frac{1}{2}\|Kx-y\|_2^2 + \frac{\rho}{2}\|x - z^k + u^k\|_2^2$ 正规方程： $(K^T K + \rho I) x = K^T y + \rho(z^k - u^k)$
- 场景 A (去模糊)：$K$ 是卷积。同样使用 FFT，分母为 $ \mathcal{F}(K) ^2 + \rho$。
- 场景 B (压缩感知/修复)：$K$ 是随机矩阵或掩膜，$K$ 只有 $M$ 行 ($M \ll N$)。直接求逆 $N \times N$ 的 $(K^T K + \rho I)$ 不可能。
Sherman-Morrison-Woodbury (SMW) 技巧： $(K^T K + \rho I)^{-1} = \frac{1}{\rho}I - \frac{1}{\rho^2} K^T (I + \frac{1}{\rho} K K^T)^{-1} K$ 现在我们只需要对 $M \times M$ 的矩阵 $(I + \frac{1}{\rho} K K^T)$ 求逆。如果是掩膜（Masking），$KK^T$ 甚至是把对角阵，求逆极其简单。

7.5 Douglas-Rachford (DR) 分裂：理论渊源

Douglas-Rachford 分裂法历史更悠久，它主要用于寻找两个最大单调算子之和的零点。 $0 \in (\partial f + \partial g)(x)$

ADMM 实际上是 DR 算法在对偶问题上的等价形式。了解 DR 有助于理解几何直觉。

7.5.1 反射算子 (Reflection Operator)

定义 Prox 算子为 $P_{\gamma f} = (I + \gamma \partial f)^{-1}$。定义 反射算子 为： $R_{\gamma f} = 2P_{\gamma f} - I$ 几何上，如果不考虑非线性，这就像是关于某个平面的镜像反射。

7.5.2 DR 迭代

DR 算法迭代变量 $t$（它不是原问题的解，而是辅助变量）： $t^{k+1} = \frac{1}{2} t^k + \frac{1}{2} R_{\gamma g} ( R_{\gamma f} (t^k) )$ 这个公式极具美感：先做 $f$ 的反射，再做 $g$ 的反射，最后取平均。 如果 $f$ 和 $g$ 对应的集合有交集，这种“反射平均”最终会稳定在交集附近。

虽然我们在工程中更多直接写 ADMM 代码，但在分析收敛性证明（见附录B）时，DR 的框架提供了强大的非扩张算子（Non-expansive Operator）不动点理论支持。

7.6 关键参数 $\rho$ 的选择与调试

参数 $\rho$ 是 ADMM 的“步长”，它不改变解的位置，但决定了收敛速度。

7.6.1 现象分析

$\rho$ 太小：惩罚项 $\frac{\rho}{2}|Ax-z|^2$ 弱。$x$ 和 $z$ 各自最小化得很好，但 $Ax$ 和 $z$ 差距很大（原始残差大）。
$\rho$ 太大：惩罚项太强。$Ax$ 被死死锁在 $z$ 上，但目标函数 $f(x)+g(z)$ 难以下降（对偶残差大）。

7.6.2 残差平衡策略 (Residual Balancing)

一种非常实用的自适应策略是让原始残差 $r$ 和对偶残差 $s$ 保持在同一数量级。

计算残差：

原始残差：$r^{k+1} = Ax^{k+1} - z^{k+1}$
对偶残差：$s^{k+1} = \rho A^T (z^{k+1} - z^k)$

更新 $\rho$（每隔 10 或 20 次迭代检测一次）：

如果 $|r|_2 > 10 |s|_2$，说明约束太松：$\rho \leftarrow 2\rho$。
如果 $|s|_2 > 10 |r|_2$，说明约束太紧：$\rho \leftarrow \rho / 2$。
注意：一旦改变 $\rho$，必须同时更新 $u \leftarrow u / 2$ 或 $u \leftarrow 2u$，以保持拉格朗日乘子 $y=\rho u$ 的连续性。

7.7 本章小结

解耦哲学：ADMM 不是一个新的优化算法，而是一个框架。它把 $f(x)+g(Ax)$ 拆解为三个步骤：求解 $x$（通常是线性反演）、求解 $z$（通常是 Prox）、同步 $u$。
子问题易解性：ADMM 的效率完全取决于子问题是否易解。如果是 FFT 可解的卷积，或是有解析解的软阈值，ADMM 就是神器。如果是复杂的非线性算子，ADMM 可能会很慢。
缩放形式：工程实现中请务必使用 Scaled Form ($u$-variable)，它简化了公式，减少了出错概率。
收敛特征：ADMM 通常能很快收敛到中等精度的解（视觉上看不出区别），但收敛到高精度（机器精度）非常慢。对于图像处理，这通常足够了。

7.8 练习题

基础题（熟悉流程）

Q1 (ADMM 迭代推导) 考虑问题 $\min_x \frac{1}{2}|x-b|_2^2 + \mathbb{I}_C(x)$，其中 $C$ 是凸集 ${x \mid Lx \le c}$。请引入变量 $z$，写出 ADMM 的标准迭代步骤。

Hint: 令 $z = Lx$，约束变为 $z \le c$。$g(z)$ 是什么？

点击查看答案思路

* 模型：$\min \frac{1}{2}\|x-b\|^2 + g(z)$ s.t. $Lx-z=0$。 * $g(z)$ 是指示函数：如果 $z \le c$ 值为 0，否则为 $+\infty$。 * $x$-update: $\min \frac{1}{2}\|x-b\|^2 + \frac{\rho}{2}\|Lx - z + u\|^2$。线性方程：$(I + \rho L^T L)x = b + \rho L^T(z-u)$。 * $z$-update: $\min \mathbb{I}_{(-\infty, c]}(z) + \frac{\rho}{2}\|Lx - z + u\|^2$。这是向区域 $(-\infty, c]$ 的投影。即 $z_i = \min((Lx+u)_i, c_i)$。

Q2 (软阈值的性质) 证明软阈值算子 $\text{soft}(v, \tau)$ 对于 $\tau$ 是连续的，但不可微。这对 z-update 的数值稳定性有什么影响？

Hint: 关注 $v = \tau$ 和 $v = -\tau$ 处。

挑战题（深入思考）

Q3 (Color Image Denoising) 假设处理彩色图像 $x \in \mathbb{R}^{H \times W \times 3}$。我们希望正则化项鼓励“颜色梯度的对齐”（即边缘处 RGB 同时变化），使用 Group Lasso 形式的 TV： $\|\nabla x\|_{2,1} = \sum_{i,j} \sqrt{(\nabla x_{R})_{ij}^2 + (\nabla x_{G})_{ij}^2 + (\nabla x_{B})_{ij}^2}$ 请推导该模型下 z-update 的闭式解。这被称为 Vector Soft Thresholding。

Hint: $z$ 现在是一个向量场。Prox 问题变为 $\min_z \lambda |z|_2 + \frac{\rho}{2}|v - z|_2^2$（这里的范数是对每个像素的 vector 做）。