凸优化算法实践：自动微分、牛顿法与拟牛顿法的工程实现

简介：面向北京大学信息科学技术学院智能科学与技术系“凸分析与优化方法”课程的编程作业合集，适用于正在系统学习自动微分、梯度类优化算法与神经网络参数调优的学生和研究者。压缩包共63个文件，以23个Python脚本为核心，包含自动微分模块、各类优化算法实现及测试代码，并配以32张收敛曲线与迭代过程图像、3个JSON配置、1份说明文件和1份附赠资源文档，整体仅1.03MB，便于快速下载与本地运行。已有92人学习浏览。内容覆盖自动微分、梯度下降、最陡下降、阻尼牛顿法、高斯牛顿法、共轭梯度法、DFP/BFGS/L-BFGS拟牛顿法，以及LADMAP、BCD等进阶方法，每个作业文件夹内都有实验对比图和可执行脚本，能够帮助读者直观理解不同优化器的收敛行为与适用场景。通过动手复现这些算法，能显著提升数值优化实现能力，并为后续科研或工程应用打下扎实基础。

1. 凸分析与优化方法：这份代码集能直接跑出牛顿法、共轭梯度与拟牛顿法的收敛曲线

北京大学信息科学技术学院智能科学与技术系的凸分析与优化方法课程，编程作业难度一直不低。这份作业集合把从自动微分、梯度下降、回溯线搜索，到阻尼牛顿法、高斯牛顿法、共轭梯度法、DFP、BFGS、L-BFGS，再到块坐标下降与随机梯度法全部串在一起，每个 HW 目录下都有可运行的 Python 源文件和对应的收敛曲线图。对正在啃凸优化教材的研究生来说，它是现成的算法骨架；对做机器学习参数调优的工程师来说，它能让你把“梯度下降不够快”换成“试试阻尼牛顿或 BFGS”。这份代码的价值在于，每类算法的实现都不是课本伪代码的直译，而是带线搜索、带收敛性检查、能跑出可视化曲线的工程化版本。

2. 自动微分到阻尼牛顿法：从反向传播到底层优化，这一章最关键

2.1 autoDiff.py 的自动微分实现逻辑

HW3 里的 autoDiff.py 和 grad4Net.py 是整套作业的地基。自动微分要解决的核心问题是：给定一个由基本运算组成的复合函数，如何精确计算它对每个输入的偏导。常见做法是构建计算图，每个节点记录操作类型和子节点，前向传播时保存中间结果，反向传播时利用链式法则逐层求导。课程作业里没有直接调用 PyTorch 的 autograd，而是手写了一个最小实现，这样你才能真正看到梯度是怎么算出来的。

# 参考 HW3 常见实现框架：节点定义与前向反向 class Node: def __init__(self, value, children=(), op=None): self.value = value self.children = children self.op = op self.grad = 0.0 def backward(self, grad=1.0): self.grad += grad if self.op == 'mul': left, right = self.children left.backward(grad * right.value) right.backward(grad * left.value) elif self.op == 'add': for child in self.children: child.backward(grad) elif self.op == 'sigmoid': # 链式法则：sigmoid'=sigmoid*(1-sigmoid) (child,) = self.children child.backward(grad * self.value * (1.0 - self.value)) else: raise TypeError(f"unknown op: {self.op}")

逻辑说明：每个 Node 持有一个 forward 之后的值，减少重复计算；backward 递归发放梯度到子节点，父节点梯度累加到 self.grad，这对应多路径共享的自动微分语义。grad4Net.py 把这个机制接到一个两层全连接网络上，通过 grad4Net.json 描述网络结构，用自动微分计算损失对每个权重的梯度，再交给梯度下降去迭代。关键参数是网络层数与激活函数类型，HW3 的测试脚本 grad4NetTest.py 会对比数值差分，验证梯度计算的正确性。

2.2 梯度下降与回溯线搜索：L_alpha.png 和 L_beta.png 蕴含步长密码

HW11 的 gradDescent.py 是标准梯度下降，HW12 的 steepDescent.py 是最陡下降。两者的区别在于最陡下降每一步沿着负梯度方向，并且通常配合精确线搜索或回溯线搜索。作业图里的 L_alpha.png、L_beta.png 是经典的可视化：横轴分别是步长 alpha 和衰减因子 beta，纵轴是目标函数值，用来标记线搜索的行为。

def backtracking_line_search(f, df, x, direction, max_alpha=1.0, beta=0.7, c=1e-4): # f：目标函数，df：梯度函数，direction：下降方向，通常取负梯度或牛顿方向 alpha = max_alpha grad_norm = df(x).dot(direction) f_x = f(x) while f(x + alpha * direction) > f_x + c * alpha * grad_norm: alpha *= beta if alpha < 1e-10: break return alpha

