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. 最速下降法实现与分析

最速下降法（梯度下降法）是最基础的一阶优化算法，其核心思想是在每一步沿着当前点梯度方向的负方向进行搜索。算法步骤如下：

1. 初始化起点 x₀ 和收敛阈值 ε
2. 计算当前梯度 ∇f(xₖ)
3. 如果 ‖∇f(xₖ)‖ < ε，停止迭代
4. 确定步长 αₖ（精确或非精确线搜索）
5. 更新 xₖ₊₁ = xₖ - αₖ∇f(xₖ)
6. 返回步骤 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)

牛顿法虽然收敛速度快（二阶收敛），但存在几个问题：

1. 需要计算和存储 Hessian 矩阵，计算量大
2. Hessian 矩阵可能不正定，导致搜索方向不是下降方向
3. 需要求解线性方程组，可能数值不稳定

阻尼牛顿法通过引入线搜索和改进 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 算法选择指南

根据问题特点选择合适算法：

5.4 实际应用建议

1. 问题规模：对于高维问题（n > 1000），梯度下降类方法更实用；低维问题可考虑牛顿法。

2. 信息可用性：

   • 只有函数值：使用无导数优化（如Nelder-Mead）
   • 有一阶导数：梯度下降、共轭梯度法
   • 有二阶导数：牛顿法、拟牛顿法

3. 约束处理：

   • 简单边界约束：投影梯度法
   • 线性约束：有效集方法
   • 非线性约束：罚函数法、增广拉格朗日法

4. 非凸问题：考虑多起点策略或全局优化方法（如模拟退火、遗传算法）

5. 实现技巧：

   • 使用自动微分计算精确梯度
   • 对病态问题考虑预处理
   • 监控收敛性（梯度范数、函数值变化、迭代点变化）

# 示例：使用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}")