1. 近端梯度法解 LASSO，先搞清楚问题本身

LASSO（Least Absolute Shrinkage and Selection Operator）这个问题，搞稀疏优化、压缩感知、特征选择的人基本都会碰到。形式上讲，它在做的是求解这样一个优化问题：

$$ \min_{x} ; \frac{1}{2} | Ax - b |_2^2 + \lambda | x |_1 $$

这里 $A$ 通常是 $m \times n$ 的测量矩阵，$b$ 是我们的观测向量，$x$ 是待恢复的稀疏系数。近端梯度法要解决的就是这个问题，而且解法很干净：把目标函数拆成光滑的二次项加上非光滑的 $L_1$ 正则项，然后用近端算子（proximal operator）处理非光滑部分。

网上讲这个的数学推导不少，但很多资料要么只给结论，要么代码跑起来一堆问题。我结合自己实现的经验，把原理、推导、代码、坑位一次性捋清楚。

1.1 为什么不能用普通梯度下降

很多人一开始看到 LASSO 会想：这不就是一个带正则的线性回归吗？直接梯度下降不就行了？

问题出在 $|x|_1$ 这一项。$L_1$ 范数在零点是不可导的，具体来说，它在 $x_i = 0$ 处的次微分是一个区间 $[-λ, λ]$，而不是一个确定的梯度值。普通梯度下降需要计算梯度，那零点处到底取多少？没有唯一的答案。

光这一点还不够。更现实的问题是：LASSO 的解通常希望在零点处精确等于 0，而不是一个很小很小的数。如果你用带 $L_2$ 正则的梯度下降，比如岭回归的那种，得到的解基本都是非零的密集向量，稀疏性完全体现不出来。加上普通梯度下降对不可导点只能做次梯度下降，而次梯度法收敛速度是 $O(1/\sqrt{k})$，慢得让人发慌。实测下来，同一个数据集，次梯度法可能要跑几万步才能有个像样的结果，近端梯度法几百步就收敛了。

近端梯度法的高明之处在于把"光滑部分求梯度 + 非光滑部分做近端映射"拆开，两步交替迭代，理论收敛速度是 $O(1/k)$，加了 NesterOV 加速后能到 $O(1/k^2)$。这个速度差异在实际项目中是决定性的。

1.2 LASSO 的本质是稀疏恢复

LASSO 为什么值得专门写？因为它干的事很有实用价值：

• 高维数据特征选择：$n$ 很大（特征很多），$m$ 很小（样本很少）的时候，模型希望只用少数几个特征做预测。
• 压缩感知信号恢复：从欠定测量中恢复原始稀疏信号，这个在图像重建、医学成像里是核心方法。
• 噪声鲁棒性：$L_1$ 正则天然比 $L_2$ 正则对离群点更鲁棒，因为它对误差的惩罚是线性的，不像平方误差那样把异常值放大。

而且近端梯度法这个名字听起来唬人，实际算起来核心就一个公式。下面这两节我把数学原理和代码对应着讲清楚，你只要理解了"近端算子"这个概念，写代码就是水到渠成的事。

2. 近端梯度法的核心原理拆解

2.1 从"先梯度下降，再近端映射"说起

近端梯度法解决的是这样一类复合优化问题：

$$ \min_x ; f(x) + g(x) $$

其中 $f$ 是光滑凸函数（比如 $\frac{1}{2}|Ax - b|_2^2$），$g$ 是凸但可能非光滑的函数（比如 $\lambda |x|_1$）。

核心想法是：在迭代的每一步，我们不直接对整个目标函数求梯度，而是先只对光滑部分 $f$ 做一步梯度下降，得到临时点 $z$：

$$ z = x^{(k)} - t_k \nabla f(x^{(k)}) $$

然后把这个临时点 $z$ 投影到"考虑 $g$ 的形状"后应该去的地方。这个投影操作近端算子定义：

$$ \operatorname{prox}_{t g}(z) = \arg\min_x ; \frac{1}{2} | x - z |_2^2 + t , g(x) $$

也就是说，找一个点 $x$，它既要离临时点 $z$ 近，又要让 $g(x)$ 尽量小。两步加起来就是：

$$ x^{(k+1)} = \operatorname{prox}_{t_k g}\left( x^{(k)} - t_k \nabla f(x^{(k)}) \right) $$

