附录 A：近端算子与投影算子完全手册 (The Proximal Operator Handbook)

摘要：本附录是变分图像处理算法实现的“字典”。在 ADMM、原始-对偶 (Primal-Dual)、PnP (Plug-and-Play) 等算法中，核心步骤往往归结为求解某个特定函数的近端算子。

使用指南：

• 遇到 $\min \frac{1}{2}|x-v|^2 + \lambda f(x)$ 子问题时，查阅此表。
• 包含 Python/NumPy 参考实现，可直接用于 def 函数编写。

1. 基础理论与几何直觉

1.1 定义

对于闭凸函数 $f: \mathbb{R}^n \rightarrow \mathbb{R} \cup \{+\infty\}$ 和参数 $\gamma > 0$（通常对应步长或惩罚参数的倒数），近端算子定义为：

$$ \text{prox}_{\gamma f}(v) = \mathop{\arg\min}_{x} \left( \underbrace{f(x)}_{\text{正则项}} + \underbrace{\frac{1}{2\gamma}|x - v|_2^2}_{\text{惯性项/数据项}} \right) $$

1.2 几何直觉：后向欧拉步骤

梯度下降 $x_{k+1} = x_k - \gamma \nabla f(x_k)$ 可以看作是前向欧拉离散化。而近端算子实际上是隐式梯度下降（后向欧拉）的一步： $$ x = v - \gamma \nabla f(x) \iff x = (I + \gamma \nabla f)^{-1}(v) = \text{prox}_{\gamma f}(v) $$ 即使 $f$ 不可微（没有梯度），近端算子依然定义良好，这正是它在图像处理中处理稀疏（$\ell_1$）和边缘（TV）的关键原因。

1.3 莫罗分解 (Moreau Decomposition)

这是一个强大的工具，通过函数的共轭 $f^*$ 来计算 $f$ 的 Prox： $$ v = \text{prox}_{\gamma f}(v) + \gamma \text{prox}_{\frac{1}{\gamma} f^*}\left(\frac{v}{\gamma}\right) $$ 应用场景：当 $f$ 本身很难求 Prox，但其对偶范数的投影球很简单时（例如 TV 正则化求解中，我们将问题转化为向对偶球的投影）。

2. 稀疏与元素级算子 (Element-wise Operators)

这类算子处理图像的每个像素或变换域系数，彼此独立，极易并行化。

2.1 $\ell_1$ 范数：软阈值 (Soft Thresholding)

模型：稀疏编码、小波去噪、压缩感知 $$ f(x) = |x|_1 = \sum_i |x_i| $$ 算子公式： $$ [\text{prox}_{\gamma |\cdot|_1}(v)]_i = \mathcal{S}_\gamma(v_i) = \text{sign}(v_i) \cdot \max(|v_i| - \gamma, 0) $$ 物理含义： 将数值向零收缩。如果信号幅度小于噪声水平 $\gamma$，则认为是纯噪声并置零；否则保留信号但减去 $\gamma$ 偏差。

Python 实现：

def soft_threshold(v, gamma):
    # 利用 numpy 的广播机制，处理复数需小心
    return np.sign(v) * np.maximum(np.abs(v) - gamma, 0)

2.2 $\ell_0$ 伪范数：硬阈值 (Hard Thresholding)

模型：非凸稀疏恢复、背景扣除 $$ f(x) = |x|_0 $$ 算子公式： $$ [\text{prox}_{\gamma |\cdot|_0}(v)]_i = \begin{cases} v_i & |v_i| > \sqrt{2\gamma} \\ 0 & |v_i| \le \sqrt{2\gamma} \end{cases} $$ 注意：这里的阈值门限是 $\sqrt{2\gamma}$。如果目标函数写作 $\frac{1}{2}|x-v|^2 + \lambda |x|_0$，则门限为 $\sqrt{2\lambda}$。

Python 实现：

def hard_threshold(v, gamma):
    threshold = np.sqrt(2 * gamma)
    return v * (np.abs(v) > threshold)

2.3 弹性网 (Elastic Net)

模型：相关特征选择、稳定稀疏恢复 $$ f(x) = \alpha |x|_1 + \beta |x|_2^2 $$ 算子公式： 先做软阈值，再做缩放。 $$ \text{prox}_{\gamma f}(v) = \frac{1}{1 + 2\gamma\beta} \mathcal{S}_{\gamma\alpha}(v) $$

3. 结构化稀疏与向量算子 (Structured Sparsity)

处理像素组、彩色通道或梯度向量。

3.1 $\ell_{2,1}$ 范数 (Group Lasso)

模型：彩色图像去噪（RGB 联合稀疏）、各向同性 TV $$ f(x) = \sum_{g \in \mathcal{G}} |x_g|_2 $$ 其中 $x_g$ 是分组向量（例如一个像素的 RGB 三个分量，或梯度的 $(D_x, D_y)$ 分量）。

算子公式（组软阈值）： $$ \text{prox}_{\gamma |\cdot|_{2,1}}(v)_g = \max\left( 0, 1 - \frac{\gamma}{|v_g|_2} \right) v_g $$ 物理含义： 不是单独处理每个分量，而是看整个“组”的能量。如果组能量太小，整个组被置零；否则，组内向量方向不变，模长收缩。这保证了彩色去噪时的色相（Hue）不变性。

Python 实现：

def group_soft_threshold(v, gamma):
    # 假设 v 的形状为 (Height, Width, Channels)，按 Channel 分组
    # norm 的形状为 (Height, Width, 1)
    norm = np.linalg.norm(v, axis=-1, keepdims=True)
    # 避免除以 0：如果 norm < epsilon，系数设为 0 即可
    scale = np.maximum(0, 1 - gamma / (norm + 1e-10))
    return v * scale

