第5章 鲁棒统计与非凸优化：从 M-估计到 IRLS/MM

本章摘要 在理想的教科书世界里，噪声服从高斯分布，图像先验是完美的拉普拉斯分布，优化问题总是凸的。然而，现实世界的成像充满了“意外”：传感器坏点、动态场景中的遮挡、非高斯的重尾噪声，以及 $L_1$ 正则化带来的固有偏差（Bias）。

本章我们将跨越“凸优化”的安全边界，进入非凸（Non-convex）与鲁棒（Robust）的荒原。我们将揭示如何设计“甚至能忽略 50% 错误数据”的损失函数，如何用非凸稀疏势函数恢复出比 $L_1$ 更锐利的边缘，并掌握 MM（Majorization-Minimization）、IRLS 和 HQS 这三把利剑，将棘手的非凸优化问题驯化为我们熟悉的一系列线性子问题。

5.1 为什么要鲁棒？当 $L_2$ 和 $L_1$ 失效时

5.1.1 最小二乘的“民主”与脆弱

最小二乘法（$L_2$ Loss）基于高斯噪声假设，其本质是“民主”的：它试图让所有数据点的误差平方和最小。这意味着，一个巨大的异常值（Outlier）可以绑架整个模型。

  • 场景：你正在计算一组像素的平均值作为背景。数值为 [100, 102, 98, 100]。
  • 异常：突然飞进来一只鸟，其中一个值变成了 2000。
  • $L_2$ 结果：平均值被拉高到 480。背景变成了灰色，完全错误。
  • $L_1$ 结果：中位数为 100。背景恢复正确。

看起来 $L_1$ 赢了？但在更复杂的线性回归或层析成像（Tomography）中，当异常值比例很高或具有结构性（如大面积遮挡）时，$L_1$ 的“软阈值”策略依然会受到牵连，产生收缩效应。我们需要更激进的策略——M-估计。

5.1.2 M-估计（M-estimator）与影响函数

M-估计将问题建模为 $\min_x \sum_i \rho(r_i)$，其中 $r_i$ 是残差。 衡量鲁棒性的核心指标是影响函数（Influence Function）：$\psi(r) = \rho'(r)$。它描述了单个数据点的误差如何“拉动”最终的解。

常见势函数家族谱系

我们通过 ASCII 图来直观对比三种势函数及其导数行为：

1. L2 (Quadratic)      2. L1 (Absolute)       3. Robust (e.g., Tukey/Geman)
   "民主派"              "中立派"               "精英派/排他派"

Loss function rho(r):
      |  /                  |  /                    |________
      | /                   | /                    /
      |/                    |/                    /
    --+--                  --+--                 --+-- 
   (Parabola)            (V-shape)             (Saturates at const)

Influence function psi(r) = rho'(r):
      |  /                  |__                   |  _
      | /                   |  |                  | / \
______|/______         _____|__|_____         ____|/   \________
      /                     |                     Redescending!
     /                      |__                   (Returns to 0)

   Linear Growth        Constant Pull          Zero Pull at infinity
(Outliers dominate)   (Outliers influence)    (Outliers ignored)

  1. 凸集家族 (Convex)：

    • $L_2$：$\psi(r) = r$。误差越大，拉力越大（无上限）。
    • $L_1$：$\psi(r) = \text{sign}(r)$。误差再大，拉力恒定为 1。
    • Huber：小误差是 $L_2$，大误差是 $L_1$。它是可微且凸的各种“好”性质的集合体。

  2. 红降家族 (Redescending / Non-convex)：

    • Geman-McClure：$\rho(r) = \frac{r^2}{r^2 + \sigma^2}$。
    • Tukey's Biweight：超过阈值后，$\psi(r) = 0$。
    • 特性：红降（Redescending）——当误差 $r \to \infty$ 时，$\psi(r) \to 0$。这意味着模型会完全切断与该数据点的联系，视其为无效测量。这是处理遮挡、高脉冲噪声的终极武器。

5.2 非凸稀疏：超越 $L_1$ 的偏差

在正则化项（Prior）中，非凸性同样重要。$L_1$ 范数虽然能诱导稀疏，但它是有代价的。

5.2.1 $L_1$ 的收缩偏差 (Shrinkage Bias)