用生活化的类比来说：想象你在一个山坡上（光滑的 $f$），想往低处走，但脚下有块玻璃挡板（非光滑的 $g$），你不能直接穿过去，只能先沿着坡面走一步，再被挡板弹回允许走的位置。近端算子就是那个"弹回"的动作。

2.2 软阈值算子：LASSO 的专属近端算子

对 LASSO 来说，$g(x) = \lambda |x|_1$，$t g(x) = t \lambda \sum_i |x_i|$。由于 $L_1$ 范数是分量可分离的，这个近端算子可以逐个坐标独立计算，每个坐标的结果就是著名的软阈值（soft-thresholding）公式：

$$ \operatorname{prox}_{t\lambda |\cdot|_1}(z)_i = \operatorname{sign}(z_i) \cdot \max\left( |z_i| - t\lambda, 0 \right) $$

换成更直观的写法：

def soft_threshold(z, threshold): return np.sign(z) * np.maximum(np.abs(z) - threshold, 0)

它做的事是：如果 $z_i$ 大于阈值 $\theta = t\lambda$，就把它往零方向压缩 $\theta$；如果小于阈值，直接置零。这个"压一压、不够就清零"的操作，就是 LASSO 稀疏性的来源。

我把这个操作叫做"绝对值上的平移门"：大于门宽的留下，小于门宽的归零。这也是为什么近端梯度法解 LASSO 看起来格外优雅——近端算子把原本很棘手的非光滑问题，变成了一个逐元素操作的闭式解。

整个迭代过程总结如下：

1. 初始化 $x^{(0)}$，比如零向量。
2. 对第 $k$ 步：
   • 计算梯度 $\nabla f(x^{(k)}) = A^T(Ax^{(k)} - b)$。
   • 做梯度下降临时更新：$z = x^{(k)} - t_k \nabla f(x^{(k)})$。
   • 应用软阈值算子：$x^{(k+1)} = S_{t_k\lambda}(z)$。
3. 重复直到收敛。

3. Python 实现代码与实验对比

3.1 基础版实现（ISTA）

迭代软阈值算法（ISTA, Iterative Shrinkage-Thresholding Algorithm）就是上面三步的直白翻译。我用 numpy 写了一个完整实现：

import numpy as np def soft_threshold(z, threshold): """软阈值算子：向量逐元素操作""" return np.sign(z) * np.maximum(np.abs(z) - threshold, 0) def ista(A, b, lam, max_iter=1000, tol=1e-8, x_init=None): """ ISTA 解 LASSO： min 0.5 ||Ax - b||^2 + lam * ||x||_1 参数 ---- A : m x n 测量矩阵 b : m 维观测向量 lam : 正则化系数 """ m, n = A.shape if x_init is None: x = np.zeros(n) else: x = x_init.copy() # 估算 Lipschitz 常数：最大的特征值可以用幂迭代法或直接取最大奇异值平方近似 # 这里为了稳定，直接用最大奇异值平方 L = np.linalg.svd(A, compute_uv=False)[0] ** 2 t = 1.0 / L # 固定步长 history = [] for k in range(max_iter): grad = A.T @ (A @ x - b) z = x - t * grad x_new = soft_threshold(z, t * lam) # 收敛判断：用相对变化量 change = np.max(np.abs(x_new - x)) history.append(change) x = x_new if change < tol: break return x, history

有个细节很多人会踩坑：$L$ 的估算。理论上步长 $t$ 要满足 $t \le 1/L$，其中 $L$ 是 $\nabla f$ 的 Lipschitz 常数，也就是 $A^TA$ 的最大特征值。直接用奇异值平方是标准做法，但每次调用np.linalg.svd在大矩阵上是额外开销。实际工程中，如果矩阵结构特殊（比如测量矩阵是部分傅里叶矩阵、随机投影矩阵），可以直接用一个已知上界。

如果不确定 $L$ 准不准，可以加一步回溯线搜索：每次迭代先选一个较大的步长尝试，如果 $f(x_{k+1}) + g(x_{k+1})$ 没有比线性化模型减少足够多，就把步长缩小（乘以一个小于 1 的系数）。这就是是"回溯近端梯度"做法，代码上量也不大，后面会给出示例。

3.2 加速版实现（FISTA）

基础版 ISTA 的收敛速度是 $O(1/k)$。什么意思？误差大概以 1 除以迭代次数的速度在降。想要误差到 $10^{-4}$，可能要一万步。加速版 FISTA（Fast Iterative Shrinkage-Thresholding Algorithm）引入了一个外推动量项，把收敛速度提升到 $O(1/k^2)$，迭代六百步效果就差不多能赶上。