逻辑说明：Armijo 条件要求目标函数下降量至少要达到 c * alpha * 梯度方向内积的预期，不满足就把步长乘 beta 缩小，重复直到满足。c 通常取 1e-4，beta 常见取值 0.5 到 0.8，beta 越大搜索越精细但每步检查次数越多。固定步长的梯度下降在病态目标上极易震荡，而回溯线搜索把每个迭代步的 alpha 用同样的准则动态调整，这也是最陡下降法优于朴素梯度下降的关键。HW12 的 converge.png 对比了这两种方法的收敛曲线，典型的图是固定步长在峡谷目标上左右反弹，回溯线搜索则平滑下降。

2.3 阻尼牛顿法：Newton.py 里的参数调节机制

阻尼牛顿法在纯牛顿法的基础上加了一个阻尼项，常见形式是把迭代写成 x_{k+1} = x_k - alpha_k * (H_k + lambda_k I)^{-1} * g_k，其中 lambda_k 就是阻尼系数。小 lambda 接近纯牛顿法，收敛快但容易震荡；大 lambda 逼近梯度下降方向，鲁棒性增强。作业里 converge_damped.png 展示的是阻尼牛顿法在不同 lambda 下的迭代曲线，有一个明显的观察点：lambda 从 0.01 调到 1 后，前期收敛速度虽然慢了一些，但整体稳定。

# HW12 Newton.py 常见实现：阻尼牛顿法 import numpy as np def damped_newton(f, grad, hess, x0, max_iter=100, tol=1e-6, damp=0.1): x = x0.copy() history = [] for _ in range(max_iter): g = grad(x) H = hess(x) dx = np.linalg.solve(H + damp * np.eye(len(x)), -g) x = x + dx history.append(f(x)) if np.linalg.norm(g) < tol: break return x, np.array(history)

参数说明：damp 是阻尼系数，np.eye 加到 Hessian 矩阵上保证求解线性系统时矩阵可逆。这里也可以换成回溯线搜索找步长，而不是直接使用 dx 全长；作业里通常两种情况都实现，对比纯阻尼牛顿和阻尼牛顿加线搜索的差异。高斯牛顿法（HW12 的 converge_gauss.png）本质上是针对非线性最小二乘的近似牛顿法，用 J^T J 近似 Hessian，适用于残差较小的优化问题，代码上就是把 hess 换成残差的雅可比内积，其余结构完全一致。

3. 拟牛顿法与共轭梯度：DFP、BFGS、L-BFGS、conjGrad 到底差多少

3.1 DFP 与 BFGS：公式只差两行，行为差了一个时代

HW13 的 DFP.py 和 HW14 的 BFGS.py 都解决同一个问题：不想显式计算 Hessian 矩阵，而是用相邻两次迭代的梯度变化和位移量来近似。DFP 是对拟牛顿矩阵 B_k 做秩二校正，BFGS 则直接对近似 Hessian 的逆矩阵做秩二更新。两者的更新公式很接近，但 DFP 在处理某些非凸目标时容易产生数值不稳定的逆近似，BFGS 在多数机器学习问题上表现更稳健。HW14 里同时出现 BFGS.py 和 L-BFGS.py，正是为了让学习者对比完整拟牛顿与内存受限版本。

# BFGS 常见更新的核心片段 def bfgs_update(H_inv, s, y): rho = 1.0 / (y.dot(s)) Hy = H_inv.dot(y) # 单秩二更新：新 H_inv = (I - rho*s*y^T) H_inv (I - rho*y*s^T) + rho*s*s^T return H_inv + rho * np.outer(s, s) - (np.outer(Hy, s) + np.outer(s, Hy)) * rho

逻辑说明：s 是参数增量 x_{k+1} - x_k，y 是梯度增量 g_{k+1} - g_k。更新必须满足曲率条件 y.dot(s) > 0，不满足时说明步长没有让目标充分下降，本次迭代跳过更新是常见保护策略。BFGS.py 用线性代数库直接维护稠密矩阵，参数量在几千以下时没问题；一旦到几万参数，存不下也只能放下。L-BFGS.py 的策略是只保留最近 m 组 s、y 对，默认 m=10 到 20，计算方向时用双循环递归，不需要显式构造矩阵。

3.2 共轭梯度法：conjGrad.py 在大规模问题上的特殊价值

HW13 里的 conjGrad.py 实现的共轭梯度法，最核心的性质是：对 n 维二次目标，理论上至多 n 步收敛。它不存储 Hessian，只依赖前一个搜索方向 p_{k-1} 和当前负梯度的线性组合，组合系数有 Fletcher-Reeves 和 Polak-Ribiere 两种计算方式。课程作业通常两者都写，对比它们在不同条件下的收敛轮数。