$L_1$ 的近端算子是软阈值（Soft Thresholding）： $$S_\lambda(y) = \text{sign}(y) \max(|y| - \lambda, 0)$$ 观察公式可知，对于那些非常大的真实信号（$|y| \gg \lambda$），它们也被减去了 $\lambda$。

  • 后果：图像去噪后，强边缘的对比度下降；MRI 重建中，高亮区域的亮度被低估。

5.2.2 非凸稀疏势函数

为了修正偏差，我们需要一种“在这个数值很小时惩罚它，在它很大时放过它”的函数。

  1. $L_p$ 范数 ($0 < p < 1$)：

    • 形式：$|x|_p^p = \sum |x_i|^p$。
    • 几何：单位球是内凹的“星形”。
    • 效果：对小系数的抑制比 $L_1$ 更狠，对大系数的收缩比 $L_1$ 更小。

  2. SCAD 与 MCP：

    • 这两种函数是统计学界精心设计的（Folded Concave Penalties）。
    • MCP (Minimax Concave Penalty) 的导数设计为：开始是 $\lambda$，线性下降，最后变为 0。
    • 直觉：只要信号足够强，正则化项就“关掉”，实现无偏估计（Unbiasedness）。

5.3 主要化-最小化 (MM) 框架

如何优化上述非凸函数？梯度下降法容易陷入局部极小且步长难调。MM 框架提供了一种将非凸问题转化为凸问题序列的通用方法论。

5.3.1 核心哲学

不要试图一步登天直接最小化复杂的 $F(x)$。我们寻找一个简单的函数 $G(x, x_k)$（通常是二次函数），它像一个“碗”一样扣在 $F(x)$ 上方，并与其在当前点 $x_k$ 相切。

5.3.2 几何直觉 (ASCII Demo)

      G(x, x_k) [Surrogate: Quadratic "Bowl"]
          \       |       /
           \      |      /
            \     |     /
             \    |    /
              \   V   /
F(x) [Original] \_|_ /
  (Wiggly)      / | \
 ______________/  .  \_______________
              /   .   \
             x_k  .    x_{k+1}
                  .
        "Descent Guarantee"
  1. 主要化（Majorize）：$G(x, x_k) \ge F(x)$ 对所有 $x$ 成立。
  2. 相切（Tangent）：$G(x_k, x_k) = F(x_k)$。
  3. 最小化（Minimize）：算出 $G$ 的谷底 $x_{k+1}$。
  4. 结论：$F(x_{k+1}) \le G(x_{k+1}, x_k) \le G(x_k, x_k) = F(x_k)$。函数值单调下降。

5.4 迭代重加权最小二乘 (IRLS)

IRLS 是 MM 算法在图像处理中最广泛的应用形式，专门用于解决 $L_p$ 范数 ($0 < p \le 1$) 和 M-估计问题。

5.4.1 从梯度到权重

考虑问题 $\min_x \sum_i \rho((Ax-y)_i)$。 令其梯度为 0：$A^T \rho'(Ax-y) = 0$。 利用恒等式 $\rho'(r) = \frac{\rho'(r)}{r} \cdot r$，我们定义权重 $w(r) = \frac{\rho'(r)}{r}$。 方程变为加权正规方程（Weighted Normal Equation）： $$ A^T W(x) (Ax - y) = 0 $$ 其中 $W(x)$ 是对角矩阵，$\text{diag}(W)_i = w((Ax-y)_i)$。

5.4.2 常见权重速查

| 模型 | 势函数 $\rho(r)$ | 权重 $w(r)$ | 物理含义 |

模型 势函数 $\rho(r)$ 权重 $w(r)$ 物理含义
Gaussian ($L_2$) $r^2$ $1$ 所有数据平等
Laplace ($L_1$) $\lvert r\rvert$ $1/\lvert r\rvert$ 误差越大，权重越小
$L_p$ ($0<p<1$) $\lvert r\rvert^p$ $1/\lvert r\rvert^{2-p}$ 极度抑制大误差/稀疏化
Geman-McClure $\frac{r^2}{r^2+\sigma^2}$ $\frac{2\sigma^2}{(r^2+\sigma^2)^2}$ 红降权重（快速衰减至0）

5.4.3 算法实现细节

IRLS 算法流程：

  1. 初始化 $x^{(0)}$，设定正则化参数 $\epsilon$（极小正数）。
  2. For $k = 0, 1, \dots$:
    • 计算残差向量：$r = Ax^{(k)} - y$。
    • 更新权重向量：$w_i = ((r_i)^2 + \epsilon^2)^{\frac{p-2}{2}}$ （以 $L_p$ 为例）。
    • 求解线性系统：$x^{(k+1)} = (A^T \text{diag}(w) A)^{-1} A^T \text{diag}(w) y$。
      • 注：如果包含正则项，公式相应变为 $(A^T W_{data} A + \lambda W_{reg})^{-1} \dots$

