研究生阶段的优化方法课，上机作业一直是个让人又爱又恨的存在。这学期拿到作业要求的时候，我反复看了三遍才确认，不是简单调个现成优化库就能交差的——从最速下降法、牛顿法、共轭梯度法到BFGS，再到带约束的罚函数法，全部要求手写实现，还要在标准测试函数上对比收敛行为。更关键的是，老师会在提交后随机抽查代码细节，问你为什么这里用Armijo回溯而不是精确线搜索，罚因子序列为什么要按这个比例递增。这门课的上机作业，本质上不是在考你会不会用工具，而是在逼你把每种优化算法背后的推导和数值行为彻底吃透。

这篇文章把我这次上机作业从读题、设计、编码到调试对比的完整过程记录下来，包括算法实现的核心细节、参数选择的依据，以及我在实际运行中踩过的坑。不管是正在修这门课、需要应付类似作业，还是单纯想搞懂这些经典优化算法到底怎么落地，这篇内容都应该能给你省下不少时间。

1. 上机作业的核心考点与整体设计思路

1.1 作业要求到底在考察什么

很多同学刚看到题目时会觉得，"优化方法"嘛，无非就是拿scipy.optimize.minimize把函数丢进去，设置一下method参数就完事了。但真正动手做才发现，作业要求写得很清楚：不允许直接调用现成的优化算法库，所有迭代过程都要自己实现。这意味着你得知道每一步迭代在干什么，为什么这样更新，终止条件怎么定，以及算法在什么情况下会失效。

从考察点来看，作业其实拆成了三个层次。第一层是基本实现能力，比如最速下降法的梯度方向、牛顿法的Hessian矩阵修正、共轭梯度法里β系数的两种常见公式能不能写对；第二层是数值实验能力，要求在同一组测试函数上对比不同算法的迭代次数、收敛精度和稳定性；第三层则是约束优化部分，通常需要你用外点罚函数法或增广拉格朗日法把带约束问题转成无约束问题，再观察罚因子变化对最优解和计算过程的影响。这三个层次分别对应了算法理解、工程实现和数值分析三种能力，也是为什么这份作业在研究生课程里一直保留至今的原因。

我自己在动手之前做了一个拆分表格，把每种算法的输入输出、迭代格式、终止条件、最坏情况全部列出来，再对照作业要求逐项打钩。这样做的好处是，编码的时候不会漏掉关键步骤，调试时也能快速定位是算法写错了还是参数没调对。

1.2 为什么选Python加NumPy而不是C++或MATLAB

选实现语言这件事，表面上是个人偏好，实际上会影响你整个上机作业的推进速度。我最后选了Python加NumPy，核心原因有三个。

第一，向量化计算让梯度表达式非常简洁。像最速下降法里x = x - alpha * grad这种更新，在NumPy里就是一行的事，而且不管变量维度是二维还是五十维，代码结构完全不用改。第二，画收敛曲线直接用matplotlib，几行代码就能生成目标函数值随迭代次数下降的图，这对实验分析太重要了。第三，Python的交互式环境方便边写边测，我可以把每个算法封装成函数，再写一个统一的实验脚本来调用，出了问题不用重新编译。

当然Python也有缺点，比如循环迭代慢。但这个作业里测试函数规模通常不大，维度一般控制在几十维以内，就感受不到性能瓶颈。如果你非要用C++去实现，那大部分时间可能都耗在矩阵运算和内存管理上，反而偏离了作业考察的核心。对我们这个场景来说，Python的快速验证价值远远大于执行效率上的损失。

代码组织上我也走了点捷径，但这条捷径值得推荐：用一个test_functions.py文件存放Rosenbrock函数、二次函数、多峰函数及其梯度；一个algorithms.py存放所有优化器的实现；一个experiment.py负责跑实验并输出结果表格和收敛图。这样每个文件职责单一，后面改参数、加对比算法都很顺手。

2. 无约束优化算法：从最速下降到BFGS的实现细节

2.1 最速下降法：步长选择不简单

最速下降法是最容易写出来的算法，迭代公式就一句话：x_new = x - alpha * grad。但真让它在Rosenbrock函数上跑起来，你会发现它简直慢到怀疑人生。原因在于，最速下降的相邻两步搜索方向是正交的，相当于在狭窄山谷里来回折返，形成典型的锯齿现象。所以这个算法在理论上最直观，在数值上却最容易暴露步长策略的问题。

