带约束线性方程组的Lasso问题求解——软阈值迭代法（ISTA）算法详解
字数 3595 2025-12-24 12:43:29

带约束线性方程组的Lasso问题求解——软阈值迭代法（ISTA）算法详解

一、问题描述

我们考虑一个在机器学习、信号处理和统计学中经常遇到的线性代数问题：带L1正则化（Lasso）的线性回归问题，也称为“基追踪去噪”问题。它通常用于特征选择或获得稀疏解。

给定一个数据矩阵 \(A \in \mathbb{R}^{m \times n}\)（通常 \(m \neq n\)，\(m\) 为样本数，\(n\) 为特征数），观测向量 \(b \in \mathbb{R}^{m}\)，以及一个正则化参数 \(\lambda > 0\)，我们的目标是求解以下凸优化问题：

\[\min_{x \in \mathbb{R}^{n}} \left\{ F(x) = \frac{1}{2} \| Ax - b \|_2^2 + \lambda \| x \|_1 \right\} \]

其中：

  • \(\frac{1}{2} \| Ax - b \|_2^2\) 是数据拟合项（最小二乘误差）。
  • \(\| x \|_1 = \sum_{i=1}^{n} |x_i|\) 是L1范数，用来促使解向量 \(x\) 变得稀疏（即许多分量为零）。
  • \(\lambda\) 用来平衡拟合精度和解的稀疏性。

这个问题没有像普通最小二乘那样的闭式解，因为目标函数不可微（在 \(x_i = 0\) 处不可导）。软阈值迭代法（Iterative Shrinkage-Thresholding Algorithm, ISTA）是解决此类非光滑优化问题的一种经典且有效的一阶方法。

二、算法推导与核心思想

2.1 问题转化与近似

我们的目标函数 \(F(x)\) 由两部分组成：

  • 光滑部分 \(f(x) = \frac{1}{2} \| Ax - b \|_2^2\)。
  • 非光滑部分 \(g(x) = \lambda \| x \|_1\)。

ISTA属于“近端梯度法”（Proximal Gradient Method）的范畴。核心思想是：在每次迭代中，先沿着光滑部分的负梯度方向走一步，然后对这个结果施加一个非光滑项对应的“近端算子”（对于L1范数来说，就是软阈值函数）。

2.2 梯度与Lipschitz常数

首先计算光滑部分 \(f(x)\) 的梯度：

\[\nabla f(x) = A^T (Ax - b) \]

梯度函数 \(\nabla f(x)\) 是Lipschitz连续的。其Lipschitz常数 \(L\) 满足：

\[\| \nabla f(x) - \nabla f(y) \|_2 \le L \| x - y \|_2, \quad \forall x, y \]

可以证明，\(L\) 的上界是 \(\| A^T A \|_2\)，即矩阵 \(A^T A\) 的最大特征值（谱范数）。实际操作中，通常通过幂法或直接计算 \(A^T A\) 的特征值来估计 \(L\)，或者设置一个可调节的步长参数 \(\eta = 1/L\)。

2.3 软阈值算子

软阈值算子 \(S_{\tau}\) 是解决L1正则化问题的关键。它对向量的每个分量独立操作：

\[[S_{\tau}(z)]_i = \text{sign}(z_i) \cdot \max \{ |z_i| - \tau, 0 \} \]

其作用是：若 \(|z_i|\) 小于阈值 \(\tau\)，则将该分量置为0；否则，将其向零的方向“压缩” \(\tau\) 个单位，符号保持不变。这正是L1范数诱导稀疏性的直接体现。

2.4 ISTA迭代公式

ISTA的标准迭代格式（固定步长版本）如下：

\[x^{(k+1)} = S_{\eta \lambda} \left( x^{(k)} - \eta \nabla f(x^{(k)}) \right) \]

代入梯度的具体表达式，得到：

\[x^{(k+1)} = S_{\eta \lambda} \left( x^{(k)} - \eta A^T (A x^{(k)} - b) \right) \]

