Python 实现最优化 5 大经典算法：梯度下降、牛顿法、罚函数法实战与收敛性对比

Python 实现最优化 5 大经典算法梯度下降、牛顿法、罚函数法实战与收敛性对比在机器学习和科学计算领域最优化算法扮演着至关重要的角色。无论是训练神经网络还是求解复杂的物理方程找到目标函数的最小值点都是核心任务。本文将深入探讨五种经典的最优化算法实现并通过 Python 代码展示它们在实际问题中的应用表现。1. 最优化算法基础与环境准备在开始实现算法之前我们需要明确一些基本概念并准备好开发环境。最优化问题通常可以表述为寻找使目标函数 f(x) 最小化的变量 x 的值。根据问题的不同可能还包含等式或不等式约束条件。对于本次实验我们将使用以下 Python 库import numpy as np import matplotlib.pyplot as plt from scipy.optimize import line_search from mpl_toolkits.mplot3d import Axes3D为了评估算法性能我们选择 Rosenbrock 函数作为测试基准。这个被称为香蕉函数的测试函数因其弯曲的谷状结构而闻名是检验优化算法性能的经典选择def rosenbrock(x): Rosenbrock函数实现 return 100 * (x[1] - x[0]**2)**2 (1 - x[0])**2 def rosenbrock_grad(x): Rosenbrock函数的梯度 return np.array([ -400 * x[0] * (x[1] - x[0]**2) - 2 * (1 - x[0]), 200 * (x[1] - x[0]**2) ]) def rosenbrock_hess(x): Rosenbrock函数的Hessian矩阵 return np.array([ [1200 * x[0]**2 - 400 * x[1] 2, -400 * x[0]], [-400 * x[0], 200] ])为了直观理解这个函数我们可以绘制其在[-2, 2]×[-1, 3]区域内的三维图像x np.linspace(-2, 2, 100) y np.linspace(-1, 3, 100) X, Y np.meshgrid(x, y) Z rosenbrock([X, Y]) fig plt.figure(figsize(12, 8)) ax fig.add_subplot(111, projection3d) ax.plot_surface(X, Y, Z, cmapviridis, alpha0.8) ax.set_xlabel(x) ax.set_ylabel(y) ax.set_zlabel(f(x,y)) plt.title(Rosenbrock Function) plt.show()Rosenbrock 函数在 (1,1) 处有全局最小值 0但到达这个点需要沿着一条狭窄弯曲的山谷前进这使得许多优化算法收敛缓慢。2. 最速下降法实现与分析最速下降法梯度下降法是最基础的一阶优化算法其核心思想是在每一步沿着当前点梯度方向的负方向进行搜索。算法步骤如下初始化起点 x₀ 和收敛阈值 ε计算当前梯度 ∇f(xₖ)如果 ‖∇f(xₖ)‖ ε停止迭代确定步长 αₖ精确或非精确线搜索更新 xₖ₊₁ xₖ - αₖ∇f(xₖ)返回步骤 2以下是 Python 实现代码def gradient_descent(f, grad, x0, max_iter10000, tol1e-6, alpha_methodexact): 最速下降法实现 参数: f: 目标函数 grad: 梯度函数 x0: 初始点 max_iter: 最大迭代次数 tol: 收敛阈值 alpha_method: 步长选择方法 (exact或armijo) x x0.copy() trajectory [x0.copy()] grad_norms [] for k in range(max_iter): g grad(x) g_norm np.linalg.norm(g) grad_norms.append(g_norm) if g_norm tol: break # 步长选择 if alpha_method exact: alpha line_search(f, grad, x, -g)[0] if alpha is None: alpha 0.001 # 备用步长 elif alpha_method armijo: alpha 1.0 beta 0.5 sigma 0.1 while f(x - alpha * g) f(x) - sigma * alpha * g_norm**2: alpha * beta # 更新点 x x - alpha * g trajectory.append(x.copy()) return x, np.array(trajectory), np.array(grad_norms)我们可以测试这个方法在 Rosenbrock 函数上的表现x0 np.array([-1.5, 2.5]) x_opt, trajectory, grad_norms gradient_descent(rosenbrock, rosenbrock_grad, x0) plt.figure(figsize(12, 6)) plt.semilogy(grad_norms) plt.xlabel(Iteration) plt.ylabel(Gradient norm) plt.title(Convergence of Gradient Descent) plt.grid() plt.show()最速下降法的主要优缺点优点实现简单计算量小仅需一阶导数内存消耗低适合大规模问题保证收敛到局部极小点在适当条件下缺点收敛速度线性在病态条件下非常缓慢在 Rosenbrock 函数等具有峡谷形状的函数上表现不佳步长选择对性能影响很大提示在实际应用中精确线搜索往往计算代价太高通常采用 Wolfe 条件或 Armijo 规则等非精确线搜索方法。3. 牛顿法及其改进实现牛顿法利用二阶导数信息构建二次模型能够提供更快的收敛速度。基本牛顿法的迭代公式为xₖ₊₁ xₖ - [∇²f(xₖ)]⁻¹∇f(xₖ)以下是基本牛顿法的 Python 实现def newton_method(f, grad, hess, x0, max_iter100, tol1e-6): 基本牛顿法实现 参数: f: 目标函数 grad: 梯度函数 hess: Hessian矩阵函数 x0: 初始点 max_iter: 最大迭代次数 tol: 收敛阈值 x x0.copy() trajectory [x0.copy()] grad_norms [] for k in range(max_iter): g grad(x) H hess(x) g_norm np.linalg.norm(g) grad_norms.append(g_norm) if g_norm tol: break try: d np.linalg.solve(H, -g) except np.linalg.LinAlgError: # Hessian矩阵奇异使用梯度下降方向 d -g # 线搜索确定步长 alpha line_search(f, grad, x, d)[0] if alpha is None: alpha 1.0 x x alpha * d trajectory.append(x.copy()) return x, np.array(trajectory), np.array(grad_norms)牛顿法虽然收敛速度快二阶收敛但存在几个问题需要计算和存储 Hessian 矩阵计算量大Hessian 矩阵可能不正定导致搜索方向不是下降方向需要求解线性方程组可能数值不稳定阻尼牛顿法通过引入线搜索和改进 Hessian 处理来解决这些问题def damped_newton(f, grad, hess, x0, max_iter100, tol1e-6, epsilon1e-6): 阻尼牛顿法实现 参数: epsilon: 用于保证Hessian矩阵正定性的小常数 x x0.copy() trajectory [x0.copy()] grad_norms [] for k in range(max_iter): g grad(x) H hess(x) g_norm np.linalg.norm(g) grad_norms.append(g_norm) if g_norm tol: break # 修正Hessian矩阵确保正定性 min_eigval np.min(np.linalg.eigvals(H)) if min_eigval epsilon: H H (epsilon - min_eigval) * np.eye(len(x)) try: d np.linalg.solve(H, -g) except np.linalg.LinAlgError: d -g # Armijo线搜索 alpha 1.0 beta 0.5 sigma 0.1 while f(x alpha * d) f(x) sigma * alpha * np.dot(g, d): alpha * beta x x alpha * d trajectory.append(x.copy()) return x, np.array(trajectory), np.array(grad_norms)我们可以比较基本牛顿法和阻尼牛顿法的性能x0 np.array([-1.5, 2.5]) # 基本牛顿法 x_opt1, traj1, gnorms1 newton_method(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0) # 阻尼牛顿法 x_opt2, traj2, gnorms2 damped_newton(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0) # 绘制收敛曲线 plt.figure(figsize(12, 6)) plt.semilogy(gnorms1, labelBasic Newton) plt.semilogy(gnorms2, labelDamped Newton) plt.xlabel(Iteration) plt.ylabel(Gradient norm) plt.title(Comparison of Newton Methods) plt.legend() plt.grid() plt.show()牛顿类算法的特点优点二阶收敛速度在接近解时收敛极快对于二次函数一步即可收敛到最优解能够处理曲率信息在病态条件下表现良好缺点需要计算和存储 Hessian 矩阵计算复杂度高每次迭代需要求解线性方程组可能数值不稳定远离最优解时 Hessian 可能不正定需要修正4. 罚函数法处理约束优化问题前面讨论的方法适用于无约束优化问题。对于约束优化问题罚函数法是一种常用的处理方式。我们重点介绍外点罚函数法和内点罚函数法。4.1 外点罚函数法外点罚函数法通过将约束违反量加入目标函数将约束问题转化为一系列无约束问题。考虑问题min f(x) s.t. g_i(x) ≥ 0, i ∈ I h_j(x) 0, j ∈ E外点罚函数定义为P(x, σ) f(x) σ∑[min(0, g_i(x))]² σ∑h_j(x)²Python 实现def exterior_penalty(f, ineq_constraints, eq_constraints, x0, sigma01.0, sigma_max1e6, max_iter100, tol1e-6): 外点罚函数法实现 参数: ineq_constraints: 不等式约束列表[g1, g2,...] eq_constraints: 等式约束列表[h1, h2,...] sigma0: 初始罚参数 sigma_max: 最大罚参数 x x0.copy() sigma sigma0 trajectory [x0.copy()] constraint_violations [] def penalty_function(x, sigma): penalty 0.0 for g in ineq_constraints: penalty min(0, g(x))**2 for h in eq_constraints: penalty h(x)**2 return f(x) sigma * penalty def penalty_gradient(x, sigma): grad grad_f(x) for g in ineq_constraints: if g(x) 0: grad 2 * sigma * min(0, g(x)) * grad_g(x) for h in eq_constraints: grad 2 * sigma * h(x) * grad_h(x) return grad for k in range(max_iter): # 求解无约束子问题 res minimize(penalty_function, x, args(sigma,), jacpenalty_gradient, methodBFGS) x res.x trajectory.append(x.copy()) # 计算约束违反量 violation 0.0 for g in ineq_constraints: violation abs(min(0, g(x))) for h in eq_constraints: violation abs(h(x)) constraint_violations.append(violation) if violation tol: break # 增加罚参数 sigma min(10 * sigma, sigma_max) return x, np.array(trajectory), np.array(constraint_violations)4.2 内点罚函数法内点罚函数法障碍函数法要求初始点在可行域内部并通过障碍函数阻止迭代点接近边界。对于不等式约束问题min f(x) s.t. g_i(x) ≥ 0, i ∈ I内点罚函数定义为B(x, μ) f(x) - μ∑ln(g_i(x))Python 实现def interior_point(f, ineq_constraints, x0, mu01.0, mu_min1e-6, max_iter100, tol1e-6): 内点罚函数法实现 参数: mu0: 初始障碍参数 mu_min: 最小障碍参数 x x0.copy() mu mu0 trajectory [x0.copy()] constraint_violations [] def barrier_function(x, mu): barrier 0.0 for g in ineq_constraints: barrier - np.log(g(x)) return f(x) mu * barrier def barrier_gradient(x, mu): grad grad_f(x) for g in ineq_constraints: grad - mu * grad_g(x) / g(x) return grad for k in range(max_iter): # 求解无约束子问题 res minimize(barrier_function, x, args(mu,), jacbarrier_gradient, methodBFGS) x res.x trajectory.append(x.copy()) # 计算约束违反量 violation 0.0 for g in ineq_constraints: violation abs(min(0, g(x))) constraint_violations.append(violation) if mu mu_min and violation tol: break # 减小障碍参数 mu max(mu / 10, mu_min) return x, np.array(trajectory), np.array(constraint_violations)4.3 罚函数法应用示例考虑以下约束优化问题min (x₁ - 2)² (x₂ - 1)² s.t. x₁ x₂ ≤ 2 x₁² - x₂ ≤ 0我们可以定义约束和目标函数def f(x): return (x[0] - 2)**2 (x[1] - 1)**2 def g1(x): return 2 - x[0] - x[1] # 转换为g(x)≥0形式 def g2(x): return x[1] - x[0]**2 ineq_constraints [g1, g2] eq_constraints [] # 外点罚函数法求解 x0 np.array([0.5, 0.5]) x_opt_ext, traj_ext, viol_ext exterior_penalty(f, ineq_constraints, eq_constraints, x0) # 内点罚函数法求解 x0_feasible np.array([0.5, 0.8]) # 必须在可行域内部 x_opt_int, traj_int, viol_int interior_point(f, ineq_constraints, x0_feasible) # 绘制结果 plt.figure(figsize(12, 6)) plt.semilogy(viol_ext, labelExterior Penalty) plt.semilogy(viol_int, labelInterior Point) plt.xlabel(Iteration) plt.ylabel(Constraint Violation) plt.title(Constraint Violation Reduction) plt.legend() plt.grid() plt.show()罚函数法的优缺点比较方法优点缺点外点罚函数法可从不可行点开始实现简单解序列仅在极限处可行需要σ→∞内点罚函数法解序列始终可行数值稳定需要可行初始点处理等式约束困难5. 算法性能对比与实战建议现在我们将五种算法在 Rosenbrock 函数上的表现进行系统对比。除了收敛速度我们还将考虑计算复杂度和实际应用中的注意事项。5.1 收敛速度对比我们统一从点 (-1.5, 2.5) 出发比较各算法的梯度范数下降情况x0 np.array([-1.5, 2.5]) # 运行各算法 _, traj_gd, gnorms_gd gradient_descent(rosenbrock, rosenbrock_grad, x0) _, traj_newton, gnorms_newton newton_method(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0) _, traj_damped, gnorms_damped damped_newton(rosenbrock, rosenbrock_grad, rosenbrock_hess, x0) # 绘制收敛曲线 plt.figure(figsize(12, 6)) plt.semilogy(gnorms_gd, labelGradient Descent) plt.semilogy(gnorms_newton, labelNewton Method) plt.semilogy(gnorms_damped, labelDamped Newton) plt.xlabel(Iteration) plt.ylabel(Gradient Norm (log scale)) plt.title(Convergence Comparison on Rosenbrock Function) plt.legend() plt.grid() plt.show()5.2 迭代路径可视化通过绘制各算法在二维平面上的搜索路径可以直观了解它们的搜索行为# 绘制等高线图 x np.linspace(-2, 2, 100) y np.linspace(-1, 3, 100) X, Y np.meshgrid(x, y) Z rosenbrock([X, Y]) plt.figure(figsize(12, 10)) plt.contour(X, Y, Z, levelsnp.logspace(0, 3, 20), cmapviridis) plt.plot(1, 1, r*, markersize15, labelOptimum) # 绘制搜索路径 plt.plot(traj_gd[:, 0], traj_gd[:, 1], o-, labelGradient Descent) plt.plot(traj_newton[:, 0], traj_newton[:, 1], s-, labelNewton Method) plt.plot(traj_damped[:, 0], traj_damped[:, 1], d-, labelDamped Newton) plt.xlabel(x) plt.ylabel(y) plt.title(Optimization Paths on Rosenbrock Function) plt.legend() plt.colorbar() plt.show()5.3 算法选择指南根据问题特点选择合适算法算法类型适用场景注意事项梯度下降法大规模问题一阶信息易得内存有限学习率选择关键可能收敛慢牛顿法中小规模问题二阶信息可用需要快速收敛计算Hessian代价高可能数值不稳定阻尼牛顿法Hessian不正定时需要更鲁棒收敛比标准牛顿法更稳定仍需要二阶信息外点罚函数法约束优化初始点不可行最终解可能轻微违反约束罚参数选择重要内点罚函数法约束优化可行初始点可得保持可行性处理等式约束需特殊技巧5.4 实际应用建议问题规模对于高维问题n 1000梯度下降类方法更实用低维问题可考虑牛顿法。信息可用性只有函数值使用无导数优化如Nelder-Mead有一阶导数梯度下降、共轭梯度法有二阶导数牛顿法、拟牛顿法约束处理简单边界约束投影梯度法线性约束有效集方法非线性约束罚函数法、增广拉格朗日法非凸问题考虑多起点策略或全局优化方法如模拟退火、遗传算法实现技巧使用自动微分计算精确梯度对病态问题考虑预处理监控收敛性梯度范数、函数值变化、迭代点变化# 示例使用scipy的优化器进行比较 from scipy.optimize import minimize x0 np.array([-1.5, 2.5]) methods [CG, BFGS, Newton-CG, trust-krylov] labels [Conjugate Gradient, BFGS, Newton-CG, Trust-Region] results [] for method in methods: res minimize(rosenbrock, x0, methodmethod, jacrosenbrock_grad, hessrosenbrock_hess if method in [Newton-CG, trust-krylov] else None) results.append(res) # 比较迭代次数和函数调用次数 print(Method\t\tIterations\tFunction Evals) for i, res in enumerate(results): print(f{labels[i]:15}{res.nit:10}\t{res.nfev:12})