线性问题

所谓线性问题，是指变量之间满足线性关系的数学模型，一般可以表示为：

y=Ax+ω(1.1)y = A x + \omega \qquad (1.1)

其中，AA 是已知的系数矩阵，维度通常为 m×nm \times n，xx 是待求解的未知向量；yy 是已知的观测向量；ω\omega 是观测噪声。 方程组的本质是要求在给定 yy 与 AA 的条件下，找到最符合关系的 xx。

在线性回归中，AA 由输入特征构成，xx 是回归系数，yy 是响应值或观测输出。

正常情况下

当 AA 满秩且条件良好时，最常见的求解方法是最小二乘法。其思想是最小化误差平方和：

min⁡x∥Ax−y∥22(1.2)\min_x \| A x - y \|_2^2 \qquad (1.2)

对上式求导并令其等于零，可以得到封闭解：

x=(ATA)−1ATy(1.3)x = (A^T A)^{-1} A^T y \qquad (1.3)

此时，若 AA 条件数较小，求得的 xx 稳定且精度高。

非正常情况下

在实际问题中，AA 可能是病态矩阵，即其条件数非常大。这种情况下，最小二乘解对观测噪声极其敏感，微小的扰动都会导致解 xx 出现巨大偏差。

为了解决这个问题，人们引入了正则化（Regularization）思想，即在优化目标中加入约束项，以限制解的复杂度并提高稳定性。最早的正则化方法是Tikhonov 正则化，由苏联数学家Andrey Tikhonov在1943年提出。这一方法在信号处理、计算机视觉、统计学习等领域都有广泛应用。

正则化方法

L2正则化

L2正则化，也称为岭回归（Ridge Regression），其优化目标是在最小二乘误差基础上加入二范数惩罚项：

min⁡x∥Ax−y∥22+λ∥x∥22(2.1)\min_x \| A x - y \|_2^2 + \lambda \| x \|_2^2 \qquad (2.1)

其中 λ>0\lambda > 0 是正则化系数。对该目标求导并令其等于零，可以得到闭式解：

x=(ATA+λI)−1ATy(2.2)x = (A^T A + \lambda I)^{-1} A^T y \qquad (2.2)

可以看到，相较于最小二乘解，L2正则化解在矩阵中多了 λI\lambda I 项，这引入了偏差，但能够显著降低解的方差，从而提高数值稳定性。

L1正则化

L1正则化，也称为Lasso 回归（Least Absolute Shrinkage and Selection Operator），其优化目标在误差平方和基础上加入一范数惩罚：

min⁡x∥Ax−y∥22+λ∥x∥1(2.3)\min_x \| A x - y \|_2^2 + \lambda \| x \|_1 \qquad (2.3)

与L2正则化不同，L1正则化无法得到简单的闭式解，需要使用迭代算法求解，如内点法。L1范数的特点是能够产生稀疏解，即很多系数被压缩为零，这正好可以用于处理压缩感知问题，使得在观测数据有限的情况下依然可以有效恢复信号。

L2与L1正则化
图 1：L2正则化和L1正则化的物理意义
左图：L1正则化对应解向量的菱形约束（产生稀疏解），右图：L2正则化对应解向量的圆形约束（均匀收缩）

迭代收缩软阈值算法

迭代收缩阈值算法（ISTA, Iterative Shrinkage-Thresholding Algorithm）是一种用于求解稀疏线性问题的迭代优化方法。该算法特别适用于带有 L1 正则化的优化问题，例如 LASSO 回归和压缩感知问题。

ISTA 的基本思想是通过迭代不断逼近稀疏解，每次迭代由两步组成：梯度下降和软阈值收缩。其迭代格式可表示为：

xk=Sλt(xk−1−2tAT(Axk−1−y))(3.1)x^k = S_{\lambda t} \Big( x^{k-1} - 2 t A^T (A x^{k-1} - y) \Big) \qquad (3.1)

