第2章 数学与算子预备：线性代数、离散化与成像算子

2.1 开篇与学习目标

在现代图像处理（Image Processing）和计算机视觉（Computer Vision）的文献中，你经常会看到精简优雅的公式，例如 $\hat{x} = \arg\min_x \frac{1}{2}|Ax-y|^2 + \lambda |x|_{TV}$。然而，当试图将这些公式转化为代码时，巨大的鸿沟便会出现：

  1. 维度的诅咒：图像 $x$ 在纸上是二维矩阵，在公式里往往变成了列向量，在内存中又是多维数组。
  2. 算子的黑盒：$A$ 到底是什么？对于去模糊，它是卷积；对于 MRI，它是傅里叶变换。我们如何用统一的代数语言描述它们？
  3. 伴随的缺失：几乎所有的优化算法（梯度下降、ADMM、Primal-Dual）都要求我们计算 $A^T$（或 $A^*$）。如果不理解其物理含义，就无法写出正确的反向投影代码。

本章是全书的“工具箱”。我们不追求数学系般的严谨证明，而是侧重于建立“算子视角”（Operator Perspective）。

学习目标：

  1. 思维转换：从“像素操作”转向“向量空间中的算子操作”。
  2. 隐式算子：学会处理无法显式存储的大型矩阵，理解 Linear Map 的编程范式。
  3. 离散微分：深入理解梯度、散度、拉普拉斯算子的离散形式及其伴随关系（Adjointness）。
  4. 几何直觉：掌握 $\ell_1, \ell_2, \ell_{2,1}, |\cdot|_*$ 等范数在优化中的几何行为（稀疏、平滑、群稀疏、低秩）。
  5. 张量初步：理解高维数据（视频、高光谱）的矩阵化（Unfolding）表示。

2.2 向量化、矩阵化与“把图像当向量”的代价

2.2.1 字典序堆叠（Lexicographic Ordering）

在变分法推导中，为了使用线性代数工具（如 $Ax=b$），我们通常假设图像 $X \in \mathbb{R}^{H \times W}$ 被拉直为一个列向量 $x \in \mathbb{R}^{N}$（其中 $N=HW$）。

$$ \text{vec}(X) = x = [X_{1,1}, \dots, X_{H,1}, X_{1,2}, \dots, X_{H,W}]^T $$ 这种表示法虽然在推导时方便，但在工程上带来了两个假象：

  1. 邻域丢失假象：向量 $x$ 中，元素 $x_i$ 和 $x_{i+1}$ 可能在原图中并不相邻（例如跨行时）。因此，定义空间正则项（如 TV）时，必须回到 2D 索引去理解。
  2. 存储爆炸假象：如果 $x$ 是 $10^6$ 维向量，那么任何线性变换 $A$ 理论上都是 $10^6 \times 10^6$ 的矩阵。直接存储 $A$ 需要数 TB 内存。

2.2.2 显式矩阵 vs. 隐式算子 (Explicit vs. Implicit)

这是从教科书到实战的最重要跨越。

  • 显式矩阵 (Explicit Matrix): 真的在内存中构建一个 (N, N) 的数组。仅适用于极小规模问题（如 $32 \times 32$ 图像）。
  • 隐式算子 (Implicit Operator / Linear Map): 将 $A$ 视为一个函数。我们只需要实现两个操作：
    • Forward mode (A(x) or fwd): 计算 $y = Ax$。
    • Adjoint mode (At(y) or adj): 计算 $\hat{x} = A^T y$。

Rule of Thumb 2.1: 在编写优化求解器（Solver）时，应该设计一个接受 input 和 mode（前向/伴随）的通用接口，或者使用面向对象的 LinearOperator 类。永远不要尝试对全尺寸图像构建密集矩阵。

2.3 卷积、Toeplitz 与频域加速

对于去模糊（Deblurring）问题，观测模型为 $y = k * x + n$。

2.3.1 卷积的矩阵形式：Toeplitz 矩阵

在一维情况下，离散卷积等价于一个 Toeplitz（托普利茨） 矩阵与向量相乘。Toeplitz 矩阵的特点是对角线元素恒定。

向量 x: [x1, x2, x3, x4, x5]
卷积核 k: [k_c, k_r] (Example: [0.5, 0.5])

对应矩阵 A (Zero Boundary Condition):
[ k_c  k_r   0    0    0  ]
[  0   k_c  k_r   0    0  ]
[  0    0   k_c  k_r   0  ]
[  0    0    0   k_c  k_r ]
[  0    0    0    0   k_c ]

在二维情况下，矩阵 $A$ 具有 BTTB (Block Toeplitz with Toeplitz Blocks) 结构。