步长选择我采用了Armijo回溯线搜索而不是精确线搜索。精确线搜索虽然听起来更"精确"，但每步都要解一个一维极小化问题，计算量大，而且对于非凸函数反而不稳定。Armijo的思路很简单：先取一个比较大的步长，如果目标函数下降量不满足充分下降条件，就把步长乘以一个小于1的比例缩小，直到条件满足。标准形式是：

def armijo_backtracking(f, x, grad, alpha=1.0, rho=0.5, c1=1e-4): while f(x - alpha * grad) > f(x) - c1 * alpha * np.dot(grad, grad): alpha *= rho if alpha < 1e-12: break return alpha

参数c1通常取1e-4，rho取0.5左右。c1的作用是控制"充分下降"的宽松程度，取太小会让步长容易通过，取太大又会让步长被迫变得过小。我自己实验下来，1e-4是个比较稳的经验值。

这里有个关键点要提醒：如果梯度为零向量，理论上已经到达稳定点，必须提前终止，否则后面会除以零或者步长乱跳。还有就是在循环里要限制最大迭代次数，防止算法在平坦区域无限跑下去。我在最速下降法上设置的是1000次上限，实测在Rosenbrock函数上经常撞到上限还没达到1e-6的梯度精度，这本身就说明最速下降法不适合作为高维强非线性的主力算法，但它作为理解优化过程的起点，价值是不可替代的。

2.2 牛顿法与阻尼牛顿法：二阶信息的双刃剑

牛顿法的迭代格式是x_new = x - H^{-1} * grad，其中H是目标函数的Hessian矩阵。它利用了二阶曲率信息，所以在二次函数上一步就能收敛，在Rosenbrock这类非二次函数上也远比最速下降法快。但好景不长，牛顿法第一个坑就是Hessian矩阵可能奇异或者非正定。

我在写代码的时候用的是阻尼牛顿法，也就是每次更新时加上一个步长系数：

def damped_newton(f, grad, hess, x0, tol=1e-6, max_iter=100): x = np.array(x0, dtype=float) for k in range(max_iter): g = grad(x) if np.linalg.norm(g) < tol: break H = hess(x) # 修正非正定Hessian try: d = np.linalg.solve(H, -g) except np.linalg.LinAlgError: H += 1e-8 * np.eye(len(x)) d = np.linalg.solve(H, -g) # 线搜索，防止步长过大跑飞 alpha = armijo_backtracking(f, x, g, alpha=1.0) x = x + alpha * d return x, k

为什么加线搜索？因为纯牛顿法在非凸区域可能朝着目标函数上升的方向走。你算出的牛顿方向虽然是基于局部二次模型的，但二次模型本身就可能是错的，如果不做步长控制，一步就可能跳到十万八千里之外。我一开始没加线搜索，直接在Rosenbrock函数上跑，某次迭代后目标值从几十直接变成几千万，那场面相当酸爽。

关于Hessian矩阵正定性，实测中最简单的处理方式是加一个小的单位阵修正，也就是所谓的LM修正。更专业的做法是修正Cholesky分解，但作业场景下加正则项已经够用。还有一点值得注意：高阶解析求导很麻烦，如果Hessian是数值差分近似的，一定要控制差分步长。我偷了个懒，对每个测试函数都手推了梯度和Hessian表达式，这样实验更可控，也方便对照算法输出的中间结果判断错误。

2.3 共轭梯度法：方向的记忆与重置

共轭梯度法在第一次迭代时方向取负梯度，之后每次方向由当前负梯度与前一个方向的线性组合构成，系数β有Fletcher-Reeves（FR）和Polak-Ribiere-Polyak（PRP）两种常见形式。FR的公式是beta = (g_new · g_new) / (g · g)，PRP则是beta = (g_new · (g_new - g)) / (g · g)。我两个都写了，实测PRP在非二次函数上表现通常更好，因为它在接近最优值时能自动调整方向，减轻锯齿问题。

非线性共轭梯度一个容易忽略的细节是重置策略。在理论上，β非负且方向必须保证是下降方向。实际写代码时，如果计算出来β为负值，或者当前梯度范数已经很小，就直接把搜索方向重置为负梯度方向：

