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, projection='3d') ax.plot_surface(X, Y, Z, cmap='viridis', alpha=0.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_iter=10000, tol=1e-6, alpha_method='exact'): """ 最速下降法实现 参数: 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_iter=100, tol=1e-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_iter=100, tol=1e-6, epsilon=1e-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, label='Basic Newton') plt.semilogy(gnorms2, label='Damped 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, sigma0=1.0, sigma_max=1e6, max_iter=100, tol=1e-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,), jac=penalty_gradient, method='BFGS') 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, mu0=1.0, mu_min=1e-6, max_iter=100, tol=1e-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,), jac=barrier_gradient, method='BFGS') 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, label='Exterior Penalty') plt.semilogy(viol_int, label='Interior 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, label='Gradient Descent') plt.semilogy(gnorms_newton, label='Newton Method') plt.semilogy(gnorms_damped, label='Damped 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, levels=np.logspace(0, 3, 20), cmap='viridis') plt.plot(1, 1, 'r*', markersize=15, label='Optimum') # 绘制搜索路径 plt.plot(traj_gd[:, 0], traj_gd[:, 1], 'o-', label='Gradient Descent') plt.plot(traj_newton[:, 0], traj_newton[:, 1], 's-', label='Newton Method') plt.plot(traj_damped[:, 0], traj_damped[:, 1], 'd-', label='Damped 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, method=method, jac=rosenbrock_grad, hess=rosenbrock_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}")