2.3.2 边界条件与矩阵结构的对应

边界处理不仅影响视觉效果，还决定了矩阵 $A$ 的性质和能否快速求逆。

| 边界条件 (Boundary) | 矩阵结构 (Matrix Structure) | 对角化工具 (Diagonalization) | 计算复杂度 |

边界条件 (Boundary) 矩阵结构 (Matrix Structure) 对角化工具 (Diagonalization) 计算复杂度
零填充 (Zero) Toeplitz 无简单对角化 $O(N \log N)$ (需补零FFT)
周期/循环 (Circular) 循环矩阵 (Circulant) FFT $O(N \log N)$
反射/镜像 (Reflective) Toeplitz + Hankel DCT (离散余弦变换) $O(N \log N)$
复制 (Replicate) 一般稀疏矩阵 无 迭代求解

2.3.3 FFT 与光学传递函数 (OTF)

在优化中，我们最常用循环边界条件，因为根据卷积定理，循环矩阵可以被傅里叶基对角化： $$ K = F^H \cdot \text{diag}(\hat{k}) \cdot F $$ 这意味着 $Kx$ 可以通过 FFT 计算：

  1. 计算核的 FFT：$\hat{k} = \mathcal{F}(k)$。这里的 $\hat{k}$ 被称为 光学传递函数 (OTF)。
    • 注意：$k$ 通常很小（如 15x15），必须补零（pad）到与 $x$ 相同大小再做 FFT。
  2. 计算 $Ax = \mathcal{F}^{-1}(\hat{k} \odot \mathcal{F}(x))$。
  3. 计算 $A^T y = \mathcal{F}^{-1}(\overline{\hat{k}} \odot \mathcal{F}(y))$。注意这里是对 OTF 取复共轭 (Complex Conjugate)。

2.4 离散微分算子：全变分的核心

全变分（Total Variation, TV）去噪的核心在于对图像梯度的稀疏约束。准确理解梯度的离散化及其伴随是实现 TV 求解器的关键。

2.4.1 前向差分 (Forward Difference) $\nabla$

定义一维离散梯度算子 $D$ 为前向差分： $$ (Dx)_i = x_{i+1} - x_i, \quad i=1,\dots,N-1 $$ 边界点 $N$ 通常定义为 $0$ (Neumann 边界) 或 $x_1 - x_N$ (周期边界)。

2.4.2 伴随算子：负散度 (Negative Divergence) $\nabla^*$

这是初学者最容易出错的地方。我们需要找到 $D^T$ 使得 $\langle Dx, y \rangle = \langle x, D^T y \rangle$。 利用“分部求和”（Summation by parts）： $$ \sum_{i=1}^{N-1} (x_{i+1}-x_i)y_i = \sum_{i=1}^{N-1} x_{i+1}y_i - \sum_{i=1}^{N-1} x_i y_i $$ 令 $j=i+1$，第一项变为 $\sum_{j=2}^N x_j y_{j-1}$。整理各项系数，你会发现 $x_i$ 的系数是 $y_{i-1} - y_i$。

结论：前向差分的伴随是后向差分的负数。 在物理上，梯度的伴随是负散度： $$ \nabla^* = -\text{div} $$

2.4.3 2D 梯度与拉普拉斯

对于图像 $X$，梯度算子 $\nabla$ 输出两个分量（水平和垂直）： $$ \nabla X = (\mathbf{g}_x, \mathbf{g}_y) \in \mathbb{R}^{H \times W \times 2} $$ 其中 $\mathbf{g}_x = D_h X, \mathbf{g}_y = D_v X$。

  • 伴随算子（输入是两个矩阵 $p, q$）： $$ \nabla^*(p, q) = D_h^T p + D_v^T q = -\text{div}(p, q) $$ 它将向量场映射回标量图像。

  • 拉普拉斯算子： $$ \Delta X = \text{div}(\nabla X) = -\nabla^* (\nabla X) $$ 对应的离散掩膜通常是 $\begin{bmatrix} 0 & 1 & 0 \\ 1 & -4 & 1 \\ 0 & 1 & 0 \end{bmatrix}$。

Rule of Thumb 2.2: 调试时，必须验证你的 grad 和 div 函数满足点积测试（Dot Product Test）：sum(grad(x) * p) ≈ sum(x * div(p))（注意 div 定义的正负号）。如果误差大于 $10^{-10}$，你的算法不可能收敛。

2.5 范数王国：从几何形状看先验

优化目标的形式决定了解的性质。