def conjugate_gradient(f, grad, x0, max_iter=100, tol=1e-6, use_pr=False): x = x0.copy() g0 = grad(x) p = -g0 for _ in range(max_iter): # 步长：二次问题上一步除出精确解，这里用线搜索近似 alpha = backtracking_line_search(f, grad, x, p) x = x + alpha * p g1 = grad(x) beta = g1.dot(g1) / g0.dot(g0) # Fletcher-Reeves if use_pr: beta = g1.dot(g1 - g0) / g0.dot(g0) # Polak-Ribiere p = -g1 + beta * p g0 = g1 if np.linalg.norm(g0) < tol: break return x

参数说明：max_iter 通常设为参数维度数，不要在迭代到 n 步后才考虑重启；对于非二次目标，共轭方向会退化，胡乱继续不如直接重置为负梯度方向，这就是所谓“每 n 步重启”的工程习惯。conjGrad.py 里可视化曲线往往显示，在中等规模问题上，它用比 BFGS 更少的内存，达到相近的收敛精度。

3.3 拟牛顿法选型对比

HW17 的 compare.py 专门做这种对比实验，把同一目标函数丢给不同求解器，记录每轮函数值，出 compare.png。这个图能直观看到，BFGS 和 L-BFGS 曲线几乎重合，DFP 掉在后面，共轭梯度在中途有几次重启的跳变。这种差异不是理论上的“谁更优”，而是数值实现细节决定的——所以才有必要自己跑一遍。

4. 约束优化与高级变体：惩罚法、Frank-Wolfe、块坐标下降、随机梯度

4.1 惩罚函数法与 Frank-Wolfe：两种处理约束的思路

HW18 的 penalty.py 把带约束问题转成带惩罚项的无约束问题，常见形式是目标函数加上 sigma * max(0, g(x))^2 或者等式约束的平方项，然后对增广目标做无约束优化，再逐步增大 sigma。这个思路简单，但 sigma 逼近无穷时 Hessian 条件数恶化，数值上容易翻车。与之对应，FW.py 的 Frank-Wolfe 算法走的是另一条路：不投影回可行域，而是每次在可行域上求解一个线性优化子问题，用凸组合的方式更新点。Frank-Wolfe 特别适合可行域是单纯形、谱范数球这类线性优化好解的结构，比如矩阵补全和稀疏学习。

4.2 块坐标下降：BCD4LR.py 与 BCD4Dict.py 的应用差异

HW19 的 BCD4LR.py 是把块坐标下降用到逻辑回归上。逻辑回归目标函数是凸的，坐标下降每次只更新一个或一组参数，固定其它参数，子问题带正则项时还有解析解。对比实验里，它在特征维度较高的数据集上，收敛速度往往比全批量梯度下降快，因为每次更新只算一个特征的损失变化。BCD4Dict.py 则是字典学习，目标带稀疏约束项，块坐标下降更新字典原子和稀疏编码交替进行。这两份代码放一起，能让读者看清坐标下降的适用面：要么子问题有解析解，要么子问题维度足够小可以快速求解。

4.3 随机梯度、加速梯度与增广惩罚：HW20、HW21 的现代武器

HW20 的 Stochastic.py 实现的是随机梯度下降与随机坐标下降的变体，它对应大数据场景，每次只用一个样本估计梯度，方向噪声很大但计算成本低。课程作业里通常会画 1.png 这类迭代曲线，随机梯度下降的曲线波动明显，但整体下降趋势比全批量更早出现，这正好反映“牺牲精度换速度”的工程取舍。HW21 的 AGD.py 是 Nesterov 加速梯度，它只在前向梯度基础上加了一个动量项，却能把凸光滑函数的收敛率从 O(1/k) 提升到 O(1/k^2)；APenalty.py 是增广拉格朗日惩罚法，属于前面 penalty.py 的升级版，引入对偶变量并用交替方向策略逼近约束。

5. 避坑与排查：五个血泪经验，照着检查能少走一整天弯路

5.1 固定步长梯度下降迟迟不收敛

现象：HW11 的 gradDescent.py 跑出的 converge.png 不是一条平稳下降曲线，而是一条锯齿形折线，500 次迭代后函数值还在同一个量级震荡。原因：固定步长在目标函数曲率方向差异大时，沿着陡峭方向的步长偏大，峡谷壁上来回弹；减步长能稳但慢，加步长直接发散，这就是梯度下降对步长极其敏感的根本原因。解决：换成回溯线搜索，把步长决策交给 Armijo 条件，HW12 里 steepDescent.py 就是干这个的。

5.2 阻尼牛顿法前期迭代反而变慢