Rule of Thumb: IRLS 的收敛速度通常是超线性的（比一阶方法快），但每次迭代需要解线性方程组。如果 $A$ 是卷积，可用 FFT 快速求解；如果是稀疏矩阵，用预条件共轭梯度（PCG）。

5.5 半二次分裂 (Half-Quadratic Splitting, HQS)

当目标函数包含难以处理的非凸正则项 $\rho(x)$ 时，HQS 提供了一种将其与数据项解耦的策略。它也是现代 Plug-and-Play (PnP) 方法的前身。

5.5.1 原理：变分形式

许多势函数 $\rho(x)$ 可以写成下确界形式（基于凸共轭理论的推广）： $$ \rho(x) = \inf_{z} \left( \frac{\mu}{2} (x - z)^2 + \psi(z) \right) $$ 这里引入了辅助变量 $z$ 和惩罚参数 $\mu$。 原始问题 $\min_x |Ax-y|^2 + \lambda \rho(x)$ 被转化为： $$ \min_{x, z} |Ax-y|^2 + \frac{\lambda \mu}{2} |x - z|^2 + \lambda \psi(z) $$

5.5.2 求解：交替最小化

我们将原来的难问题拆解为两个子问题：

  1. $x$-子问题（去卷积/逆问题）： $$ x_{k+1} = \arg\min_x |Ax-y|^2 + \frac{\lambda \mu}{2} |x - z_k|^2 $$
  • 这是一个标准的 $L_2$ 正则化最小二乘问题。
    • 解法：若 $A$ 为卷积，FFT 闭式解；否则用 CG 求解。

  1. $z$-子问题（去噪/近端映射）： $$ z_{k+1} = \arg\min_z \frac{\mu}{2} |x_{k+1} - z|^2 + \psi(z) $$

    • 这是一个逐点（Point-wise）优化问题。
    • 等价于求解 $Prox_{\psi/\mu}(x_{k+1})$。
    • 关键点：即使 $\rho$ 是非凸的，只要它的近端算子容易算（如硬阈值 Hard Thresholding），这一步就很简单。

5.5.3 Continuation 策略

在 HQS 中，$\mu$ 通常不是固定的。

  • 开始时，设 $\mu$ 较小，允许 $x$ 和 $z$ 差异较大（探索全局）。
  • 随着迭代，逐渐增大 $\mu$ ($\mu \to \infty$)，强制 $x$ 逼近 $z$。 这种策略不仅加速收敛，还能帮助算法跳出浅层的局部极小值。

5.6 非凸优化的陷阱与生存指南

5.6.1 局部最优（Local Minima）

非凸函数像一个充满坑洼的山谷，梯度下降很容易卡在半山腰的坑里。 生存指南：渐进非凸性 (Graduated Non-Convexity, GNC)

  • 不要一开始就优化非凸目标。
  • Step 1: 构造一个凸的近似函数（如用 $L_2$ 或 $L_1$），求得初始解 $x_{init}$。
  • Step 2: 稍微增加非凸性（如 Huber 的阈值变小，或 $L_p$ 的 $p$ 从 1 降到 0.9），以 $x_{init}$ 为起点优化。
  • Step 3: 重复此过程，直到目标函数变成真正的非凸形式。

5.6.2 参数敏感性

  • $\epsilon$ (IRLS): 防止除零的平滑参数。太大 $\to$ 变成 $L_2$；太小 $\to$ 数值溢出。
  • $\mu$ (HQS): 惩罚参数。增长太快 $\to$ 锁定在伪解；增长太慢 $\to$ 不收敛。

5.7 本章小结

| 概念 | 核心思想 | 适用场景 | 关键公式/算子 |

概念 核心思想 适用场景 关键公式/算子
Robust Statistics 抑制离群值，切断大误差的影响 椒盐噪声、遮挡去除、背景建模 Huber, Tukey, Welsch
$L_p$ ($p<1$) 减少 $L_1$ 的偏差，获得更稀疏的解 压缩感知、超分、去模糊 $\lvert x\rvert^p$, SCAD
IRLS 将非凸/非线性转化为加权最小二乘序列 $L_p$ 优化, M-估计 $W = \text{diag}(1/\lvert r\rvert)$
HQS 变量分裂，分离数据项与先验项 图像复原, PnP 框架 $x$-step (Inversion), $z$-step (Prox)
MM 构造凸上界函数逼近非凸函数 通用非凸优化框架 Surrogate Function