2.5.1 $\ell_p$ 范数体系

  • $\ell_2$ (Euclidean): $|x|_2 = \sqrt{\sum x_i^2}$。
    • 导数：线性 ($x$)。
    • 几何：球形。能量分散，倾向于非零的稠密解。
    • 应用：高斯噪声似然、吉洪诺夫正则（平滑）。
  • $\ell_1$ (Manhattan): $|x|_1 = \sum |x_i|$。
    • 导数：$\text{sign}(x)$。
    • 几何：菱形（多胞体）。在坐标轴处有尖点，极易将解推向 0。
    • 应用：稀疏编码、脉冲噪声去除、压缩感知。
  • $\ell_0$ (Pseudo-norm): $|x|_0 = #\{i : x_i \neq 0\}$。
    • 几何：直接定义在坐标轴上。
    • 应用：理想的稀疏性，但非凸难解（通常松弛为 $\ell_1$）。

2.5.2 混合范数 (Mixed Norms) / 结构化稀疏

当我们处理彩色图像或组稀疏时，会用到 $\ell_{2,1}$ 范数。 假设图像梯度在某点 $(i,j)$ 为向量 $g_{i,j} \in \mathbb{R}^2$。 $$ | \nabla X |_{2,1} = \sum_{i,j} \sqrt{ (D_h X)_{i,j}^2 + (D_v X)_{i,j}^2 } $$

  • 内层 $\ell_2$：耦合水平和垂直方向（各向同性 TV），希望梯度的大小作为一个整体被惩罚，而不是分量单独惩罚。
  • 外层 $\ell_1$：希望大部分位置的梯度模长为 0（分段常数）。

2.5.3 核范数 (Nuclear Norm)

$$ |X|_* = \sum \sigma_i(X) $$

  • 它是矩阵奇异值的 $\ell_1$ 范数。
  • 这是秩函数 $\text{rank}(X)$ 的凸松弛。用于低秩矩阵恢复（Robust PCA）。

2.6 张量与 Kronecker 积：高维数据的处理

当处理视频（3D）或高光谱（3D/4D）数据时，简单的向量化会破坏结构。

2.6.1 Kronecker 积 (Kronecker Product)

对于可分离算子（如 2D 高斯模糊），其大矩阵 $A$ 可以写成两个小矩阵的 Kronecker 积： $$ A = A_c \otimes A_r $$ 其中 $A_r$ 作用于行，$A_c$ 作用于列。 性质：$(A_c \otimes A_r) \text{vec}(X) = \text{vec}(A_r X A_c^T)$。 这意味着我们不需要构建大矩阵，只需对图像进行左右矩阵乘法。

2.6.2 模-n 展开 (Mode-n Unfolding)

对于张量 $\mathcal{X} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$，模-1 展开 $X_{(1)}$ 是将其按切片并排成一个宽矩阵。

  • 在视频处理中，将视频帧拉成列向量组成的矩阵，此时矩阵的低秩性对应于背景静止或线性相关。

2.7 常用成像算子库 (The Operator Zoo)

这是复现论文时的速查表。

| 任务 | 正向模型 $y = A(x)$ | 算子 $A$ 描述 | 伴随算子 $A^*(y)$ | 备注 |

任务 正向模型 $y = A(x)$ 算子 $A$ 描述 伴随算子 $A^*(y)$ 备注
去噪 $y = x + n$ Identity ($I$) Identity ($I$) 最简单
去模糊 $y = k * x$ Convolution ($K$) Correlation ($K^\dagger$) $A^*$ 卷积核翻转
下采样 $y = x_{\downarrow s}$ Subsampling ($S$) Upsampling w/ zeros ($S^T$) $A^*$ 补零而非插值
修复 (Inpainting) $y = M \odot x$ Masking ($M$) Masking ($M$) $M$ 为 0/1 矩阵，自伴随
MRI $y = P \mathcal{F} x$ Fourier + Mask $\mathcal{F}^{-1} P^T y$ 复数域共轭
CT (Radon) $y = \mathcal{R} x$ Radon Transform Back-projection 也就是反投影
压缩感知 $y = \Phi x$ Random Matrix $\Phi^T y$ 高斯或伯努利矩阵

2.8 SVD 与奇异值阈值化 (SVT)

对于低秩问题，核心算子不是卷积，而是 SVD。 $$ X = U \Sigma V^T, \quad \Sigma = \text{diag}(\sigma_i) $$ 奇异值阈值化 (SVT) 算子 $\mathcal{D}_\tau(X)$ 是核范数正则化问题的近端算子（Proximal Operator）： $$ \mathcal{D}_\tau(X) = U \cdot \text{soft}(\Sigma, \tau) \cdot V^T $$ 其中 $\text{soft}(x, \tau) = \text{sign}(x)\max(|x|-\tau, 0)$。

  • 直觉：SVD 分解出图像的主要成分（大奇异值）和细节/噪声（小奇异值）。SVT 砍掉小奇异值，强制去噪或提取低秩背景。

