在现代图像处理（Image Processing）和计算机视觉（Computer Vision）的文献中，你经常会看到精简优雅的公式，例如 $\hat{x} = \arg\min_x \frac{1}{2}|Ax-y|^2 + \lambda |x|_{TV}$。然而，当试图将这些公式转化为代码时，巨大的鸿沟便会出现：
本章是全书的“工具箱”。我们不追求数学系般的严谨证明，而是侧重于建立“算子视角”（Operator Perspective）。
学习目标：
在变分法推导中，为了使用线性代数工具（如 $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\]这种表示法虽然在推导时方便，但在工程上带来了两个假象：
这是从教科书到实战的最重要跨越。
(N, N) 的数组。仅适用于极小规模问题（如 $32 \times 32$ 图像）。A(x) or fwd): 计算 $y = Ax$。At(y) or adj): 计算 $\hat{x} = A^T y$。Rule of Thumb 2.1: 在编写优化求解器（Solver）时，应该设计一个接受
input和mode（前向/伴随）的通用接口，或者使用面向对象的LinearOperator类。永远不要尝试对全尺寸图像构建密集矩阵。
对于去模糊（Deblurring）问题，观测模型为 $y = k * x + n$。
在一维情况下，离散卷积等价于一个 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) 结构。
边界处理不仅影响视觉效果，还决定了矩阵 $A$ 的性质和能否快速求逆。
| 边界条件 (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) | 一般稀疏矩阵 | 无 | 迭代求解 |
在优化中，我们最常用循环边界条件，因为根据卷积定理，循环矩阵可以被傅里叶基对角化： \(K = F^H \cdot \text{diag}(\hat{k}) \cdot F\) 这意味着 $Kx$ 可以通过 FFT 计算：
全变分（Total Variation, TV）去噪的核心在于对图像梯度的稀疏约束。准确理解梯度的离散化及其伴随是实现 TV 求解器的关键。
定义一维离散梯度算子 $D$ 为前向差分： \((Dx)_i = x_{i+1} - x_i, \quad i=1,\dots,N-1\) 边界点 $N$ 通常定义为 $0$ (Neumann 边界) 或 $x_1 - x_N$ (周期边界)。
这是初学者最容易出错的地方。我们需要找到 $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}\)
对于图像 $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}$，你的算法不可能收敛。
优化目标的形式决定了解的性质。
| $\ell_1$ (Manhattan): $|x|_1 = \sum | x_i | $。 |
当我们处理彩色图像或组稀疏时，会用到 $\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 }\)
\(\|X\|_* = \sum \sigma_i(X)\)
当处理视频（3D）或高光谱（3D/4D）数据时，简单的向量化会破坏结构。
对于可分离算子（如 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)$。 这意味着我们不需要构建大矩阵，只需对图像进行左右矩阵乘法。
对于张量 $\mathcal{X} \in \mathbb{R}^{n_1 \times n_2 \times n_3}$，模-1 展开 $X_{(1)}$ 是将其按切片并排成一个宽矩阵。
这是复现论文时的速查表。
| 任务 | 正向模型 $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$ | 高斯或伯努利矩阵 |
对于低秩问题，核心算子不是卷积，而是 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)$。
fftshift 的滥用与遗忘psf2otf 或 ifftshift，导致结果整体平移了半张图。A' 或 A.conj().T 时，务必确保处理了虚部。如果只做转置不做共轭，梯度方向会出错。ifft2(fft2(x)) 精确等于 x。有些库需要手动乘除归一化因子。检查你的算子库是否满足 $A^* A \approx I$ （如果是正交变换）或能量守恒。