FISTA 的迭代公式和 ISTA 的主要区别是：每一步不是用当前点 $x_k$ 做梯度更新，而是用一个"外推点" $y_k$：

$$ y_k = x_k + \frac{\theta_{k-1} - 1}{\theta_k}(x_k - x_{k-1}) $$

$$ x_{k+1} = \operatorname{prox}_{t g}\left( y_k - t \nabla f(y_k) \right) $$

其中 $\theta_k$ 递推：$\theta_{k+1} = \frac{1 + \sqrt{1 + 4\theta_k^2}}{2}$，初值 $\theta_0 = 1$。

def fista(A, b, lam, max_iter=1000, tol=1e-8, x_init=None): """加速近端梯度法解 LASSO""" m, n = A.shape if x_init is None: x = np.zeros(n) else: x = x_init.copy() L = np.linalg.svd(A, compute_uv=False)[0] ** 2 t = 1.0 / L y = x.copy() theta = 1.0 history = [] for k in range(max_iter): grad = A.T @ (A @ y - b) z = y - t * grad x_new = soft_threshold(z, t * lam) theta_new = (1 + np.sqrt(1 + 4 * theta**2)) / 2 # 外推系数 beta = (theta - 1) / theta_new y = x_new + beta * (x_new - x) change = np.max(np.abs(x_new - x)) history.append(change) x = x_new theta = theta_new if change < tol: break return x, history

注意一个坑：FISTA 中 y 在外推时如果步长 $t$ 选得偏大（即 $L$ 估计偏小），动量外推会让迭代发散。这也是很多从文献里抄代码的人跑崩的根本原因。稳健的做法是要么保守估计 $L$，要么加上重启策略（后面细说）。

3.3 数值实验：ISTA vs FISTA

我构造了一个随机稀疏恢复问题来对比两种算法的表现：真解 $x_{\text{true}}$ 是一个含 30 个非零元、幅值服从标准正态的 500 维向量；测量矩阵 $A$ 是 $200 \times 500$ 的随机高斯矩阵；观测 $b = A x_{\text{true}} + \epsilon$，噪声方差设为 0.01。正则系数 $\lambda$ 取 $0.1 |A^T b|_\infty$，这是实践中一个常用的启发式取值。

np.random.seed(42) m, n = 200, 500 k_true = 30 A = np.random.randn(m, n) / np.sqrt(m) x_true = np.zeros(n) idx = np.random.choice(n, k_true, replace=False) x_true[idx] = np.random.randn(k_true) b = A @ x_true + 0.01 * np.random.randn(m) lam = 0.1 * np.max(np.abs(A.T @ b)) x_ista, hist_ista = ista(A, b, lam, max_iter=5000) x_fista, hist_fista = fista(A, b, lam, max_iter=5000)

评估指标用两个：目标函数值和相对误差 $|x^{(k)} - x^| / |x^|$。下图的核心结论是：

• FISTA 在 200 步左右就已经收敛到稳定解，ISTA 需要 2000 步以上。
• 两者的最终解几乎一致（相对误差小于 $10^{-6}$），说明收敛到了同一个最优点。
• 稀疏模式恢复准确：非零元位置与真解完全一致，幅值偏差在噪声水平内。

更直观的对比数据（迭代次数到某个误差阈值所需量级）：

这个差距就是加速项的威力。实际业务中如果矩阵维度过千上万，ISTA 可能跑几分钟都收不敛，FISTA 几十秒就能出结果，所以能用 FISTA 就别老老实实跑 ISTA 了。

4. 参数调优、停止准则与工程坑位

4.1 正则系数 $\lambda$ 怎么选

$\lambda$ 是 LASSO 模型里最敏感的超参数，它直接决定解的稀疏程度。$\lambda$ 太大，解会全部变成零；$\lambda$ 太小，稀疏性体现不出来。实际中最常用的选择手段是交叉验证：

• 把数据分成训练集和验证集。
• 在一组候选 $\lambda$ 上分别解 LASSO，用验证集误差挑最优的。
• 常见候选序列是从 $\lambda_{\max}$ 到 $\lambda_{\max} \times 10^{-4}$ 的对数等比数列。

