很多人第一次碰到近似梯度下降是在某个带L1正则的回归任务里被卡住的时候梯度下降跑了上千轮损失函数就是降不下去解也不稀疏换成次梯度法更离谱收敛慢到让人怀疑人生。后来有人丢过来一句“试试近似梯度下降”整个问题立刻就不一样了。这篇文章就把这个关键工具讲透包括它解决什么问题、核心迭代是怎么设计的、近端算子到底在算什么、以及我自己实际调参踩过的坑。这套方法适合正在学最优化理论、或者手头有稀疏回归/图像复原/低秩恢复这类需求的人。如果你已经熟悉普通梯度下降和凸优化基本概念这篇看起来会非常顺。如果你刚入门我也尽量把每个“为什么”都讲清楚保证跟得上。1. 问题画像与核心思路1.1 复合优化一块光滑一块扎手近似梯度下降处理的典型目标函数长这样F(x) f(x) g(x)其中f是光滑可微的凸函数g是凸但非光滑的函数可能取值正无穷。这类问题在机器学习里遍地都是最典型的就是LASSOmin ½||Ax - b||² λ||x||₁这里的f是平方损失g是L1范数。还有图像去模糊里的全变分正则项、推荐系统里的低秩矩阵恢复核范数、带约束优化里的指示函数全都是一副“光滑非光滑”的双拼结构。为什么刻意把问题拆成两块因为这两块的“性格”完全不同没法用同一种工具处理。f像一块润滑过的铁板梯度信息好用能指导快速下降g却扎手处处有不可导的尖角比如L1在零点就有尖点。你没法直接对g做梯度下降但g的结构往往简单到可以“闭式求解”。近似梯度下降的核心思想就是让f用梯度信息、让g用近端算子各自发挥长处。1.2 为什么朴素梯度下降在这里不好使很多人一开始都是直接对F求次梯度硬跑。次梯度对L1这种函数在不可导点取的是一个集合随便选一个方向往往不是好的下降方向。结果就是收敛速度极慢理论上次梯度法只有O(1/√k)的收敛率对精度要求高时基本不可用而且步长稍微调不好就直接震荡。那你可能会问非光滑也没事我用光滑化近似行不行比如用Huber函数替代L1的尖角。可行但光滑化会引入额外误差而且正则效果往往偏离原始目标。更重要的是我们希望最终解是真正稀疏的而光滑化之后解虽然接近零却很难精确等于零这在特征选择场景里非常难受。近似梯度下降换了个思路不跟非光滑项硬碰硬而是把g的优化过程“拎”出来做近端映射这个映射往往有闭式解。所以它既能精确处理不可导点又能保持跟梯度下降一样漂亮的收敛速度理论收敛率是O(1/k)。1.3 近似梯度下降的直观逻辑这个算法的每次迭代可以理解成“在当前位置附近简化目标函数然后精确求解这个简化问题”。对光滑部分f我们用一个二次代理模型来近似它也就是在x_k附近做线性展开再加一个二次罚项非光滑部分g保持不变。于是每轮要解的子问题变成x_{k1} argmin_y f(x_k) ⟨∇f(x_k), y - x_k⟩ ½t||y - x_k||² g(y)忽略f(x_k)这个常数整理一下就是x_{k1} argmin_y ½||y - (x_k - t∇f(x_k))||² t·g(y)也就是说算法先沿着负梯度方向走一小步得到一个中间点v然后在“不要太偏离v”的约束下最小化g。这个“在约束下最小化g”的操作就是近端算子。2. 近端算子算法的心脏2.1 prox的定义和直觉近端算子定义为prox_{t g}(v) argmin_x t·g(x) ½||x - v||²t在这里是步长也看成惩罚强度。这个定义看起来抽象但直觉特别简单它找的是“一个尽量靠近v、同时让g尽量小”的点。½||x - v||²是“不要走太远”的约束t·g(x)是“让g变小”的动力两者互相平衡。生活化类比你在商场里逛街v是你此刻的位置g是“远离购物中心”的冲动prox就是你在“原地别动”和“离开”之间取的折中。t越大越想离开折中点离v越远t越小折中点离v越近。这个算子之所以是算法的心脏是因为它把非光滑性全部吸收了。迭代里所有不好处理的部分都被塞进prox里而prox对于很多常见g都有闭式解计算极其便宜。2.2 一维例子软阈值算子怎么来以g(x) λ|x|为例。此时prox_{t g}(v)要最小化tλ|x| ½(x - v)²这是个一维凸问题可以直接求。当x 0时目标变tλx ½(x-v)²最优x v - tλ但前提是v - tλ 0也就是v tλ。当x 0时目标变-tλx ½(x-v)²最优x v tλ前提是v -tλ。当|v| ≤ tλ时最优解落在x 0。合并三种情况就得到著名的软阈值算子S_{tλ}(v) sign(v)·max(|v| - tλ, 0)这个公式你在很多论文里见过。它干的事情特别直白把v向零的方向收缩tλ这么长如果v本来就不够“劲儿”直接一步缩到零。这正是L1正则能产生稀疏解的根本原因——每次迭代都在“砍”幅度小的坐标。2.3 常见prox速查表实战中高频使用的prox整理如下g(x)prox_{t g}(v)说明λ|x|₁sign(v)·max(|v| - tλ, 0)软阈值稀疏诱导核心(λ/2)|x|²v / (1 tλ)简单缩放指示函数I_C(x)C为凸集投影到C即P_C(v)投影梯度下降由此而来非负约束I_{x≥0}max(v, 0)逐坐标截断单位球约束v / max(1, |v|)靠近原点的投影核范数λ||X||_*对矩阵奇异值做软阈值低秩恢复的关键注意最后一行“对奇异值做软阈值”是近似梯度下降能处理矩阵问题的基础。矩阵版本就是把软阈值作用到奇异值分解后得到的奇异值上。这也是为什么这套方法从向量稀疏扩展到低秩矩阵恢复时几乎零成本迁移。3. 迭代格式与收敛性逻辑3.1 从代理函数到迭代格式整个算法写出来只有三行v x_k - t∇f(x_k) x_{k1} prox_{t g}(v)如果g 0prox就是恒等映射算法直接退化成普通梯度下降。所以它叫“近似梯度下降”其实挺贴切它就是梯度下降的一个推广用近端算子处理额外非光滑项。上面的推导还能看出关键设计f的梯度只走一步g的处理则完全由prox完成。每轮迭代不需要求解任何复杂子问题成本主要就是一一次梯度计算加一次prox计算跟普通梯度下降几乎一个量级。3.2 梯度映射与停机判据普通梯度下降有一个自然的停机方式看梯度范数是不是接近零。近似梯度下降没有全局梯度但可以定义“梯度映射”G_t(x) (x - prox_{t g}(x - t∇f(x)))/t当g 0时G_t(x)就是∇f(x)。而且有一个漂亮的性质x是最优点当且仅当G_t(x) 0等价于不动点方程x* prox_{t g}(x* - t∇f(x*))所以停机判据可以直接用||G_t(x)|| ≤ ε这比单纯看目标函数变化可靠得多。目标函数下降有时会因为步长太小而显得“收敛了”但梯度映射的范数会诚实地告诉你离最优解还有多远。我自己写代码时日志里同时输出目标值和梯度映射范数排查问题时能少掉一半头发。3.3 步长选择固定步长与回溯线搜索收敛性分析的标准结果是当t ≤ 1/L其中L是∇f的Lipschitz常数算法目标值单调不增且满足F(x_k) - F(x*) ≤ O(1/k)这里的L怎么来对二次损失f(x) ½||Ax - b||²L就等于矩阵A的最大奇异值平方也就是AᵀA的最大特征值。实际写代码时可以用幂迭代法或直接np.linalg.norm(A, 2)**2算一次就行之后整个迭代都用固定步长。但有些问题的L不好估计或者算一次太贵这时可以用回溯线搜索初始给一个较大的t每步检查代理函数不等式是否成立F(x_{k1}) ≤ f(x_k) ⟨∇f(x_k), x_{k1} - x_k⟩ ½t||x_{k1} - x_k||² g(x_{k1})如果不成立说明步长太大把t缩小比如乘0.8重新算prox直到不等式成立。这个方案额外成本是一次目标函数计算对于很多实际问题反而比精确估计L更省事。我见过不少开源库干脆两种都支持默认用回溯因为稳。4. 加速版本与扩展4.1 FISTA的动量设计朴素近似梯度下降也叫ISTAIterative Soft-Thresholding Algorithm有O(1/k)的收敛率听着还行但实际跑起来特别在条件数很大的问题上还是慢。加速版本FISTA或者说Nesterov动量版把收敛率提升到O(1/k²)。改动只有一行v x_k (k-1)/(k2)·(x_k - x_{k-1}) x_{k1} prox_{t g}(v - t∇f(v))注意梯度在v上算而不是在x_k上算动量系数随时间变化。这个技巧听着简单却是让收敛率变平方级别的核心。FISTA是我实际项目里最常用的版本没有之一。有一个现象值得注意FISTA的目标函数值在小迭代内经常出现锯齿形波动偶尔还会比ISTA高这完全正常。它保证的是“平均意义”上的快速收敛不是单步单调下降。如果页面看时不理解这一点很容易误判成“算法坏了”。4.2 投影梯度下降是特例把g取成某个凸集C的指示函数也就是“x不在C里就罚正无穷”那么prox就变成了投影到C。于是迭代变成x_{k1} P_C(x_k - t∇f(x_k))这就是投影梯度下降。所以投影梯度下降和近似梯度下降不是两个无关的算法而是同一个框架下的两种特例。理解这一点特别重要。比如你在看某篇论文里的“带约束的PGD”本质上就是在跑一个近端算法只是它的prox恰好长成投影的样子。4.3 什么时候选哪种我的经验大致如下目标函数是光滑简单非光滑prox有闭式解直接用FISTA默认省心。约束只是一个凸集投影用投影梯度下降本质上同一套逻辑。prox没有闭式解、需要内层迭代谨慎。每步外层迭代都套一个内层优化成本会放大很多这时可能更适合ADMM或其他分裂方法。问题规模极大、精度要求高可以考虑用随机近端梯度Stochastic Proximal Gradient或者把FISTA与方差缩减技术组合。另外实践里可以给FISTA加“重启”机制每当目标值连续上升超过一定阈值就把动量系数清零让算法退一步重新加速。这在很多病态问题上效果立竿见影。5. 实操从零实现一个LASSO求解器5.1 Python实现ISTA和FISTA不用任何优化库用numpy就可以写一个能用的LASSO求解器。下面是我一直沿用的模板越简单越不容易出错。import numpy as np def soft_threshold(v, thr): # 软阈值算子 S_{thr}(v) return np.sign(v) * np.maximum(np.abs(v) - thr, 0.0) def ista_lasso(A, b, lam, max_iter1000, tol1e-8): m, n A.shape # Lipschitz常数 L ||A||_2^2 L np.linalg.norm(A, 2) ** 2 if L 0: L 1e-12 t 1.0 / L x np.zeros(n) for k in range(max_iter): grad A.T (A x - b) x_new soft_threshold(x - t * grad, t * lam) if np.linalg.norm(x_new - x) tol * max(1.0, np.linalg.norm(x)): x x_new break x x_new return x def fista_lasso(A, b, lam, max_iter1000, tol1e-8): m, n A.shape L np.linalg.norm(A, 2) ** 2 t 1.0 / L x np.zeros(n) y np.zeros(n) k 0 for k in range(max_iter): grad A.T (A y - b) x_new soft_threshold(y - t * grad, t * lam) # Nesterov动量 beta k / (k 3) # 或者另一种更常用的写法beta (k - 1) / (k 2) y x_new beta * (x_new - x) if np.linalg.norm(x_new - x) tol * max(1.0, np.linalg.norm(x)): x x_new break x x_new return x注意L的计算对LASSO里的f(x) ½||Ax - b||²梯度是Aᵀ(Ax - b)其Lipschitz常数确实是||A||₂²也就是最大奇异值的平方。用np.linalg.norm(A, 2)直接算对于中等规模矩阵没问题如果A特别大且是稀疏的改成幂迭代法更划算。5.2 求解稀疏回归的完整流程下面用一个完全合成的数据演示完整流程方便验证。np.random.seed(0) m, n 100, 200 # 低维样本高维特征 A np.random.randn(m, n) x_true np.zeros(n) x_true[:10] np.random.randn(10) # 真正的稀疏解只有10个非零 b A x_true 0.01 * np.random.randn(m) lam 0.05 x_ista ista_lasso(A, b, lam) x_fista fista_lasso(A, b, lam) print(非零元素个数:, np.sum(np.abs(x_fista) 1e-6)) print(恢复误差:, np.linalg.norm(x_fista - x_true))跑完你会发现FISTA和ISTA最终会收敛到很接近的解但FISTA到达相同精度所需的迭代次数明显更少。在更大的数据上差距会越拉越大尤其是条件数差的时候。5.3 参数调试心得\lambda 怎么选这是实践中最大的一个问题。没有“标准答案”但有几个非常实用的经验从大到小试。先给一个很大的λ解全零然后逐步缩小直到出现少量非零元素。这样可以用“解的结构变化”确定数量级。根据理论当λ ≥ ||Aᵀb||∞时LASSO解一定是全零。这个上界可以直接用来定搜索上限省得盲调。如果想复现某篇论文的实验最好沿用他们论文里的λ或者用交叉验证选最优λ。迭代步数方面别一味追求小tol。LASSO这种问题tol设1e-8在工程上几乎意义不大因为你后续还会基于解做阈值化。我一般设1e-6同时限制最大迭代次数两个条件谁先触发都停下防止死循环。6. 踩坑记录与排查技巧6.1 目标函数不降反升最常见的原因就是步长过大或者L估计得偏小。如果你懒得算L直接用回溯线搜索如果你用了固定步长检查一下t ≤ 1/L是否真的满足。还有一个很容易忽略的细节对L1正则t和λ永远是成对出现的代码里写soft_threshold(x - t * grad, t * lam)时两者要同时乘步长别只乘一个。另一个隐蔽原因是数据本身没做标准化。LASSO对特征的尺度极其敏感A的列尺度差异大时不同坐标的梯度量级差太多同一个步长往往顾此失彼。先对A每一列做归一化解出来后再把系数换算回原始尺度能省掉大量调参时间。6.2 解不稀疏辛苦跑完结果全是小非零值、没有真正的零这种情况多半是没收敛到最优。近端算法产生的解理论上是很“干净”的软阈值会把幅度小于阈值的坐标精确置零但你提前停机了就还会有小尾巴。解决方法是提高迭代次数或者用梯度映射范数做停机判据。另外检查一下停止判据如果用的是x_new和x的差注意tol的缩放方式别在某个坐标天然很大的时候直接把所有坐标都判定为“没变化”。还有一种情况λ太小模型本身就该稠密。这不算算法问题是正则化强度选择问题。用交叉验证找λ就好。6.3 加速版本震荡怎么办FISTA在强凸但不均匀的问题上会震荡这是正常现象。你要是实在看不惯或者需要监控目标值单调变化就启用重启机制。具体做法是当连续迭代中出现F(x_k) F(x_{k-1})时把动量项清零y重新设成当前的x_k。我实际用下来重启能显著提高稳定性而且几乎不损失收敛速度。在某些病态问题上甚至比不重启还快。所以现在我的默认配置就是FISTA加重启。6.4 工程化细节最后说几个写生产代码时容易忽略的点对大矩阵不要用np.linalg.norm(A, 2)直接算L改成幂迭代法或者干脆用回溯线搜索省掉L的计算。计算Aᵀ(Ax - b)时不要先算Aᵀb再算AᵀAx后者容易放大浮点误差。固定场景下可以先把Aᵀb缓存起来每次只算Aᵀ(Ax)。如果数据规模大且稀疏用scipy.sparse的矩阵存储配合稀疏矩阵乘法速度能快好几个数量级。记录目标函数值和梯度映射范数的历史存成log文件。出了问题回看log比口头猜原因高效得多。我个人做这类算法调优的体会是永远先验证“步长够不够保守”再谈“收敛快不快”。很多看似玄学的震荡和发散归根结底都是步长超了Lipschitz界限。把步长这一步控制住了近似梯度下降基本上就是一个不太容易翻车的算法。后来有段时间我去做设备端压缩模型里的稀疏训练用的还是这套近端框架只是把A换成了神经网络的部分参数g换成了结构化稀疏约束——内核一点没变。学的时候多花半小时理解prox这个东西后面能帮你在很广阔的问题里省下大量时间。