5.8 练习题

基础题

  1. 权重计算： 给定残差向量 $r = [0.1, 0.001, 5.0]$。 (a) 计算 $L_2$ 权重。 (b) 计算 $L_1$ 权重（设 $\epsilon=0$）。 (c) 计算 Tukey Bisquare 权重（设阈值 $c=4.0$）。哪个数据点被“丢弃”了？

    Hint (a) 全为 1。(b) $10, 1000, 0.2$。(c) 5.0 的残差超过 4.0，权重为 0。

  2. HQS 解析解： 在 HQS 的 $x$-更新步中，如果不使用 FFT，而是针对逐像素问题（如仅有数据项和恒等变换），目标函数为 $\min_x (x-y)^2 + \mu (x-z)^2$。求 $x$ 的解析解。

    Hint 这是加权平均：$x = \frac{y + \mu z}{1 + \mu}$。这展示了 $x$ 是观测值 $y$ 和辅助变量 $z$ 的折中。

  3. 非凸几何： 画出 $f(x) = |x|^{0.5}$ 的图像。说明为什么在该函数上做梯度下降，当 $x$ 接近 0 时，梯度会变得无穷大？这不仅是数值问题，也是其强稀疏诱导能力的来源。

    Hint $x^{0.5}$ 在原点处是垂直切线。巨大的梯度会将微小的非零值迅速推向 0，或者反弹（需要步长控制）。

挑战题

  1. MM 构造： 对于非凸函数 $f(x) = \log(1 + x^2)$（Cauchy/Lorentzian loss），请利用不等式 $\log(1+x^2) \le \log(1+x_k^2) + \frac{x^2 - x_k^2}{1+x_k^2}$ 构造一个二次替代函数 $Q(x, x_k)$。证明这个替代函数满足 MM 的条件。

    Hint 基于函数的凹凸性或泰勒展开的上界。你需要证明该抛物线在 $x_k$ 处与原函数相切，且开口曲率足以包络原函数。

  2. IRLS 的死锁： 如果使用 $L_1$ 范数的 IRLS 求解 $\min |x|$，初始值 $x_0 = 0$。计算第一次迭代的权重 $w$。会发生什么？这说明初始化需要注意什么？

    Hint $w \approx 1/\epsilon$。如果 $\epsilon$ 极小，权重极大，数值不稳定。更重要的是，如果初值恰好在“坑”里，算法可能无法跳出。通常不建议从 0 初始化 IRLS。

5.9 常见陷阱与错误 (Gotchas)

  1. 混淆鲁棒数据项与稀疏正则项：

    • 数据项 (Data Fidelity) 用 M-estimator（如 Huber）：是为了抗拒错误的观测 $y$。
    • 正则项 (Regularization) 用非凸势（如 $L_p$）：是为了诱导解 $x$ 的某些性质（如锐利边缘）。
    • 不要搞反了。在正则项上用 Geman-McClure 可能会导致图像出现奇怪的阶梯或断裂，除非你确实想要分段常数解。

  2. IRLS 的 $\epsilon$ 依赖症： 很多初学者写 IRLS 代码时，为了防止除以零，将 $\epsilon$ 设为固定值（如 $1e-3$）。这实际上改变了目标函数，把 $L_1$ 变成了 Huber。如果你发现结果不够稀疏，尝试使用动态 $\epsilon$ 策略：每次迭代令 $\epsilon_{k+1} = \max(\epsilon_k / \eta, 10^{-8})$。

  3. 非凸优化的“随机性”： 如果在非凸优化中发现结果每次跑都不一样，或者有很多噪点，通常是因为陷入了局部极小。必须使用 Warm Start（热启动）：先跑 50 次 $L_2$ 或 $L_1$ 迭代，再切换到非凸 Loss。

  4. HQS 中的 FFT 陷阱： 在 HQS 的 $x$-update 步骤中，通常形式是 $(A^T A + \mu I)^{-1}$。

    • 只有当 $A$ 是循环卷积（对应周期边界条件）时，才能用 FFT 对角化。
    • 如果 $A$ 包含掩膜（Masking/Inpainting）或非周期边界，直接用 FFT 是错误的！此时需改用 PCG 或 ADMM。