其中，xk−1x^{k-1} 是第 k−1k-1 次迭代的解向量，tt 是步长（学习率），通常与 ATAA^T A 的最大特征值有关，Sλt(⋅)S_{\lambda t}(\cdot) 是软阈值函数，用于引入稀疏性，λ\lambda 表示正则化参数，控制稀疏程度。

通过不断迭代，xkx^k 会收敛到 L1 正则化优化问题的稀疏解，从而实现信号恢复或稀疏表示。

公式推导

我们先考察一个无正则项的最优化问题，其目标函数为：

min⁡xf(x)=∥Ax−y∥22(4.1)\min_x f(x) = \| A x - y \|_2^2 \qquad (4.1)

梯度下降法

设第 kk 次迭代的解为 xk−1x^{k-1}，根据梯度下降法的迭代格式，有：

rk=xk−1−t∇f(xk−1)(4.2)r^k = x^{k-1} - t \nabla f(x^{k-1}) \qquad (4.2)

其中，t>0t > 0 为步长，∇f(xk−1)\nabla f(x^{k-1}) 是在 xk−1x^{k-1} 处的梯度。

二阶泰勒展开

在点 xk−1x^{k-1} 附近，对 f(x)f(x) 进行二阶泰勒展开近似：

f(x)≈f(xk−1)+∇f(xk−1)T(x−xk−1)+12(x−xk−1)T∇2f(xk−1)(x−xk−1)(4.3)f(x) \approx f(x^{k-1}) + \nabla f(x^{k-1})^T (x - x^{k-1}) + \frac{1}{2}(x - x^{k-1})^T \nabla^2 f(x^{k-1}) (x - x^{k-1}) \qquad (4.3)

由于 f(xk−1)f(x^{k-1}) 为常数，与优化无关，可以舍去。

步长选择

对于二次型目标函数，Hessian 矩阵为：

∇2f(x)=2ATA(4.4)\nabla^2 f(x) = 2A^T A \qquad (4.4)

由梯度的利普希茨连续性可知，∇2f(x)\nabla^2 f(x) 的最大特征值 LL 决定了梯度变化的上界。由于 ATAA^T A 对称半正定，其最大特征值等于谱半径（Spectral Radius）ρ(ATA)\rho(A^T A)，因此有：

L=2ρ(ATA)L = 2\rho(A^T A)

为保证迭代稳定性，取步长：

t=1L=12ρ(ATA)(4.5)t = \frac{1}{L} = \frac{1}{2\rho(A^T A)} \qquad (4.5)

其中 LL 为 2ATA2A^T A 的最大特征值，即 ATAA^T A 最大特征值的两倍。

公式化简

由二阶泰勒上界与利普希茨常数选取 tt，在舍去与 xx 无关的常数项后，可构造如下近似子问题：

xk=arg⁡min⁡x ∇f(xk−1)T(x−xk−1)+12t∥x−xk−1∥22(4.6)x^k = \arg\min_x\ \nabla f(x^{k-1})^T(x - x^{k-1}) + \frac{1}{2t}\|x - x^{k-1}\|_2^2 \qquad (4.6)

对上式“配方”可得：

∇f(xk−1)T(x−xk−1)+12t∥x−xk−1∥22=12t∥x−(xk−1−t∇f(xk−1))∥22+const,(4.7)\nabla f(x^{k-1})^T(x - x^{k-1}) + \frac{1}{2t}\|x - x^{k-1}\|_2^2 = \frac{1}{2t}\big\|x - \big(x^{k-1} - t\nabla f(x^{k-1})\big)\big\|_2^2 + \text{const}, \qquad (4.7)

因此，

xk=arg⁡min⁡x 12t∥x−rk∥22 ⇒ xk=rk=xk−1−t ∇f(xk−1).(4.8)x^k = \arg\min_x\ \frac{1}{2t}\|x - r^k\|_2^2 \ \Rightarrow\ x^k = r^k = x^{k-1} - t\,\nabla f(x^{k-1}). \qquad (4.8)

软阈值函数

令带正则的目标为

