第17章 工程实现与复现实验：从“公式”到“可跑代码”

本章摘要：在论文中，作者常说“我们使用 ADMM 求解公式 (5)”。然而，从公式 (5) 到一套稳健、快速的代码，中间隔着无数的工程细节。本章不教你写 Hello World，而是教你如何构建隐式算子库、如何用幂法自动确定步长、如何利用矩阵求逆引理加速子问题求解，以及当结果全黑或发散时，如何像外科医生一样进行调试诊断。

17.1 线性算子编程范式：从矩阵到“动作”

在图像复原与重建中，观测模型 $y = Ax + n$ 中的 $A$ 往往过于巨大，无法显式存储。我们需要将其抽象为线性算子（Linear Operator）。

17.1.1 显式矩阵 vs. 隐式算子

  • 显式矩阵 (Explicit Matrix)：
    • 存储方式：二维数组 float A[M][N] 或稀疏矩阵 scipy.sparse。
    • 优点：数学直观，可直接调用 inv(A) 或 eig(A)。
    • 缺点：对于 $512 \times 512$ 的图像，卷积矩阵大小约为 $2.6 \times 10^{11}$ 个元素（TB 级别内存），完全不可用。
  • 隐式算子 (Implicit Operator / Matrix-free)：
    • 存储方式：两个函数句柄（Function Handles）或一个类对象。
      • A(x)：实现前向操作（如 FFT $\to$ 乘核 $\to$ IFFT $\to$ 降采样）。
      • At(y)：实现伴随操作（如 上采样 $\to$ FFT $\to$ 乘共轭核 $\to$ IFFT）。
    • 优点：内存占用极低（仅需存储图像和核）。
    • 缺点：无法直接求逆，无法直接查看特征值，必须通过迭代法求解。

17.1.2 核心公理：点积测试 (The Dot-Product Test)

这是验证隐式算子是否正确的唯一真理。如果你编写了 A 和 At，但没有通过点积测试，后续的梯度下降方向将是错误的，FISTA/ADMM 等算法绝对不会收敛。

测试原理： 对于任意线性算子 $A: \mathbb{R}^N \to \mathbb{R}^M$，其伴随算子 $A^*: \mathbb{R}^M \to \mathbb{R}^N$ 必须满足： $$ \langle A(x), y \rangle_{\mathbb{R}^M} = \langle x, A^*(y) \rangle_{\mathbb{R}^N} $$ 标准工程实现流程：

  1. 随机化：生成随机向量 $x \in \mathbb{R}^N$ 和 $y \in \mathbb{R}^M$（通常服从标准正态分布）。
  2. 双向计算：
    • $v_1 = A(x)$
    • $v_2 = A^*(y)$
  3. 内积比对：
    • $val_1 = \text{sum}(v_1 \cdot y)$
    • $val_2 = \text{sum}(x \cdot v_2)$
  4. 相对误差判定： $$ \text{error} = \frac{|val_1 - val_2|}{\max(|val_1|, |val_2|) + \epsilon} $$ 该误差应小于 $10^{-12}$（双精度下）。

常见伴随算子推导表：

  • 卷积 $h * x$ $\longrightarrow$ 相关 $h(-\cdot) * y$ (实数域) 或 $\bar{h}(-\cdot) * y$ (复数域)。
  • 降采样 (Stride $k$) $\longrightarrow$ 零填充上采样 (Zero-filling with stride $k$)。
  • 掩膜 $M \odot x$ $\longrightarrow$ 掩膜 $M \odot y$ (对角矩阵的转置是自身)。
  • 梯度 $\nabla$ $\longrightarrow$ 负散度 $-\text{div}$ (注意边界条件处理)。

17.2 复杂度与加速：数值计算的艺术

17.2.1 FFT 的工程陷阱

在去模糊和 CSC（卷积稀疏编码）中，FFT 是核心加速器。

  1. 边界伪影 (Boundary Artifacts)： FFT 隐含周期性边界条件。直接卷积会导致图像右侧边缘的内容“穿越”到左侧。

    • 解决方案：edgetaper 或 padding。在进入 FFT 前，将图像 $H \times W$ 扩展为 $(H+K) \times (W+K)$，计算完后裁剪回原尺寸。 2. 最优尺寸 (Optimal Size)： FFT 算法在处理长度为 $2^k$ 或仅包含小质因数（2, 3, 5, 7）的长度时最快。

    • Bad: fft(img, 1021) (1021 是质数，退化为 $O(N^2)$)。

    • Good: fft(img, 1024)。

    • Rule-of-Thumb: 总是使用 next_fast_len() 函数来确定 padding 大小。 3. 中心化 (Centering)： 光学 PSF 通常以中心为原点，而 FFT 假设原点在数组左上角 $(0,0)$。

    • 必须操作：在做 fft2(kernel) 之前，必须执行 ifftshift(kernel)，否则卷积结果会发生半个图像宽度的平移。

17.2.2 矩阵求逆引理 (Sherman-Morrison-Woodbury)

