1. 近端梯度法解 LASSO先搞清楚问题本身LASSOLeast 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^{(k1)} \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 看起来格外优雅——近端算子把原本很棘手的非光滑问题变成了一个逐元素操作的闭式解。整个迭代过程总结如下初始化 $x^{(0)}$比如零向量。对第 $k$ 步计算梯度 $\nabla f(x^{(k)}) A^T(Ax^{(k)} - b)$。做梯度下降临时更新$z x^{(k)} - t_k \nabla f(x^{(k)})$。应用软阈值算子$x^{(k1)} S_{t_k\lambda}(z)$。重复直到收敛。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_iter1000, tol1e-8, x_initNone): 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_uvFalse)[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_{k1}) g(x_{k1})$ 没有比线性化模型减少足够多就把步长缩小乘以一个小于 1 的系数。这就是是回溯近端梯度做法代码上量也不大后面会给出示例。3.2 加速版实现FISTA基础版 ISTA 的收敛速度是 $O(1/k)$。什么意思误差大概以 1 除以迭代次数的速度在降。想要误差到 $10^{-4}$可能要一万步。加速版 FISTAFast 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_{k1} \operatorname{prox}_{t g}\left( y_k - t \nabla f(y_k) \right) $$其中 $\theta_k$ 递推$\theta_{k1} \frac{1 \sqrt{1 4\theta_k^2}}{2}$初值 $\theta_0 1$。def fista(A, b, lam, max_iter1000, tol1e-8, x_initNone): 加速近端梯度法解 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_uvFalse)[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, replaceFalse) 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_iter5000) x_fista, hist_fista fista(A, b, lam, max_iter5000)评估指标用两个目标函数值和相对误差 $|x^{(k)} - x^| / |x^|$。下图的核心结论是FISTA 在 200 步左右就已经收敛到稳定解ISTA 需要 2000 步以上。两者的最终解几乎一致相对误差小于 $10^{-6}$说明收敛到了同一个最优点。稀疏模式恢复准确非零元位置与真解完全一致幅值偏差在噪声水平内。更直观的对比数据迭代次数到某个误差阈值所需量级误差阈值ISTA 所需迭代FISTA 所需迭代$10^{-2}$约 400约 30$10^{-4}$约 1900约 90$10^{-6}$约 4600约 220这个差距就是加速项的威力。实际业务中如果矩阵维度过千上万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_uvFalse)[0] ** 2 # 方法二幂迭代法只算最大特征值大矩阵下更快 def power_iteration(A, num_iter100): 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, axis0) A_norm A / col_norms此时目标函数会有点变化等价于对 $x$ 做了缩放但稀疏模式不变。列归一化后 $L$ 的上界更好估计数值也更稳定。对实测中测量矩阵列尺度差异极大的情况这步基本是必须的。如果不想要固定步长可以上回溯近端梯度backtracking proximal gradientdef proximal_gradient_backtracking(A, b, lam, max_iter1000, tol1e-8): m, n A.shape x np.zeros(n) L np.linalg.svd(A, compute_uvFalse)[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 x4.3 停止准则别只看目标函数很多初学者判断是否收敛只看目标函数值不再下降。这样做很容易被平坦区骗到目标函数在最优解附近确实很平差一点数值差很小但解本身可能离最优解还很远。我建议用两套指标结合相对变化量$|x^{(k1)} - x^{(k)}|_\infty$ 小于某个阈值比如 $10^{-8}$。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. 常见问题速查表与排错思路我把做近端梯度法过程中最常见的问题整理成了一张速查表方便遇到问题快速定位现象可能原因排查方法迭代不收敛/目标函数震荡步长过大L 低估增大 L 或加回溯线搜索FISTA 加重启收敛极慢L 估算过大、步长过小用幂迭代重新估计 L检查是否列归一化解全变成零$\lambda$ 太大缩小 lambda参考 $\lambda_{\max}$ 边界解稀疏性不足$\lambda$ 太小增大 lambda结合交叉验证选参软阈值后解有很小的非零噪声阈值设定太小 / 未收敛加强停止准则检查代码是否对绝对值做了对称处理FISTA 比 ISTA 还慢动量参数设置不当/振荡加重启检查步长是否稳定目标函数平滑但解不对很可能停止条件过早改用 KKT 残差停止准则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$ 越接近真正的最优解。实际使用中我把tol1e-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_1L_2$、全变分ADMM 是把问题拆成更小子问题并行求解的利器。但如果你只是想快速解一个常规大小的 LASSOFISTA 是最具性价比的选择代码量不到 30 行理论性质好实践稳。这也是它在优化工具包里出现频率这么高的原因。我在实际调参中还有个体会不要迷信固定步长除非矩阵设计得非常规整。真实数据里的矩阵条件数往往不理想最省心的是配合一个简单的预警机制——每 20 步检查目标函数是否比前 50 步的最小值高出 1%如果高出了就强行缩小步长或重启动量。这个小监控可以在不牺牲效率的前提下避免项目进入半夜跑完发现崩了的窘境。最后再分享一个小技巧近端梯度类算法本质上干的事情就是展开一个稀疏性约束的递归结构。如果你从深度学习视角去看 FISTA 的迭代式会发现它和一个带残差连接的神经网络非常相似——每一层的结构就是线性变换 阈值非线性。这个视角在做 learned ISTA、即把迭代展开成网络结构的时候非常有用。我有一次做图像稀疏重建就是初始化网络参数为 FISTA 的 $A^T$ 和步长网络训练速度比随机初始化快了一个数量级。把这层联系记在心里近端梯度法对你的价值就远不止解一个 LASSO 那么简单了。