4. 正则化方法入门

反卷积是反问题中最典型的病态问题之一。设一维卷积模型为 y = h * x + n，其中 h 是已知的核（点扩散函数），x 是待恢复的真实信号，n 是噪声。直接对卷积矩阵求逆会把噪声放大成剧烈振荡，解完全不可用。正则化的作用，是在“贴合数据”和“保持解的平滑/稳定”之间引入一个可控的折中，而正则化参数的选择，正是决定这个折中点落在哪里的关键。

一、Tikhonov正则化的形式与直觉

Tikhonov正则化把原问题改写为一个带惩罚项的最小二乘问题：

min_x  ||H x - y||²  +  λ ||L x||²

第一项是数据残差，衡量解与观测的吻合程度；第二项是正则项，L 通常取单位矩阵（零阶，抑制幅值）或一阶差分矩阵（抑制剧烈跳变）。λ > 0 是正则化参数。对上述目标求导并令其为零，可得闭式解：

x_λ = (HᵀH + λ LᵀL)⁻¹ Hᵀ y

理解 λ 的作用，可以看两个极端：λ → 0 时解退化为最小二乘解，残差最小但噪声被放大，解出现高频振荡；λ → ∞ 时正则项主导，解被强烈压向零（或压向常数，取决于 L），残差很大，信号细节被抹平。因此 λ 是一个“偏差—方差”权衡旋钮：增大它降低方差、提高稳定性，同时引入偏差、损失分辨率。

二、L曲线：用几何形状选参

L曲线法把不同 λ 下的解残差范数 ||H x_λ - y|| 与正则解范数 ||L x_λ|| 分别作为横纵坐标，在双对数坐标下描点。随着 λ 从很小增大到很大，曲线通常呈现明显的“L”形：竖直段对应 λ 太小、残差小但解范数急剧增大（噪声放大区）；水平段对应 λ 太大、解范数小但残差迅速上升（过平滑区）。拐角处是两者相对平衡的位置，通常取曲率最大点作为推荐参数。

L曲线的优点是直观、不需要噪声方差先验；缺点是拐角有时并不尖锐，曲率计算对采样密度和数值差分敏感，需要在对数域上做平滑后再求曲率，否则容易选到抖动造成的伪拐点。

三、广义交叉验证（GCV）：用预测能力选参

GCV 的思想是“留一法”的解析近似：对每个 λ，考察去掉第 i 个观测后模型对它的预测误差，再对所有 i 平均。利用正则解的线性性质，可以避免真正做 N 次重解，得到闭式指标：

GCV(λ) = ||H x_λ - y||² / (N - trace(A_λ))²
A_λ = H (HᵀH + λ LᵀL)⁻¹ Hᵀ

其中 A_λ 称为影响矩阵（帽子矩阵），trace(A_λ) 是有效自由度，反映解对数据的拟合能力。分母中的 (N - trace(A_λ))² 是对留一法误差的修正：λ 越小，模型越“灵活”，trace(A_λ) 越接近 N，分母变小，惩罚加重，从而抑制过拟合。取使 GCV(λ) 最小的 λ 即为推荐值。

GCV 的优点是无需人工干预、有统计解释；缺点是在某些问题上曲线可能很平坦甚至出现多个局部极小，且需要计算 trace(A_λ)，对大规模问题需借助随机迹估计或利用矩阵结构加速。

四、实现一个简单的一维反卷积求解

下面给出一个可直接运行的 Python 示例，用高斯核构造卷积矩阵，对含噪观测做 Tikhonov 正则化，并同时用 L 曲线与 GCV 选参。示例参数为教学假设：信号长度 N=200，高斯核标准差 σ=3，噪声标准差约为信号幅值的 2%。

import numpy as np

N, sigma, lam_grid = 200, 3.0, np.logspace(-6, 1, 60)

# 1) 构造高斯核与卷积矩阵 H（Toeplitz）
k = np.arange(-15, 16)
h = np.exp(-k**2 / (2*sigma**2)); h /= h.sum()
H = np.zeros((N, N))
for i in range(N):
    for j, kj in enumerate(k):
        col = i + kj
        if 0 <= col < N:
            H[i, col] = h[j]

# 2) 构造真实信号与含噪观测（示例数据）
x_true = np.zeros(N)
x_true[60:90] = 1.0; x_true[120:150] = 0.7
y = H @ x_true + 0.02 * np.random.randn(N)

# 3) 一阶差分正则矩阵 L
L = np.diff(np.eye(N), axis=0)

# 4) 扫描 λ，记录残差、解范数与 GCV
res, reg, gcv = [], [], []
for lam in lam_grid:
    A = H.T @ H + lam * (L.T @ L)
    x = np.linalg.solve(A, H.T @ y)
    r = H @ x - y
    res.append(np.linalg.norm(r))
    reg.append(np.linalg.norm(L @ x))
    # 影响矩阵迹：trace(H A⁻¹ Hᵀ)
    tr = np.trace(H @ np.linalg.solve(A, H.T))
    gcv.append(np.sum(r**2) / (N - tr)**2)

res, reg, gcv = map(np.array, (res, reg, gcv))
lam_lcurve = lam_grid[np.argmax(np.gradient(np.gradient(np.log(reg))) /
                                np.gradient(np.log(res)))]  # 简化曲率近似
lam_gcv = lam_grid[np.argmin(gcv)]
print("L曲线 λ*:", lam_lcurve, " GCV λ*:", lam_gcv)

运行后应观察三件事：其一，res 随 λ 增大单调上升，reg 单调下降，二者在对数坐标下构成 L 形；其二，gcv 曲线通常先降后升，最小值位置给出推荐 λ；其三，用两个推荐 λ 分别重建 x，与 x_true 对比，应看到 λ 过小时解在平顶区出现高频毛刺，λ 过大时两个矩形脉冲的边沿被明显圆滑、幅值被压低。若两种准则给出的 λ 相差一个数量级以上，说明问题对参数较敏感，应结合残差是否接近噪声水平（本例约 0.02·√N ≈ 0.28）做交叉判断，而不是盲信单一准则。

五、选参时的实用判断

无论用哪种准则，都应回到“残差是否与噪声水平匹配”这一物理判据上做校验。若最优 λ 对应的残差范数远大于估计的噪声范数，说明正则化过强、信号被过度压制；若远小于噪声范数，说明仍在拟合噪声。此外，L 的选择会改变 L 曲线的形状与 GCV 的最优点，零阶正则倾向于把解压向零、适合稀疏小幅值信号，一阶正则倾向于平滑、适合连续缓变信号。实践中可先固定 L 扫描 λ，再换 L 比较重建结果与残差，形成“正则算子 + 正则参数”的联合选择，而不是只调一个标量。

把这一节的方法落到一维反卷积上，核心就是三步：写出 H 与 L，在 λ 网格上批量求解并记录残差、解范数与 GCV 指标，最后用 L 曲线拐角或 GCV 最小值给出推荐参数，并用噪声水平做合理性校验。掌握这套流程后，二维图像复原、谱反卷积等问题只是把矩阵结构换掉，选参逻辑完全一致。