d = -g_new + max(beta, 0) * d

这里的max(beta, 0)是一个简单而有效的保护。如果你不处理负β，方向可能变得不下降，后面就全乱套了。

收敛速度方面，共轭梯度法在Rosenbrock函数上表现明显优于最速下降法，但不如Newton法。它的优势在于不需要存储和求逆Hessian矩阵，所以在高维问题上内存优势很大。写作业时建议你同时输出每一轮迭代的方向内积np.dot(d, g)，如果这个值长时间为正，说明方向计算有bug，需要立刻回头查β公式。

2.4 BFGS拟牛顿法：Hessian近似的工程选择

BFGS是这次作业里我最喜欢的算法，因为它几乎结合了所有优点：不需要精确Hessian、收敛速度快、数值稳定性好。它的核心是用迭代过程中梯度和位移的信息去近似Hessian的逆，更新公式不复杂但容易写错。

经典BFGS更新公式里，需要定义s = x_new - x和y = g_new - g，然后更新矩阵B：

def bfgs_update(B, s, y): rho = 1.0 / np.dot(y, s) I = np.eye(len(s)) # 注意：如果 y*s 接近0，rho会爆炸，需要保护 term = I - rho * np.outer(s, y) B_new = term @ B @ term.T + rho * np.outer(s, s) return B_new

初始的B一般取单位矩阵，这样第一步相当于最速下降步。BFGS有个很好的性质：只要初始B正定，并且每次线搜索满足Armijo条件，更新后的B始终保持正定。这个性质是它稳定性的基石，所以线搜索不是可有可无的，而是算法修正的一部分。

我在比较实验里统计过，在Rosenbrock函数上，同样从(-1.2, 1.0)出发，达到梯度范数小于1e-6，最速下降法超过500次迭代甚至经常撞上限，共轭梯度法大约需要几十次，BFGS大约二十几次，阻尼牛顿法十次以内就可以完成。这个对比能直观说明，利用二阶或者近似二阶信息对优化效率的提升是非常显著的。

这个表不是绝对标准，因为结果会受初始点和线搜索参数影响，但它能反映各类算法的基本画像。

3. 约束优化实战：罚函数法与增广拉格朗日法的取舍

3.1 外点罚函数法的实现与数值困难

作业的约束优化部分，最基础的要求是用外点罚函数法求解一个带等式约束或不等式约束的最小化问题。核心思路是把约束条件变成一个惩罚项加到目标函数上，然后对无约束问题反复求解。比如等式约束c(x)=0，惩罚形式就是目标函数加上mu * ||c(x)||^2，罚因子mu从小到大递增。

实现时有个非常重要的细节：罚因子的递增策略。如果每次只增加一个很小的倍数，外层循环会收敛得很慢；但如果你一次把罚因子调得特别大，无约束子问题就会变得病态，目标函数里惩罚项占据绝对主导，数值计算精度下降。我在作业里选了mu = 10 * mu的比例，配合从可行/不可行初始点分别做了测试。

另一个问题是初始点的选择。外点罚函数法的名字已经说明了，它允许从不可行点出发，算法自己去把点"拉"回可行域。这既是优点也是隐患：初始点离可行域太远时，惩罚项梯度会非常大，步长控制稍有不慎就会导致数值溢出。我在试验一个二维问题min x1^2 + x2^2, s.t. x1 + x2 = 1时，从(100, -100)出发，罚因子没控制好，目标值直接跑到1e10量级，画出来的收敛曲线惨不忍睹。

所以实操上我建议：外点罚函数法的外层循环不要贪心，罚因子倍增控制在5到20倍之间，内层无约束优化器的精度可以先用1e-4粗跑，再逐步提高到1e-6。这样整体计算量比较均衡，收敛曲线也漂亮。

3.2 内点法与增广拉格朗日法的对比

如果约束是不等式约束g(x) <= 0，常见做法是内点法，也叫障碍函数法。思路是在目标函数后面加上对数障碍项-mu * sum(log(-g(x)))，使得迭代点始终保持在可行域内部。和外点法最大的区别是，内点法的迭代点始终可行，因此很多应用场景更喜欢它。