$\lambda_{\max}$ 的确定有一个很有用的闭合表达式：当 $\lambda \ge |A^T b|\infty$ 时，LASSO 的最优解就是全零向量。证明不复杂，KKT 条件的对偶可行性直接给出这个界。所以实践中 $\lambda{\max} = |A^T b|_\infty$，然后向下扫描。

从经验上看，一个比较稳妥的操作是：先用 $\lambda_{\max}$ 算出"全都稀疏"的起点，然后逐步缩小 $\lambda$，观察非零元素个数和交叉验证误差。注意不要把 $\lambda$ 压得太低，否则模型过拟合噪声，稀疏解的优势就没了。

4.2 步长与 $L$ 的估算细节

固定步长 $t = 1/L$ 是最省事的做法，但前提是 $L$ 估得准。$L$ 本质上就是 $A^TA$ 的最大特征值。在 numpy 中：

# 方法一：直接奇异值计算（稳定但小贵） L = np.linalg.svd(A, compute_uv=False)[0] ** 2 # 方法二：幂迭代法，只算最大特征值，大矩阵下更快 def power_iteration(A, num_iter=100): m, n = A.shape v = np.random.randn(n) v /= np.linalg.norm(v) for _ in range(num_iter): v = A.T @ (A @ v) v /= np.linalg.norm(v) return np.dot(v, A.T @ (A @ v)) L_approx = power_iteration(A)

注意：如果 $A$ 是列归一化后的随机矩阵，$L$ 通常会小于某个常数；但如果不做任何归一化，$L$ 可能很大，步长会非常小，收敛就慢。因此建议在做实验前把 $A$ 的每一列做归一化处理：

col_norms = np.linalg.norm(A, axis=0) A_norm = A / col_norms

此时目标函数会有点变化，等价于对 $x$ 做了缩放，但稀疏模式不变。列归一化后 $L$ 的上界更好估计，数值也更稳定。对实测中测量矩阵列尺度差异极大的情况，这步基本是必须的。

如果不想要固定步长，可以上回溯近端梯度（backtracking proximal gradient）：

def proximal_gradient_backtracking(A, b, lam, max_iter=1000, tol=1e-8): m, n = A.shape x = np.zeros(n) L = np.linalg.svd(A, compute_uv=False)[0] ** 2 # 初始步长给大一点，后面不满足条件再缩 t = 2.0 / L eta = 0.8 # 缩小系数 for _ in range(max_iter): grad = A.T @ (A @ x - b) while True: z = soft_threshold(x - t * grad, t * lam) # Armijo 型条件：目标函数下降足够多 f_x = 0.5 * np.linalg.norm(A @ x - b)**2 + lam * np.abs(x).sum() f_z = 0.5 * np.linalg.norm(A @ z - b)**2 + lam * np.abs(z).sum() model = f_x + grad @ (z - x) + 0.5 / t * np.linalg.norm(z - x)**2 if f_z <= model: break t *= eta x = z return x

4.3 停止准则：别只看目标函数

很多初学者判断是否收敛只看目标函数值不再下降。这样做很容易被"平坦区"骗到：目标函数在最优解附近确实很平，差一点数值差很小，但解本身可能离最优解还很远。

我建议用两套指标结合：

1. 相对变化量：$|x^{(k+1)} - x^{(k)}|_\infty$ 小于某个阈值（比如 $10^{-8}$）。
2. KKT 残差：检查近端梯度迭代的"不动点条件"是否满足，即 $|x - \operatorname{prox}_g(x - t\nabla f(x))|$ 是否足够小。这个比目标函数下降更可靠。

自己实现的调试阶段，建议每 100 步打印一次目标函数值、相对变化、非零元素个数。有一次我调一个稀疏字典学习任务，目标函数看着已经纹丝不动了，但非零元个数还在缓慢变化，最后发现是停止准则过于宽松导致的，换成 KKT 残差后问题就暴露了。

4.4 工程实现上的常见坑位

坑一：软阈值公式写错符号。软阈值是 $\operatorname{sign}(z)\max(|z| - \theta, 0)$，不是 $\max(\operatorname{sign}(z)z - \theta, 0)$。后者会把负数方向搞乱。对称正负都要压，别只压正数。

坑二：忘记列归一化导致步长太小。测量矩阵各列尺度差异极大时，$L$ 会被最大奇异值主导，导致全局步长极小，收敛慢得离谱。归一化之后效果立竿见影。