在 ADMM 或半二次分裂（HQS）中，我们常遇到如下形式的子问题（岭回归形式）： $$ x_{k+1} = \arg\min_x \frac{\rho}{2} |x - z|^2 + \frac{\mu}{2} |Ax - y|^2 $$ 其解析解为： $$ (\rho I + \mu A^T A) x = \rho z + \mu A^T y $$

  • 直接法：直接求逆 $(\rho I + \mu A^T A)^{-1}$ 是灾难性的。
  • 频域法：如果 $A$ 是卷积，可利用 FFT 对角化求解。
  • Woodbury 加速：如果 $A$ 是“肥胖”矩阵（观测少，未知数多，如压缩感知 $M \ll N$），利用引理： $$ (\rho I + \mu A^T A)^{-1} = \frac{1}{\rho} \left( I - A^T (\frac{\rho}{\mu} I + A A^T)^{-1} A \right) $$ 此时只需对 $M \times M$ 的矩阵求逆，而非 $N \times N$。这在 MRI 重建和随机投影中能带来百倍加速。

17.2.3 步长的自动选择：幂法 (Power Method)

在使用梯度下降或 ISTA 时，步长 $\tau$ 必须满足 $\tau \le 1/L$，其中 $L = \lambda_{max}(A^T A)$ 是 Lipschitz 常数。