g(x)=f(x)+λ∥x∥1,λ>0.(4.9)g(x) = f(x) + \lambda\|x\|_1, \qquad \lambda>0. \qquad (4.9)

结合公式(4.8)可以得到如下子问题：

xk=arg⁡min⁡x 12t∥x−rk∥22+λ∥x∥1.(4.10)x^k = \arg\min_x\ \frac{1}{2t}\|x - r^k\|_2^2 + \lambda\|x\|_1. \qquad (4.10)

该问题按坐标可分解为对每个分量 xix_i 的一维优化：

xik=arg⁡min⁡u 12t(u−rik)2+λ∣u∣.(4.11)x_i^k = \arg\min_u\ \frac{1}{2t}(u - r_i^k)^2 + \lambda|u|. \qquad (4.11)

为了求解 xikx_i^k，我们先对该目标函数进行求导（注意 ∣u∣|u| 在 u=0u=0 不可导，需要用次梯度 ∂∣u∣\partial|u|）：

  • 对 u≠0u \neq 0 时，导数为

    ddu(12t(u−rik)2+λ∣u∣)=1t(u−rik)+λ sign⁡(u)\frac{d}{du} \Big( \frac{1}{2t}(u - r_i^k)^2 + \lambda|u| \Big) = \frac{1}{t}(u - r_i^k) + \lambda\,\operatorname{sign}(u)

  • 对 u=0u = 0 时，使用次梯度：

    ∂∣u∣=[−1,1],因此0∈1t(0−rik)+λ[−1,1].\partial |u| = [-1, 1], \quad \text{因此} \quad 0 \in \frac{1}{t}(0 - r_i^k) + \lambda [-1, 1].

令导数（或次梯度包含 0）等于 0，得到最优性条件：

0∈1t(u−rik)+λ ∂∣u∣.(4.12)0 \in \frac{1}{t}(u - r_i^k) + \lambda\,\partial|u|. \qquad (4.12)

根据 rikr_i^k 的大小，可以分三种情况求解：

  1. rik>λtr_i^k > \lambda t
    此时 u>0u>0，导数为 1t(u−rik)+λ=0\frac{1}{t}(u - r_i^k) + \lambda = 0，解得

    u=rik−λtu = r_i^k - \lambda t

  2. rik<−λtr_i^k < -\lambda t
    此时 u<0u<0，导数为 1t(u−rik)−λ=0\frac{1}{t}(u - r_i^k) - \lambda = 0，解得

    u=rik+λtu = r_i^k + \lambda t

  3. ∣rik∣≤λt|r_i^k| \le \lambda t
    此时 u=0u=0 可以满足次梯度条件，因为 0∈[−λ+−rikt,λ+−rikt]0 \in [-\lambda + \frac{-r_i^k}{t}, \lambda + \frac{-r_i^k}{t}]，因此

    u=0u = 0

综上，得到分段形式：

xik={rik−λt,rik>λt,0,∣rik∣≤λt,rik+λt,rik<−λt,(4.13)x_i^k= \begin{cases} r_i^k - \lambda t, & r_i^k > \lambda t,\\[2mm] 0, & |r_i^k|\le \lambda t,\\[2mm] r_i^k + \lambda t, & r_i^k < -\lambda t, \end{cases} \qquad (4.13)

这正是软阈值收缩算子 Sλt(⋅)S_{\lambda t}(\cdot)，因此

xk=Sλt(rk),Sθ(z)=sign⁡(z) max⁡(∣z∣−θ, 0).(4.14)x^k = S_{\lambda t}(r^k), \qquad S_\theta(z) = \operatorname{sign}(z)\,\max(|z|-\theta,\,0). \qquad (4.14)

将 rk=xk−1−t∇f(xk−1)=xk−1−2tAT(Axk−1−y)r^k = x^{k-1} - t\nabla f(x^{k-1}) = x^{k-1} - 2tA^T(Ax^{k-1}-y) 代入公式(4.14)即可。

ISTA 迭代格式，证毕。