坑三：FISTA 发散了还想硬调。FISTA 对步长比 ISTA 敏感得多。当出现目标函数"先降后升"的振荡时，最快排查办法就是检查 L 的估计是否是下界。最稳的方案是给 FISTA 加一个简单的重启策略：监测到目标函数上升或者外推方向与下降方向不匹配时，把 $\theta$ 重置为 1。这个技巧在理论上不损失加速效果，实际中特别适合工程落地。

坑四：稀疏性的"假恢复"。如果 $\lambda$ 选得太小，恢复出来的非零元个数比真实情况多很多，信噪比高的数据集上看起来误差很小，但换到噪声大的数据集就崩了。不要只看训练误差，要看验证集上的表现。

5. 常见问题速查表与排错思路

我把做近端梯度法过程中最常见的问题整理成了一张速查表，方便遇到问题快速定位：

5.1 KKT 残差怎么算

近端梯度迭代的不动点条件是：

$$ x^* = \operatorname{prox}_{t g}\left( x^* - t \nabla f(x^*) \right) $$

所以 KKT 残差可以定义为：

def kkt_residual(x, A, b, lam, t): grad_f = A.T @ (A @ x - b) z = x - t * grad_f return np.max(np.abs(x - soft_threshold(z, t * lam)))

这个量越小，说明 $x$ 越接近真正的最优解。实际使用中，我把tol=1e-8定为默认，对大多数中等规模问题来说已经非常精确了。如果追求工程上"足够好"的精度，1e-7到1e-6通常就够用了，没必要死磕到 1e-12，否则最后几百次迭代都在修正小数点后第 9 位，纯属浪费时间。

5.2 特征正交化与稀疏恢复的一致性

有个隐藏问题是：当 $A$ 的列高度相关时，LASSO 的恢复结果可能不稳定，也就是说 $A$ 中相关性强的特征会随机分配权重。这个不是求解算法的问题，而是模型本身的多重共线性导致非零元选择对噪声敏感。如果你在做特征选择，建议在 LASSO 解出来后做一个最小二乘精修（debias）：对选出的非零列集 $\mathcal{S}$，用 $(A_\mathcal{S}^T A_\mathcal{S})^{-1} A_\mathcal{S}^T b$ 重新估计幅值。这是 LASSO 后处理的标准操作，能显著降低幅值偏差。

因为 LASSO 的 $L_1$ 惩罚天生对幅值有压缩，这个压缩量通常在 10%~20% 的水平。精修后才是真正干净的幅值估计。

5.3 大规模场景的曲线救国方案

近端梯度法在高维大规模问题上还有几个变体：

• 随机近端梯度（Stochastic Proximal Gradient）：每次迭代只用一个小批量样本计算梯度，适合 $m$ 极大的场景，应用在流式数据或大规模样本下的在线稀疏学习。
• 近端坐标下降（Proximal Coordinate Descent）：每次只更新一个坐标，对特征数 $n$ 极大但期望稀疏解的场景往往比近端梯度更快，在电子表格级数据上表现优秀，因为它天然利用了稀疏性。
• ADMM：如果要加等式约束或者复杂复合正则（比如 $L_1+L_2$、全变分），ADMM 是把问题拆成更小子问题并行求解的利器。

但如果你只是想快速解一个常规大小的 LASSO，FISTA 是最具性价比的选择：代码量不到 30 行，理论性质好，实践稳。这也是它在优化工具包里出现频率这么高的原因。

我在实际调参中还有个体会：不要迷信固定步长，除非矩阵设计得非常规整。真实数据里的矩阵条件数往往不理想，最省心的是配合一个简单的预警机制——每 20 步检查目标函数是否比前 50 步的最小值高出 1%，如果高出了就强行缩小步长或重启动量。这个小监控可以在不牺牲效率的前提下，避免项目进入半夜跑完发现崩了的窘境。

最后再分享一个小技巧：近端梯度类算法本质上干的事情就是"展开一个稀疏性约束的递归结构"。如果你从深度学习视角去看 FISTA 的迭代式，会发现它和一个带残差连接的神经网络非常相似——每一层的结构就是"线性变换 + 阈值非线性"。这个视角在做 learned ISTA、即把迭代展开成网络结构的时候非常有用。我有一次做图像稀疏重建，就是初始化网络参数为 FISTA 的 $A^T$ 和步长，网络训练速度比随机初始化快了一个数量级。把这层联系记在心里，近端梯度法对你的价值就远不止解一个 LASSO 那么简单了。