但内点法的核心难点在于初始点必须严格可行。作业里如果给的初始点刚好卡在边界附近，对数项直接变成无穷，算法直接罢工。我一开始没注意这一点，设置了一个边界上的初始点，反复跑出NaN。解决办法有两个：一个是手动往可行域内部移动初始点；另一个是先用投影或弹性方法求一个可行点，再去跑内点循环。后者实现复杂度高一些，但更通用。

增广拉格朗日法则在罚函数法基础上引入了乘子估计，迭代格式大约是L(x, lambda) = f(x) + mu * ||c(x)||^2 + lambda * c(x)，其中lambda每次按公式更新。这个方法的优势在于罚因子不需要取得特别大，乘子项会承担一部分约束作用，从而减轻病态问题。作业里如果给了加分项，实现增广拉格朗日法是个不错的选择，因为它对数值稳定性的改善是看得见的。我在同个约束问题上对比过，外点罚函数法需要罚因子到1e6才能让约束误差降到1e-5，而增广拉格朗日法罚因子到1e3左右效果就足够了。

3.3 一个带约束的测试案例完整流程

我作业里选的测试案例是min f(x) = x1^2 + x2^2 + x1*x2，约束条件是x1^2 + (x2-1)^2 - 1 <= 0。这个问题的约束是一个圆域，目标函数是二次型，解析解能算出来，方便验证代码正确性。

用外点罚函数法时，惩罚目标函数是：

def penalty_obj(x, mu): g = x[0]**2 + (x[1]-1)**2 - 1 # 不等式约束值 return x[0]**2 + x[1]**2 + x[0]*x[1] + mu * max(0, g)**2

注意这里约束违反量用了max(0, g)，因为只要g <= 0就是可行的，惩罚项应该为0。这个细节不写对，罚函数法就会把可行点也往里面硬拉，结果肯定错。

外层循环我设置如下：

mu = 1 for k in range(50): x_opt = minimize_unconstrained(penalty_obj, mu) violation = max(0, constraint(x_opt)) if violation < 1e-5: break mu *= 10

跑完以后把迭代路径画在可行域图上，可以看到点从圆外逐步被拉进圆周附近，最后停在约束边界上的最优点附近。这张图几乎算是作业的灵魂，它比任何口头描述都能说明算法的行为。

4. 实验对比、参数调节与常见问题排查

4.1 如何设计可复现的实验并分析收敛曲线

上机作业最拉分的地方通常在实验对比部分。老师想看到的不只是代码能跑，而是你会不会系统性地对比算法性能。我建议每个测试函数固定使用同一个初始点、同一个终止条件，然后分别记录迭代次数、目标函数最终值、梯度范数、运行耗时。

收敛曲线画两张图最有说服力。第一张是目标函数值随迭代次数下降的曲线，用semilogy画，因为很多算法后期下降是线性收敛或超线性收敛，对数坐标能清楚看到斜率差异。第二张是梯度范数随迭代次数的变化，梯度范数到0的速度反映了算法接近最优点的速度。这两张图放在一起，基本就能直观区分出线性收敛、超线性收敛和二次收敛的区别。

我自己的实验脚本里有一个统一的run_experiment函数，输入算法对象、测试函数、初始点和精度，输出一个字典。循环跑完所有算法后，直接生成一个DataFrame，再分别画图。这样做的好处是后面你只需要改初始点或精度，所有数据自动更新，不会出现表格和图像对不上的情况。

4.2 上机过程中踩过的典型坑

这一部分是我最想写的，因为每个人在实际调试中遇到的坑其实非常相似，但教科书上从来不会告诉你。

第一个坑是线搜索里目标函数值没有下降。最常见的原因是梯度计算写错了，方向根本不是下降方向。我的排查思路是先单独调用梯度函数，在某个点计算数值梯度和解析梯度对比，确保梯度无误后再跑算法。

第二个坑是牛顿法Hessian奇异。这个问题在作业里几乎必然遇到，因为很多测试函数在满足终止条件之前，某一步Hessian就是不可逆的。我加了一个try ... except组件，捕捉LinAlgError后给Hessian对角线加一个小扰动。这个处理不能太随意，扰动太小没效果，太大又破坏了二阶信息的准确性，一般取1e-6到1e-8量级。