2.9 本章小结

  1. 代码即算子：抛弃 $A$ 是矩阵的观念，在代码中 $A$ 是一个函数/类。
  2. 伴随是关键：$A^*$ 类似于转置，对于卷积是翻转核，对于梯度是负散度，对于下采样是补零。
  3. 范数定乾坤：$\ell_1$ 导致稀疏，$\ell_2$ 导致平滑，$|\cdot|_*$ 导致低秩。
  4. 边界需谨慎：FFT 对应循环边界，如果不匹配会产生振铃。

2.10 练习题

习题 2.1（基础）：1D 卷积与矩阵

给定信号 $x = [1, 2, 3, 4]^T$，卷积核 $k=[1, -1]$（$k$ 的原点在第一个元素）。

  1. 写出对应“零填充边界”的卷积矩阵 $A$（输出大小 4）。
  2. 写出 $A^T$。
  3. 计算 $y=Ax$ 和 $\hat{x}=A^T y$。你会发现 $A^T$ 在做什么操作？

Hint: $A$ 是下双对角矩阵。$A^T$ 实际上是用翻转后的 $k$ 进行卷积（相关）。

习题 2.2（基础）：FFT 卷积的中心化

在 Matlab/Python 中，fft2 假设数据的原点在 (0,0)（左上角）。如果你有一个 $3 \times 3$ 的高斯核，中心在中间。直接做 fft2 会发生什么？你应该如何预处理这个核？

Hint: 思考相位平移。如果不做 ifftshift，卷积结果会发生位移。

习题 2.3（挑战）：下采样算子的伴随

定义下采样算子 $S$ 为：保留偶数索引像素，丢弃奇数索引像素。 证明：$S$ 的伴随算子 $S^T$ 是将低分辨率图像放回偶数位置，并在奇数位置填 0（而不是线性插值）。 这是一个常见的误解，很多初学者在实现超分辨率时，错误的把 $S^T$ 写成了双三次插值。

Hint: 设 $y = Sx$，写出内积 $\langle Sx, y \rangle = \sum y_k (Sx)_k$，还原回 $\sum x_i (\dots)$ 的形式。

习题 2.4（挑战）：TV 的伴随推导

对于 2D 图像，证明 $\langle \nabla u, \mathbf{p} \rangle = \langle u, -\text{div } \mathbf{p} \rangle$。你需要明确写出边界项在何种边界条件下会消失（例如周期边界或 Neumann 边界）。

Hint: 使用二维的分部求和公式，关注边界项的索引。

习题 2.5（思考）：$\ell\_0$ 与 $\ell\_1$ 的几何

画出二维平面上，方程 $x_1 + x_2 = 1$（直线）与 $\ell_1$ 球（菱形）和 $\ell_2$ 球（圆）的相切情况。解释为什么 $\ell_1$ 倾向于在坐标轴上相切（得到稀疏解），而 $\ell_2$ 倾向于在非零处相切。

2.11 常见陷阱 (Gotchas)

1. fftshift 的滥用与遗忘

  • 错误：在做频域导数或卷积时，忘记 psf2otf 或 ifftshift，导致结果整体平移了半张图。
  • 调试：输入一个脉冲图像（中间为1，其余为0），看经过 $A$ 和 $A^*$ 后脉冲的位置是否正确。

2. 伴随算子 $\neq$ 逆算子

  • 概念：$A^*$ 只是转置。$A^{-1}$ 是逆。
  • 后果：如果你在迭代公式中需要 $A^{-1}$（例如去卷积的直接解），却用了 $A^*$（匹配滤波），图像还是模糊的。但在梯度下降 $\nabla f(x) = A^T(Ax-y)$ 中，我们需要的是 $A^T$。不要混淆求解器需要的算子类型。

3. 复数域的共轭

  • 在处理 MRI 或频域数据时，内积定义为 $\langle x, y \rangle = x^H y = \sum x_i \overline{y_i}$。
  • 陷阱：在 Python/Matlab 中计算 A' 或 A.conj().T 时，务必确保处理了虚部。如果只做转置不做共轭，梯度方向会出错。

4. 归一化因子的丢失

  • FFT 和 IFFT 经常伴随着 $1/N$ 或 $1/\sqrt{N}$ 的因子。
  • Rule：确保 ifft2(fft2(x)) 精确等于 x。有些库需要手动乘除归一化因子。检查你的算子库是否满足 $A^* A \approx I$ （如果是正交变换）或能量守恒。