3.2 $\ell_2$ 球约束投影

模型：去噪中的噪声水平约束 $|Ax - y|_2 \le \varepsilon$ $$ C = \{x \mid |x|_2 \le \varepsilon\} $$ 算子公式： $$ \text{proj}_C(v) = \min\left( 1, \frac{\varepsilon}{|v|_2} \right) v $$

4. 矩阵与谱算子 (Matrix & Spectral Operators)

图像被视为矩阵 $X \in \mathbb{R}^{m \times n}$，常用于低秩建模。

4.1 核范数 (Nuclear Norm / Trace Norm)

模型：Robust PCA、矩阵补全、低秩去噪 $$ f(X) = |X|_* = \sum_{i} \sigma_i(X) $$ 算子公式（奇异值阈值 SVT）： 设 $v$ 的奇异值分解为 $U \Sigma V^T$。 $$ \text{prox}_{\gamma |\cdot|_*}(v) = U \cdot \mathcal{S}_\gamma(\Sigma) \cdot V^T $$ 即：对奇异值做软阈值处理，保持奇异向量不变。

工程陷阱： SVD 计算极其昂贵 ($O(\min(m,n)^3)$)。在实际应用（如 RPCA）中，通常只计算部分奇异值（使用 PROPACK 或 Lanczos 方法）以加速。

Python 实现：

def svt(X, gamma):
    U, S, Vt = np.linalg.svd(X, full_matrices=False)
    S_new = np.maximum(S - gamma, 0)
    # 重建矩阵，注意利用广播
    return (U * S_new) @ Vt

4.2 秩约束 (Rank Constraint)

模型：固定秩近似 $$ C = \{X \mid \text{rank}(X) \le r\} $$ 算子公式： 保留前 $r$ 个最大的奇异值，其余置零。 $$ \text{proj}_C(X) = \sum_{i=1}^r \sigma_i u_i v_i^T $$

5. 常见数据项与非标准正则项

5.1 KL 散度 (Poisson Noise)

模型：低光照图像恢复、医学成像 (PET/CT) $$ f(x) = \sum (x_i - y_i \log x_i + \text{const}) $$ 此处 $y$ 是观测到的光子计数。

算子公式： 求解方程 $x - y + \frac{1}{\gamma}(x - v) = 0$ 的正根。 $$ \text{prox}_{\gamma f}(v) = \frac{(v-\gamma) + \sqrt{(v-\gamma)^2 + 4\gamma y}}{2} $$ 注意：通常 $x$ 需限制为非负。

5.2 Huber 损失

模型：鲁棒回归、梯度平滑（减少阶梯效应） $$ L_\delta(x) = \begin{cases} \frac{1}{2}x^2 & |x| \le \delta \\ \delta(|x| - \frac{\delta}{2}) & |x| > \delta \end{cases} $$ 算子公式： 这是一个“平滑收缩”算子，比软阈值更平滑。 $$ \text{prox}_{\gamma L_\delta}(v) = \frac{\delta}{\delta + \gamma} v + \frac{\gamma}{\delta + \gamma} \mathcal{S}_{\delta}(v) $$ 或者更简单的形式：基于 $v$ 的大小分段缩放。

6. 约束集的投影 (Indicator Functions)

形式为 $\min \iota_C(x)$，对应 Prox 为 $P_C(x)$。

| 约束类型 | 集合 $C$ | 投影算子 $\text{proj}_C(v)$ | 应用场景 |

7. 高级技巧：Proximal Calculus

当你遇到复杂的函数组合时，不需要每次都用定义推导。

7.1 缩放与平移

若 $g(x) = f(ax + b)$，其 Prox 通常很难直接写出，除非 $a$ 是正交矩阵。 但若 $g(x) = \lambda f(x)$，则 $\text{prox}_{\gamma g} = \text{prox}_{(\gamma\lambda) f}$。 一定要分清算法步长 $\rho$ 和正则参数 $\lambda$。在 ADMM 代码中，通常调用 prox(v, lambda/rho)。

7.2 求解 Total Variation (TV) 的 Prox

$$ f(x) = |\nabla x|_1 $$ 陷阱：二维 TV 的 Prox 没有简单的闭式解！ 解决方案：

1. 使用 ADMM 分裂：引入 $z = \nabla x$，将问题变为求解 $z$ 的 $\ell_1$ Prox（软阈值）和 $x$ 的泊松方程（FFT 求解）。
2. 对偶方法：使用 Chambolle's Algorithm（基于对偶球投影的定点迭代），通常迭代 5-10 次即可获得高精度近似。

8. 代码调试自检清单 (Debug Checklist)

在实现上述算子时，若算法不收敛，请检查：

1. 参数缩放：确认传入 soft_threshold 的阈值是 $\lambda$ 还是 $\lambda/\rho$？对于 ADMM，这通常是错误的根源。
2. 复数处理：处理 FFT 后的数据或 MRI 数据时，abs(v) 返回模长，sign(v) 返回相位向量 $e^{i\phi}$。确保逻辑正确保留了相位。
3. 除零保护：在 Group Lasso 或 TV 投影中计算 $v / |v|_2$ 时，务必加上 eps ($10^{-10}$)。
4. 维度检查：在做 SVD 或 Group 操作时，明确是对哪个轴（axis）操作。对图像矩阵做 norm 可能会意外变成标量而不是向量。
5. SVD 精度：np.linalg.svd 默认全分解。对于大矩阵，应使用 scipy.sparse.linalg.svds (ARPACK) 仅计算前 $k$ 个奇异值。