第三个坑是共轭梯度法方向不下降。我在2.3节里提过，如果β计算出来是负的，必须重置方向。但还有一个隐藏问题：你计算β用到的梯度可能是上一轮迭代之前缓存的旧值。我踩过一次这个坑，搞了整整一下午才发现beta公式里用的梯度是旧变量，跟当前梯度差了整整一轮迭代。

第四个坑是罚函数法目标函数被惩罚项完全淹没。尤其是约束函数数量级很大时，目标函数本身的贡献可以忽略不计，导致罚函数的最优解严重偏向约束可行而不关心目标值。解决方法一是缩小初始罚因子，二是约束函数先做归一化处理，三是改用增广拉格朗日法。我最后是通过归一化约束和适度增加罚因子倍数两个手段解决。

4.3 参数调节经验：从精度到调试细节

关于参数调节，我的体会是不要迷信某组固定的数值，而要理解每个参数在算法中的角色。tol是终止容差，一般取1e-4到1e-6。对于作业来说，取1e-5是一个比较均衡的选择：既能展示出收敛趋势，又不会因为精度过高导致迭代次数过多、图像拉不开差距。max_iter设置在500到1000之间足够，但要注意最速下降法在Rosenbrock函数上很可能撞到上限，这在实验报告里反而可以作为讨论点，不必强行让每种算法都收敛到极高精度。

线搜索参数rho和c1是另一个经验集中的地方。rho=0.5是比较保守的衰减，适合大多数测试函数；如果你觉得步长衰减太慢，可以改成0.2，但不要低于0.1，否则每一步步长太小，收敛速度会很慢。c1取1e-4通常没问题，但如果你发现步长频繁被缩小到接近机器精度，可以把c1调小到1e-6，给步长更多通过空间。

调试细节上我推荐一个技巧：在优化器的核心循环里加一个可选的trace参数，每次迭代把当前变量、目标值和梯度范数打印出来。你不需要全量输出，只看最初几轮和最后几轮，就能判断算法是否正常。我曾在测试一个多峰函数时，发现算法第十轮开始来回震荡，就是因为trace输出里看到变量在两个局部极小之间反复跳跃。没有这个输出，我可能需要花更长时间才能定位问题。

5. 实操中的额外体会与调试技巧

5.1 调试优化算法的通用顺序

踩了这么多坑之后，我总结出一个调试优化算法的通用顺序，分享出来应该对后面做类似作业的同学有帮助。第一步永远是最小测试：用一个一元函数或者一个二维二次函数跑算法，因为这类问题你能手动算出理论最优解，代码错了马上能发现。第二步是梯度检验：在随机点用数值差分和解析梯度对比，误差一般要在1e-6以内，这能排除最隐蔽的梯度写错问题。第三步才是上标准测试函数，去比较不同算法的收敛行为。如果跳过前两步直接上Rosenbrock，一旦结果不对，你根本分不清是算法逻辑错了还是参数没调好。

还有一个习惯我觉得特别值得养成，就是每次跑实验之前记录当前代码版本对应的参数。有一段时间我同时在调Armijo的rho和罚因子倍数，结果跑了三次实验，三次结果都不一样，后来才发现rho从0.5改成了0.3，造成数据完全不可比。从那以后，我所有的实验脚本都会在输出结果时附带一份配置说明，把初始点、精度、最大迭代次数和线搜索参数全部打出来。这个习惯尤其在写第4章实验对比时特别有用。

5.2 这个作业还可以怎么扩展

如果时间充裕，我建议在完成基本要求之后，再做两个小扩展。第一个是给BFGS加上L-BFGS的简化版本，只用最近几轮的s和y向量来近似Hessian的逆。L-BFGS在高维问题上的优势非常明显，实现了它，你对"近似二阶方法"的理解会上升一个层次。第二个是画约束优化里可行域和目标函数等高线的叠加图。外点罚函数法的迭代路径图非常直观，能展示点如何从不可行区域逐步越过边界到达最优点。这两个扩展都能提升作业的完成度和报告的说服力，而且实现成本不算高，都是基于现有代码的增量修改。

上机作业这件事，表面上是写代码、跑实验、交报告，实际上是一次对优化算法从"听说"到"用熟"的转变。写代码的过程会逼你搞清楚每一个公式里的下标和符号，跑实验的过程会逼你去面对收敛性、稳定性和参数敏感度这些硬问题。这些体会是只看书完全不写代码根本得不到的。