其中：

  • \(\eta\) 是步长，需满足 \(0 < \eta \le 1/L\) 以保证收敛。
  • \(S_{\eta \lambda}\) 中的阈值是 \(\eta \lambda\)，这是步长与正则化参数的乘积。

三、算法步骤详解

假设我们已给定矩阵 \(A\)，向量 \(b\)，参数 \(\lambda\)，并设置好步长 \(\eta\) 或通过计算得到 \(L\)。

步骤 1：初始化

  • 设置初始解 \(x^{(0)}\)，通常可以设为零向量 \(0 \in \mathbb{R}^{n}\)。
  • 设置最大迭代次数 \(K_{max}\) 和容忍误差 \(\epsilon\)（例如 \(10^{-6}\)），用于判断收敛。
  • 计算或估计梯度Lipschitz常数 \(L\)（例如 \(L = \|A^T A\|_2\)），并令步长 \(\eta = 1 / L\)（或更保守的 \(\eta = 0.9 / L\)）。
  • 令迭代计数器 \(k = 0\)。

步骤 2：主迭代循环

当 \(k < K_{max}\) 且未达到收敛标准时，重复以下步骤：

  1. 计算梯度：
    计算当前残差和梯度：

\[ r^{(k)} = A x^{(k)} - b \]

\[ g^{(k)} = A^T r^{(k)} \]

这里的 \(g^{(k)}\) 就是 \(\nabla f(x^{(k)})\)。

  1. 梯度下降步：
    沿着负梯度方向走一步：

\[ v^{(k)} = x^{(k)} - \eta g^{(k)} \]

这个 \(v^{(k)}\) 是只考虑光滑部分的梯度下降结果。

  1. 软阈值操作（近端映射）：
    对 \(v^{(k)}\) 的每一个分量 \(v_i^{(k)}\) 应用软阈值算子：

\[ x_i^{(k+1)} = \text{sign}(v_i^{(k)}) \cdot \max \{ |v_i^{(k)}| - \eta \lambda, 0 \}, \quad i=1, \dots, n \]

这个步骤可以并行完成。它直接产生了新的迭代点 \(x^{(k+1)}\)。

  1. 收敛性检查：
    检查相邻两次迭代解的变化是否足够小：

\[ \| x^{(k+1)} - x^{(k)} \|_2 / \| x^{(k)} \|_2 < \epsilon \]

或者检查目标函数值 \(F(x^{(k+1)})\) 的变化是否小于阈值。如果满足，则提前终止循环。

  1. 更新迭代计数器：

\[ k \leftarrow k + 1 \]

步骤 3：输出结果

  • 最后一次迭代得到的 \(x^{(k)}\) 即为所求的稀疏近似解。

四、几何解释与收敛性

  • 几何解释： 每一步先做梯度下降（试图最小化拟合误差），这会使解远离稀疏性；紧接着的软阈值操作（施加L1惩罚）又将解“拉向”坐标轴（使得某些分量变为零）。最终，算法在这两种力量之间达到一个平衡点。
  • 收敛性： ISTA的收敛速度是次线性的，即目标函数值的误差以 \(O(1/k)\) 的速度下降。虽然不快，但每一步的计算成本很低（主要是两次矩阵-向量乘法和一个 \(O(n)\) 的软阈值操作），使其对于大规模问题非常实用。

五、加速版本：FISTA（Fast ISTA）

为了提高收敛速度，Nesterov提出了加速梯度方法。将其与ISTA结合，就是著名的 FISTA（Fast Iterative Shrinkage-Thresholding Algorithm）。
与ISTA的主要区别在于，它引入了动量项。迭代步骤如下：

  1. 引入辅助序列 \(\{y^{(k)}\}\)。
  2. 从 \(y^{(k)}\) 出发进行梯度计算和软阈值操作，得到 \(x^{(k+1)}\)。
  3. 按特定公式更新 \(y^{(k+1)}\)，使其包含之前迭代的“动量”。

FISTA可以达到 \(O(1/k^2)\) 的收敛速度，在实际应用中比ISTA快得多。

通过以上步骤，我们详细剖析了用于求解Lasso问题的ISTA算法。它巧妙地将梯度下降与软阈值操作结合，成为处理带L1正则化线性模型的基础工具。