如何快速求 $L$？ 不要调用 eig(A'A)。使用幂法（Power Iteration）：

Initialize b_k as random vector
Loop 10-20 times:
    b_{k+1} = A^T(A(b_k))  // 应用 A^T A
    b_{k+1} = b_{k+1} / norm(b_{k+1})
L ≈ norm(A^T(A(b_{k+1})))

此过程只需运行一次，作为预计算步骤。

17.3 数值细节：精度与归一化

17.3.1 归一化与尺度不变性

很多初学者发现代码在自己的 demo 图片上有效，换一张图就失效。通常是尺度（Scale）问题。

  • 数据归一化： 始终将输入图像转换到 $[0, 1]$ 浮点域。不要在 $[0, 255]$ 上做优化。

    • 原因：$\lambda |x|_1$ 对数值大小敏感。如果像素值放大 255 倍，数据项 $|Ax-y|^2$ 会放大 $255^2 \approx 65000$ 倍，而正则项可能只放大 255 倍。你需要重新调参 $\lambda$。

    • 参数归一化： 编写算法时，最好让参数对图像尺寸不敏感。

    • Bad: loss = sum((Ax-y).^2) + lambda * sum(abs(x))

    • Good: loss = mean((Ax-y).^2) + lambda * mean(abs(x)) 这样图像分辨率变了，$\lambda$ 依然大约适用。

17.3.2 浮点陷阱

  • NaN 的传播：如果不慎出现了 0/0 或 inf，它会像病毒一样在迭代中污染所有像素。
    • 防御：在除法操作（如 x / max(norm, eps)）中，分母始终加上 eps（如 $10^{-10}$）。
  • 累积误差：在几千次迭代中，float32 的误差会累积。对于高精度要求的科学成像（如天文、医疗），核心迭代建议使用 float64，仅在显示/存储时转回 float32。

17.4 超参数调优与自适应策略

17.4.1 惩罚参数 $\lambda$ 的选取

  • 视觉调节法：先设 $\lambda$ 非常大（全平滑），然后指数级减小（$10, 1, 0.1, 0.01$），直到噪声刚开始出现。
  • Morozov 偏差原理 (Discrepancy Principle)： 如果已知噪声水平 $\sigma$，调整 $\lambda$ 使得最终残差 $|Ax - y|_2^2 \approx M \sigma^2$。 这可以通过牛顿法在外层循环自动寻找 $\lambda$。

17.4.2 ADMM 的自适应 $\rho$

ADMM 的惩罚参数 $\rho$ 极难手动调节。工业界标准做法是使用残差平衡（Residual Balancing）：

  • 原始残差 (Primal Residual) $r^k = Ax^k + Bz^k - c$（衡量约束满足程度）。
  • 对偶残差 (Dual Residual) $s^k = \rho A^T B (z^k - z^{k-1})$（衡量解的稳定性）。

策略：

  1. 如果 $|r^k| > \mu |s^k|$（原问题收敛太慢），则 $\rho \leftarrow \tau \rho$（加大惩罚）。
  2. 如果 $|s^k| > \mu |r^k|$（对偶问题收敛太慢），则 $\rho \leftarrow \rho / \tau$（减小惩罚）。
    • 典型值：$\mu=10, \tau=2$。

17.5 收敛诊断与可视化

不要盯着控制台滚动的数字看，你需要可视化监控。

17.5.1 监控仪表盘 (Dashboard)

一个好的复现脚本应实时记录以下曲线：

  1. 目标函数值 (Objective Value)：必须单调下降（对于凸问题）。如果震荡，步长太大。
  2. 对数坐标残差 (Log Residuals)：$\log_{10}|x^{k+1} - x^k|$。这能帮你判断是处于“线性收敛”还是“次线性收敛”。
  3. 度量指标 (PSNR/SSIM)：如果是在有 Ground Truth 的实验中，PSNR 通常先上升，过拟合后下降。这也是早停 (Early Stopping) 的依据。

17.5.2 假收敛 (False Convergence)

  • 现象：残差极小，代码停止，但结果是一团模糊。
  • 原因：正则化参数 $\lambda$ 太大，或者步长太小导致原地踏步。
  • 判断：检查对偶间隙 (Duality Gap)，或者简单地查看 $|Ax-y|$ 是否合理。

17.6 “跑不对”排查清单 (The Debug Playbook)

当复现结果完全错误时，请按此流程操作：

  1. Sanity Check (健全性检查)：

    • 设 $\lambda=0$，迭代几次。结果应该是含噪图像的最小二乘解（通常充满噪声）。如果结果是黑屏或平滑的，说明数据项（Data Fidelity）代码写错了。
    • 设 $y=0$。结果应该是 $0$。如果不是，说明代码中有硬编码的偏置。

  2. Gradient Check (梯度检验)： 如果你推导了复杂的正则项梯度（如非局部先验），不要盲信公式。 使用有限差分法验证： $$ \frac{f(x + \epsilon d) - f(x - \epsilon d)}{2\epsilon} \approx \langle \nabla f(x), d \rangle $$ 选取随机方向 $d$，计算数值微分与解析梯度的余弦相似度。应该接近 1.0。

  3. Overfitting Test (过拟合测试)： 取一张极小的图片（例如 $16 \times 16$），不加噪声。运行你的算法。

    • 结果必须完美恢复原图（机器精度级误差）。
    • 如果做不到，说明优化器本身有 bug，与统计特性无关。

  4. Identity Test (恒等测试)： 将算子 $A$ 设为单位矩阵 $I$。此时去噪算法应该退化为简单的软阈值（Soft Thresholding）或 Prox 算子。验证结果是否符合预期。

17.7 本章小结

  • 实现：优先使用隐式算子 + 点积测试，这是大规模优化的基础。
  • 加速：善用 FFT（注意 padding 和 shift）与 Woodbury 矩阵恒等式。
  • 调参：使用 Power Method 自动定步长，使用残差平衡自动定 $\rho$。
  • 心态：优化代码的调试不同于逻辑代码。可视化曲线和梯度检验是你的听诊器。

17.8 练习题

基础题

  1. 点积测试实践： 编写一个 Python/MATLAB 函数，输入为两个函数句柄 A 和 At，输出一个布尔值，判断它们是否互为转置。（需考虑复数情况）。

    Hint注意复数域内积定义 $\langle u, v \rangle = u^H v$。测试应包含实部和虚部。

  2. Lipschitz 常数估计： 构造一个 $A$ 为随机矩阵。分别使用 scipy.linalg.svdvals (求最大奇异值) 和本章介绍的幂法计算 Lipschitz 常数。比较两者的计算时间与结果精度。

    Hint当矩阵尺寸 $>5000$ 时，SVD 会极慢，而幂法依然瞬间完成。

  3. FFT 移位验证： 生成一个中心为白点的 $11 \times 11$ 图像作为核。计算其 FFT。观察 fft2(kernel) 和 fft2(ifftshift(kernel)) 的相位谱（Phase Spectrum）区别。

    Hint前者会有剧烈的线性相位变化（条纹），后者相位接近零（平坦）。线性相位对应空间位移。

挑战题

  1. 调试 PnP-ADMM： 假设你实现了一个 Plug-and-Play ADMM，使用预训练的 DnCNN 作为去噪器。但在迭代中，图像并没有变清晰，而是逐渐变暗并最终全黑。请列出至少 3 种可能的原因和排查手段。

    Hint1. DnCNN 输入未归一化（期望 [0,1] 但输入了 [0,255]）。2. ADMM 步长/参数不匹配，导致去噪强度过大。3. DnCNN 训练时减去了均值，而应用时未处理。

  2. 内存极限挑战： 在只有 4GB 显存的 GPU 上，如何实现一个 3D 视频（$T=100, H=512, W=512$）的卷积稀疏编码？请设计分块（Patch-based）或在线（Online）处理策略。

    Hint不能一次性将变量载入。可以使用“梯度累积”思路，或者将视频沿时间轴切片，利用 ADMM 的可分裂性逐块更新。

17.9 常见陷阱与错误 (Gotchas)

  • Python 的 copy 问题： 在迭代 $x_{k+1} = x_k - \dots$ 时，简单的赋值 x_old = x_new 在 Python 中只是引用拷贝。修改 x_new 会同时改变 x_old，导致动量项（Momentum）计算错误。

    • Fix: 使用 x_old = x_new.copy()。
    • OpenCV vs. Scikit-Image vs. MATLAB：
    • 读取图片：OpenCV 是 BGR，其他是 RGB。
    • 范围：OpenCV 读入是 uint8 [0,255]，skimage.img_as_float 是 [0,1]。
    • 混用库是复现灾难的源头。
    • 一维向量 vs. 二维矩阵： numpy 中 (N,) 数组和 (N, 1) 数组在广播（Broadcasting）时行为不同。点积测试中常因此翻车。建议始终显式保持维度 keepdims=True。