最优化算法Python实战从梯度下降到KKT条件的代码实现期末考试临近最优化理论中的各种算法是否让你头疼不已与其反复翻看公式推导不如动手用Python将这些算法实现一遍。本文将带你用NumPy和SymPy复现梯度下降、牛顿法和KKT条件这三大最优化考题通过代码验证理论让抽象的概念变得直观可见。1. 环境准备与问题定义在开始编码之前我们需要准备好Python环境和问题定义。假设我们有以下三个考题需要解决最速下降法最小化函数 f(x₁, x₂) x₁² 2x₂² - 2x₁x₂ - 2x₂牛顿法最小化函数 f(x₁, x₂) x₁⁴ 2x₂⁴ x₁²x₂² - 3x₁x₂KKT条件带约束优化问题 min f(x₁, x₂) x₁² x₂²约束为 x₁ x₂ ≥ 1首先安装必要的库pip install numpy sympy matplotlib然后导入我们将用到的模块import numpy as np import sympy as sp from sympy import symbols, diff, solve, Eq import matplotlib.pyplot as plt from matplotlib import cm from mpl_toolkits.mplot3d import Axes3D2. 最速下降法的Python实现最速下降法梯度下降法是最基础的一阶优化算法其核心思想是沿着当前点的负梯度方向进行搜索。让我们先实现单步迭代的代码。2.1 定义目标函数和梯度def f1(x): 考题7的目标函数 return x[0]**2 2*x[1]**2 - 2*x[0]*x[1] - 2*x[1] def grad_f1(x): 考题7目标函数的梯度 return np.array([2*x[0] - 2*x[1], 4*x[1] - 2*x[0] - 2])2.2 最速下降法单步迭代def steepest_descent_step(x, grad_func, alpha0.1): 最速下降法单步迭代 :param x: 当前点 :param grad_func: 梯度函数 :param alpha: 固定步长 :return: 下一个迭代点 gradient grad_func(x) return x - alpha * gradient2.3 验证考题7的第一次迭代让我们按照考题要求从初始点x⁽⁰⁾ (0,0)开始进行第一次迭代x0 np.array([0., 0.]) gradient_at_x0 grad_f1(x0) x1 steepest_descent_step(x0, grad_f1) print(f初始点: {x0}) print(f初始梯度: {gradient_at_x0}) print(f第一次迭代后的点: {x1})输出结果应该与手算一致初始点: [0. 0.] 初始梯度: [ 0. -2.] 第一次迭代后的点: [0. 0.2]2.4 可视化迭代路径为了更直观理解算法的行为我们可以绘制迭代路径def plot_optimization_path(func, path, title): 绘制优化路径的等高线图 x np.linspace(-1, 1, 100) y np.linspace(-1, 1, 100) X, Y np.meshgrid(x, y) Z func(np.array([X, Y])) plt.figure(figsize(8, 6)) plt.contour(X, Y, Z, levels20, cmapviridis) plt.plot(path[:, 0], path[:, 1], ro-) plt.title(title) plt.xlabel(x1) plt.ylabel(x2) plt.colorbar() plt.show() # 进行多次迭代并绘制路径 path [x0] x x0.copy() for _ in range(10): x steepest_descent_step(x, grad_f1) path.append(x) path np.array(path) plot_optimization_path(f1, path, 最速下降法优化路径)3. 牛顿法的Python实现牛顿法是二阶优化算法利用了目标函数的Hessian矩阵信息通常比梯度下降法收敛更快。3.1 定义考题8的目标函数def f2(x): 考题8的目标函数 return x[0]**4 2*x[1]**4 x[0]**2*x[1]**2 - 3*x[0]*x[1] def grad_f2(x): 考题8目标函数的梯度 return np.array([4*x[0]**3 2*x[0]*x[1]**2 - 3*x[1], 8*x[1]**3 2*x[0]**2*x[1] - 3*x[0]]) def hessian_f2(x): 考题8目标函数的Hessian矩阵 return np.array([[12*x[0]**2 2*x[1]**2, 4*x[0]*x[1] - 3], [4*x[0]*x[1] - 3, 24*x[1]**2 2*x[0]**2]])3.2 牛顿法单步迭代def newton_step(x, grad_func, hessian_func): 牛顿法单步迭代 :param x: 当前点 :param grad_func: 梯度函数 :param hessian_func: Hessian矩阵函数 :return: 下一个迭代点 gradient grad_func(x) hessian hessian_func(x) delta_x np.linalg.solve(hessian, -gradient) return x delta_x3.3 验证考题8的第一次迭代从初始点x⁽⁰⁾ (1,1)开始x0_newton np.array([1., 1.]) x1_newton newton_step(x0_newton, grad_f2, hessian_f2) print(f初始点: {x0_newton}) print(f第一次迭代后的点: {x1_newton})3.4 比较牛顿法和最速下降法让我们比较两种方法在相同问题上的表现# 牛顿法迭代路径 path_newton [x0_newton] x x0_newton.copy() for _ in range(5): # 牛顿法通常收敛更快 x newton_step(x, grad_f2, hessian_f2) path_newton.append(x) path_newton np.array(path_newton) # 最速下降法迭代路径 path_sd [x0_newton] x x0_newton.copy() for _ in range(20): x steepest_descent_step(x, grad_f2, alpha0.01) # 需要更小的步长 path_sd.append(x) path_sd np.array(path_sd) # 绘制比较图 plt.figure(figsize(12, 6)) plt.subplot(1, 2, 1) plt.contour(X, Y, f2(np.array([X, Y])), levels20, cmapviridis) plt.plot(path_newton[:, 0], path_newton[:, 1], ro-) plt.title(牛顿法优化路径) plt.subplot(1, 2, 2) plt.contour(X, Y, f2(np.array([X, Y])), levels20, cmapviridis) plt.plot(path_sd[:, 0], path_sd[:, 1], bo-) plt.title(最速下降法优化路径) plt.show()4. KKT条件的自动推导与验证KKT条件是解决约束优化问题的关键。我们将使用SymPy来自动推导KKT条件并验证考题9的解。4.1 定义考题9的符号和函数# 定义符号变量 x1, x2, mu symbols(x1 x2 mu) # 定义目标函数和约束 f x1**2 x2**2 g x1 x2 - 1 # 约束改写为g(x) 0形式 # 构建拉格朗日函数 L f - mu * g4.2 自动推导KKT条件# 计算梯度条件 dL_dx1 diff(L, x1) dL_dx2 diff(L, x2) # 互补松弛条件 complementary mu * g # KKT条件方程组 kkt_conditions [ Eq(dL_dx1, 0), # 梯度条件1 Eq(dL_dx2, 0), # 梯度条件2 Eq(g, 0), # 原问题约束取等号因为mu0时g0 mu 0 # 乘子非负 ] # 求解KKT条件 solution solve(kkt_conditions, (x1, x2, mu), dictTrue) print(KKT条件的解:) for sol in solution: print(sol)4.3 验证解的最优性# 验证解是否满足二阶充分条件 def check_optimality(x_star, mu_star): 验证二阶充分条件 # 计算Hessian矩阵 H np.array([[2, 0], [0, 2]]) # 计算约束的梯度 grad_g np.array([1, 1]) # 检查二阶条件 # 对于所有满足grad_g.dot(d) 0的dd.T H d 0 # 因为H是正定矩阵这个条件自动满足 return True x_star np.array([float(solution[0][x1]), float(solution[0][x2])]) mu_star float(solution[0][mu]) if check_optimality(x_star, mu_star): print(f点 {x_star} 满足二阶充分条件是最优解) else: print(f点 {x_star} 不满足二阶充分条件)4.4 可视化约束优化问题# 绘制约束优化问题的可行域和最优解 x np.linspace(-1, 2, 100) y np.linspace(-1, 2, 100) X, Y np.meshgrid(x, y) Z X**2 Y**2 plt.figure(figsize(8, 6)) plt.contour(X, Y, Z, levels20, cmapviridis) plt.plot(x, 1 - x, r-, label约束边界 x1 x2 1) plt.fill_between(x, 1 - x, 2, colorred, alpha0.1, label可行域) plt.plot(x_star[0], x_star[1], ro, label最优解) plt.title(约束优化问题可视化) plt.xlabel(x1) plt.ylabel(x2) plt.legend() plt.colorbar() plt.show()5. 算法扩展与实用技巧在实际应用中我们还需要考虑一些实现细节和优化技巧。5.1 自适应步长的梯度下降固定步长可能导致收敛问题我们可以实现Armijo线搜索def armijo_line_search(x, grad_func, beta0.5, sigma0.1): Armijo线搜索确定步长 gradient grad_func(x) alpha 1.0 while f1(x - alpha * gradient) f1(x) - sigma * alpha * np.dot(gradient, gradient): alpha * beta return alpha def adaptive_steepest_descent(x0, grad_func, max_iter100, tol1e-6): 带自适应步长的最速下降法 x x0.copy() path [x0] for _ in range(max_iter): grad grad_func(x) if np.linalg.norm(grad) tol: break alpha armijo_line_search(x, grad_func) x x - alpha * grad path.append(x) return np.array(path) # 测试自适应步长算法 path_adaptive adaptive_steepest_descent(np.array([0., 0.]), grad_f1) plot_optimization_path(f1, path_adaptive, 自适应步长最速下降法)5.2 处理牛顿法的Hessian矩阵不正定问题当Hessian矩阵不正定时牛顿方向可能不是下降方向。我们可以通过修改Hessian矩阵来解决def modified_newton_step(x, grad_func, hessian_func, epsilon1e-6): 修正的牛顿法处理Hessian不正定情况 gradient grad_func(x) hessian hessian_func(x) # 检查Hessian是否正定 min_eigval np.min(np.linalg.eigvals(hessian)) if min_eigval 0: # 不正定添加正则化项 hessian (epsilon - min_eigval) * np.eye(hessian.shape[0]) delta_x np.linalg.solve(hessian, -gradient) return x delta_x5.3 约束优化问题的数值解法对于更复杂的约束问题我们可以实现惩罚函数法def penalty_method(f, g, x0, mu1.0, rho2.0, max_iter100, tol1e-6): 外点惩罚函数法 x x0.copy() path [x0] for _ in range(max_iter): # 定义惩罚函数 def P(x): penalty mu * min(0, g(x))**2 return f(x) penalty # 使用牛顿法最小化惩罚函数 # 这里简化实现实际应用中需要更鲁棒的优化器 grad_P ... # 惩罚函数的梯度 hessian_P ... # 惩罚函数的Hessian x modified_newton_step(x, grad_P, hessian_P) path.append(x) # 检查约束违反程度 if g(x) -tol: break # 增加惩罚系数 mu * rho return np.array(path)6. 性能分析与比较为了更深入地理解这些算法的行为让我们对它们的性能进行系统分析。6.1 收敛速度比较我们可以比较不同算法在相同问题上的收敛速度def track_convergence(algorithm, x0, *args, max_iter50, tol1e-6, **kwargs): 跟踪算法的收敛过程 x x0.copy() f_values [f1(x)] grad_norms [np.linalg.norm(grad_f1(x))] for _ in range(max_iter): x algorithm(x, *args, **kwargs) f_values.append(f1(x)) grad_norm np.linalg.norm(grad_f1(x)) grad_norms.append(grad_norm) if grad_norm tol: break return f_values, grad_norms # 测试不同算法 sd_f, sd_gn track_convergence(steepest_descent_step, np.array([0., 0.]), grad_f1, alpha0.1) newton_f, newton_gn track_convergence(newton_step, np.array([0., 0.]), grad_f1, lambda x: np.array([[2, -2], [-2, 4]])) adaptive_f, adaptive_gn track_convergence(adaptive_steepest_descent, np.array([0., 0.]), grad_f1) # 绘制收敛曲线 plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.semilogy(sd_f, b-, label最速下降法) plt.semilogy(newton_f, r-, label牛顿法) plt.semilogy(adaptive_f, g-, label自适应步长) plt.title(函数值收敛曲线) plt.xlabel(迭代次数) plt.ylabel(f(x)) plt.legend() plt.subplot(1, 2, 2) plt.semilogy(sd_gn, b-, label最速下降法) plt.semilogy(newton_gn, r-, label牛顿法) plt.semilogy(adaptive_gn, g-, label自适应步长) plt.title(梯度范数收敛曲线) plt.xlabel(迭代次数) plt.ylabel(||∇f(x)||) plt.legend() plt.show()6.2 不同初始点的行为分析算法的表现可能依赖于初始点的选择。让我们分析不同初始点下的行为initial_points [np.array([0., 0.]), np.array([1., -1.]), np.array([-1., 1.])] plt.figure(figsize(15, 5)) for i, x0 in enumerate(initial_points, 1): path adaptive_steepest_descent(x0, grad_f1) plt.subplot(1, 3, i) x np.linspace(-1.5, 1.5, 100) y np.linspace(-1.5, 1.5, 100) X, Y np.meshgrid(x, y) Z f1(np.array([X, Y])) plt.contour(X, Y, Z, levels20, cmapviridis) plt.plot(path[:, 0], path[:, 1], ro-) plt.title(f初始点: {x0}) plt.xlabel(x1) plt.ylabel(x2) plt.tight_layout() plt.show()6.3 算法复杂度分析让我们从理论和实践两个角度分析算法的复杂度算法每次迭代的计算复杂度理论收敛速度实际观察到的收敛速度最速下降法O(n)线性线性有时较慢牛顿法O(n³)因矩阵求逆二次接近二次当接近解时修正牛顿法O(n³)超线性取决于修正策略注意在实际实现中牛顿法的Hessian矩阵求逆可以通过Cholesky分解等方法来优化特别是当矩阵稀疏时。7. 实际应用案例为了加深理解让我们将这些算法应用于一些实际问题中。7.1 线性回归问题我们可以将最优化算法应用于简单的线性回归问题# 生成线性回归数据 np.random.seed(42) X np.random.rand(100, 1) y 3 * X 2 np.random.randn(100, 1) * 0.1 # 定义损失函数和梯度 def linear_regression_loss(theta): 线性回归的损失函数 return np.sum((X.dot(theta[1:2]) theta[0] - y)**2) def linear_regression_grad(theta): 线性回归的梯度 error X.dot(theta[1:2]) theta[0] - y grad_w 2 * X.T.dot(error) grad_b 2 * np.sum(error) return np.array([grad_b, grad_w[0]]) # 使用自适应最速下降法求解 theta0 np.array([0., 0.]) path_lr adaptive_steepest_descent(theta0, linear_regression_grad) # 绘制结果 plt.figure(figsize(12, 5)) plt.subplot(1, 2, 1) plt.scatter(X, y) x_plot np.array([[0], [1]]) for i, theta in enumerate(path_lr[::5]): y_plot x_plot * theta[1] theta[0] plt.plot(x_plot, y_plot, alpha0.3, labelfiter {i*5} if i 3 else None) plt.title(线性回归拟合过程) plt.legend() plt.subplot(1, 2, 2) losses [linear_regression_loss(theta) for theta in path_lr] plt.semilogy(losses) plt.title(损失函数收敛曲线) plt.xlabel(迭代次数) plt.ylabel(损失值) plt.show()7.2 逻辑回归问题我们还可以将这些算法应用于逻辑回归问题# 生成二分类数据 np.random.seed(42) X np.random.randn(100, 2) y (X[:, 0] X[:, 1] 0).astype(int) # 定义sigmoid函数 def sigmoid(z): return 1 / (1 np.exp(-z)) # 定义逻辑回归的损失函数和梯度 def logistic_loss(theta): z X.dot(theta[1:]) theta[0] return -np.sum(y * np.log(sigmoid(z)) (1 - y) * np.log(1 - sigmoid(z))) def logistic_grad(theta): z X.dot(theta[1:]) theta[0] error sigmoid(z) - y grad_w X.T.dot(error) grad_b np.sum(error) return np.array([grad_b, grad_w[0], grad_w[1]]) # 使用牛顿法求解 def logistic_hessian(theta): z X.dot(theta[1:]) theta[0] s sigmoid(z) * (1 - sigmoid(z)) H np.zeros((3, 3)) X_aug np.hstack([np.ones((100, 1)), X]) H X_aug.T.dot(s.reshape(-1, 1) * X_aug) return H theta0 np.array([0., 0., 0.]) path_logistic [] x theta0.copy() for _ in range(10): path_logistic.append(x) x newton_step(x, logistic_grad, logistic_hessian) path_logistic np.array(path_logistic) # 绘制决策边界 plt.figure(figsize(8, 6)) plt.scatter(X[:, 0], X[:, 1], cy, cmapbwr) x1_min, x1_max X[:, 0].min(), X[:, 0].max() x2_min, x2_max X[:, 1].min(), X[:, 1].max() xx1, xx2 np.meshgrid(np.linspace(x1_min, x1_max, 100), np.linspace(x2_min, x2_max, 100)) for i, theta in enumerate(path_logistic): if i % 2 0: Z sigmoid(theta[0] xx1 * theta[1] xx2 * theta[2]) 0.5 plt.contour(xx1, xx2, Z, levels[0.5], alpha0.3, colors[black]) plt.title(逻辑回归决策边界演变) plt.show()7.3 带约束的投资组合优化最后我们看一个带约束的优化问题示例# 投资组合优化问题 # 最小化风险min w^T Σ w # 约束w^T μ target_return, sum(w) 1, w 0 # 生成随机资产数据 np.random.seed(42) n_assets 5 mu np.random.rand(n_assets) * 0.1 # 预期收益 Sigma np.random.randn(n_assets, n_assets) Sigma Sigma.T Sigma # 协方差矩阵 target_return 0.05 # 使用KKT条件求解 w symbols([fw{i} for i in range(n_assets)] [lambda1, lambda2] [fmu{i} for i in range(n_assets)]) # 构建拉格朗日函数 risk sum(Sigma[i,j] * w[i] * w[j] for i in range(n_assets) for j in range(n_assets)) return_constraint sum(mu[i] * w[i] for i in range(n_assets)) - target_return budget_constraint sum(w[i] for i in range(n_assets)) - 1 nonneg_constraints [w[i] for i in range(n_assets)] L risk - w[n_assets] * return_constraint - w[n_assets1] * budget_constraint sum(w[n_assets2i] * nonneg_constraints[i] for i in range(n_assets)) # 构建KKT条件 kkt_conditions [] for i in range(n_assets): kkt_conditions.append(Eq(diff(L, w[i]), 0)) kkt_conditions.append(Eq(return_constraint, 0)) kkt_conditions.append(Eq(budget_constraint, 0)) for i in range(n_assets): kkt_conditions.append(Eq(w[n_assets2i] * w[i], 0)) kkt_conditions.append(w[n_assets2i] 0) kkt_conditions.append(w[i] 0) # 求解KKT条件实际应用中可能使用数值方法 # solution solve(kkt_conditions, w) # 对于大问题符号求解可能不可行在实际应用中对于这种规模的优化问题我们通常会使用数值优化库如scipy.optimize或专门的凸优化库如cvxpy来求解。