带约束线性方程组的Lasso问题求解——软阈值迭代法（ISTA）算法详解 一、问题描述 我们考虑一个在机器学习、信号处理和统计学中经常遇到的线性代数问题： 带L1正则化（Lasso）的线性回归问题 ，也称为“基追踪去噪”问题。它通常用于特征选择或获得稀疏解。 给定一个数据矩阵 \( A \in \mathbb{R}^{m \times n} \)（通常 \( m \neq n \)，\( m \) 为样本数，\( n \) 为特征数），观测向量 \( b \in \mathbb{R}^{m} \)，以及一个正则化参数 \( \lambda > 0 \)，我们的目标是求解以下凸优化问题： \[ \min_ {x \in \mathbb{R}^{n}} \left\{ F(x) = \frac{1}{2} \| Ax - b \|_ 2^2 + \lambda \| x \|_ 1 \right\} \] 其中： \( \frac{1}{2} \| Ax - b \|_ 2^2 \) 是数据拟合项（最小二乘误差）。 \( \| x \| 1 = \sum {i=1}^{n} |x_ i| \) 是L1范数，用来促使解向量 \( x \) 变得稀疏（即许多分量为零）。 \( \lambda \) 用来平衡拟合精度和解的稀疏性。 这个问题没有像普通最小二乘那样的闭式解，因为目标函数不可微（在 \( x_ i = 0 \) 处不可导）。软阈值迭代法（Iterative Shrinkage-Thresholding Algorithm, ISTA）是解决此类非光滑优化问题的一种经典且有效的一阶方法。 二、算法推导与核心思想 2.1 问题转化与近似 我们的目标函数 \( F(x) \) 由两部分组成： 光滑部分 \( f(x) = \frac{1}{2} \| Ax - b \|_ 2^2 \)。 非光滑部分 \( g(x) = \lambda \| x \|_ 1 \)。 ISTA属于“近端梯度法”（Proximal Gradient Method）的范畴。核心思想是： 在每次迭代中，先沿着光滑部分的负梯度方向走一步，然后对这个结果施加一个非光滑项对应的“近端算子”（对于L1范数来说，就是软阈值函数） 。 2.2 梯度与Lipschitz常数 首先计算光滑部分 \( f(x) \) 的梯度： \[ \nabla f(x) = A^T (Ax - b) \] 梯度函数 \( \nabla f(x) \) 是Lipschitz连续的。其Lipschitz常数 \( L \) 满足： \[ \| \nabla f(x) - \nabla f(y) \|_ 2 \le L \| x - y \|_ 2, \quad \forall x, y \] 可以证明，\( L \) 的上界是 \( \| A^T A \|_ 2 \)，即矩阵 \( A^T A \) 的最大特征值（谱范数）。实际操作中，通常通过幂法或直接计算 \( A^T A \) 的特征值来估计 \( L \)，或者设置一个可调节的步长参数 \( \eta = 1/L \)。 2.3 软阈值算子 软阈值算子 \( S_ {\tau} \) 是解决L1正则化问题的关键。它对向量的每个分量独立操作： \[ [ S_ {\tau}(z)]_ i = \text{sign}(z_ i) \cdot \max \{ |z_ i| - \tau, 0 \} \] 其作用是：若 \( |z_ i| \) 小于阈值 \( \tau \)，则将该分量置为0；否则，将其向零的方向“压缩” \( \tau \) 个单位，符号保持不变。这正是L1范数诱导稀疏性的直接体现。 2.4 ISTA迭代公式 ISTA的标准迭代格式（固定步长版本）如下： \[ x^{(k+1)} = S_ {\eta \lambda} \left( x^{(k)} - \eta \nabla f(x^{(k)}) \right) \] 代入梯度的具体表达式，得到： \[ x^{(k+1)} = S_ {\eta \lambda} \left( x^{(k)} - \eta A^T (A x^{(k)} - b) \right) \] 其中： \( \eta \) 是步长，需满足 \( 0 < \eta \le 1/L \) 以保证收敛。 \( S_ {\eta \lambda} \) 中的阈值是 \( \eta \lambda \)，这是步长与正则化参数的乘积。 三、算法步骤详解 假设我们已给定矩阵 \( A \)，向量 \( b \)，参数 \( \lambda \)，并设置好步长 \( \eta \) 或通过计算得到 \( L \)。 步骤 1：初始化 设置初始解 \( x^{(0)} \)，通常可以设为零向量 \( 0 \in \mathbb{R}^{n} \)。 设置最大迭代次数 \( K_ {max} \) 和容忍误差 \( \epsilon \)（例如 \( 10^{-6} \)），用于判断收敛。 计算或估计梯度Lipschitz常数 \( L \)（例如 \( L = \|A^T A\|_ 2 \)），并令步长 \( \eta = 1 / L \)（或更保守的 \( \eta = 0.9 / L \)）。 令迭代计数器 \( k = 0 \)。 步骤 2：主迭代循环 当 \( k < K_ {max} \) 且未达到收敛标准时，重复以下步骤： 计算梯度 ： 计算当前残差和梯度： \[ r^{(k)} = A x^{(k)} - b \] \[ g^{(k)} = A^T r^{(k)} \] 这里的 \( g^{(k)} \) 就是 \( \nabla f(x^{(k)}) \)。 梯度下降步 ： 沿着负梯度方向走一步： \[ v^{(k)} = x^{(k)} - \eta g^{(k)} \] 这个 \( v^{(k)} \) 是只考虑光滑部分的梯度下降结果。 软阈值操作（近端映射） ： 对 \( v^{(k)} \) 的每一个分量 \( v_ i^{(k)} \) 应用软阈值算子： \[ x_ i^{(k+1)} = \text{sign}(v_ i^{(k)}) \cdot \max \{ |v_ i^{(k)}| - \eta \lambda, 0 \}, \quad i=1, \dots, n \] 这个步骤可以并行完成。它直接产生了新的迭代点 \( x^{(k+1)} \)。 收敛性检查 ： 检查相邻两次迭代解的变化是否足够小： \[ \| x^{(k+1)} - x^{(k)} \|_ 2 / \| x^{(k)} \|_ 2 < \epsilon \] 或者检查目标函数值 \( F(x^{(k+1)}) \) 的变化是否小于阈值。如果满足，则提前终止循环。 更新迭代计数器 ： \[ k \leftarrow k + 1 \] 步骤 3：输出结果 最后一次迭代得到的 \( x^{(k)} \) 即为所求的稀疏近似解。 四、几何解释与收敛性 几何解释 ： 每一步先做梯度下降（试图最小化拟合误差），这会使解远离稀疏性；紧接着的软阈值操作（施加L1惩罚）又将解“拉向”坐标轴（使得某些分量变为零）。最终，算法在这两种力量之间达到一个平衡点。 收敛性 ： ISTA的收敛速度是次线性的，即目标函数值的误差以 \( O(1/k) \) 的速度下降。虽然不快，但每一步的计算成本很低（主要是两次矩阵-向量乘法和一个 \( O(n) \) 的软阈值操作），使其对于大规模问题非常实用。 五、加速版本：FISTA（Fast ISTA） 为了提高收敛速度，Nesterov提出了加速梯度方法。将其与ISTA结合，就是著名的 FISTA（Fast Iterative Shrinkage-Thresholding Algorithm） 。 与ISTA的主要区别在于，它引入了动量项。迭代步骤如下： 引入辅助序列 \( \{y^{(k)}\} \)。 从 \( y^{(k)} \) 出发进行梯度计算和软阈值操作，得到 \( x^{(k+1)} \)。 按特定公式更新 \( y^{(k+1)} \)，使其包含之前迭代的“动量”。 FISTA可以达到 \( O(1/k^2) \) 的收敛速度，在实际应用中比ISTA快得多。 通过以上步骤，我们详细剖析了用于求解Lasso问题的ISTA算法。它巧妙地将梯度下降与软阈值操作结合，成为处理带L1正则化线性模型的基础工具。