现象：converge_damped.png 上，阻尼系数取 0.5 的曲线比取 0.01 的曲线在前 20 轮下降明显更慢，一度怀疑是自己实现错。原因：阻尼项在分母上削弱了二阶信息的作用，lambda 越大，迭代方向越接近纯梯度下降，接近牛顿法的高收敛阶被抑制。解决：阻尼系数不要从一开始就固定，建议先用小阻尼比如 0.001 跑 5 轮，如果发现函数值反弹或震荡，再逐渐增大；也可以用回溯线搜索替代固定阻尼，让步长自适应收缩。

5.3 L-BFGS 换数据集后收敛速度骤降

现象：同一套 L-BFGS.py 代码，在 A 数据集上 30 轮收敛，换到 B 数据集后 200 轮还在缓慢下降。原因：L-BFGS 保留了历史 s、y 对来近似当前 Hessian，换数据后旧梯度信息不再代表当前目标曲率，近似矩阵与实际 Hessian 严重不符。解决：切换数据集或任务时强制清空 L-BFGS 的历史缓冲区，从上一次迭代点重新积累曲率对。作业框架里通常没有自动 reset，需要自己在训练循环开头手动调 LBFGS.reset()。

5.4 共轭梯度法不重启导致收敛崩溃

现象：conjGrad.py 在多项式回归目标上跑到第 20 轮，损失突然飙升，然后恢复，反复横跳。原因：共轭梯度仅对二次目标能维持共轭方向，非二次目标上每一轮方向之间共轭性被破坏，继续沿用旧方向会引入误差，数值上积累到一定程度就产生方向退化。解决：每最多 n 步重置一次，n 为参数维度；同时用 Polak-Ribiere 的 beta 计算公式，当 g_{k+1} 与 g_k 近似正交时，beta 自动接近零，产生隐式重置效果。

5.5 数值差分验证自动微分梯度时误差无法收敛

现象：grad4NetTest.py 对比自动微分梯度和有限差分梯度，误差始终停在 1e-4 不往下走。原因：有限差分步长 h 固定取 1e-6 时，截断误差与舍入误差有冲突；h 太大一阶差分误差大，h 太小浮点舍入吃掉有效位。解决：采用中心差分 (f(x+h)-f(x-h))/(2h)，并把 h 设置在 1e-4 到 1e-6 之间分段测试；如果两种步长得到的误差量级变化剧烈，说明目标函数在某些点上不可微或者计算图中存在数值不稳定的操作。

6. 把作业里的收敛图当作验证工具：一种高效的算法自查方式

看完上面这些实现，最实用的一个技巧是把作业里附带的 PNG 图当成“标记图像”，用来验证你自己重写的算法是否正确。这种方法的原理很简单：这些图是老师提供的地面真值实现跑出来的结果，只要你的算法在同样的初始点、同样的目标函数、同样的迭代次数下，函数值序列能贴合原图，你的实现就没大问题。但如果你的最优值差了一大截，或者下降速度明显不对，通常说明代码里某个细节写错了，这是最直观的自查工具。

import numpy as np import matplotlib.pyplot as plt # 读取作业收敛数据，比如自己记录的函数值 history # 假设 history 是 (iters, values) 的 npy def plot_convergence(history, reference_img="converge.png", label="mine"): plt.plot(history[:, 0], history[:, 1], label=label, linewidth=2) plt.xlabel("iteration") plt.ylabel("objective value") plt.yscale("log") plt.grid(True, alpha=0.3) plt.legend() # 把作业原图叠加成水印对比，可以在 jupyter 里用 Image 打开参考图后手动对齐 plt.savefig(f"compare_{label}.png", dpi=150)

参数说明：plot_convergence 要求历史数组第一列是迭代序号，第二列是目标值；横轴用线性刻度，纵轴用对数刻度，因为大部分优化算法在前几轮下降很快，后面进展缓慢，线性刻度会把后期细节压扁。用对数刻度能更清楚看到收敛率是否匹配。如果你在自己的复现工程里跑算法，建议每次迭代记录三个指标：函数值、梯度范数、参数更新量的范数，三者一起画在同一张图上。函数值不降但梯度范数很小，说明目标本身需要检查；函数值下降但梯度范数不降，则可能是线搜索方向有问题。

具体操作上，我一般会把作业里的 converge.png、converge_damped.png、converge_gauss.png 全部放在一个 reference 文件夹里，自己每重写一个算法就重新跑一遍，生成的图直接和原图叠加对比。如果曲线形状走向一致，只是迭代轮数略有差别，那是线搜索参数不同导致的，属于正常范围；如果形状都变了，就要回去查更新公式。从那以后，我每次拿到新的优化器定义，都强制自己先画一条收敛曲线再谈精度，这个习惯帮我避开了很多“代码看着对、结果全是坑”的状况。希望帮到你，这套作业源码值得按上面思路亲手跑一遍。

本文还有配套的精品资源，点击获取