拟牛顿法：python代码实现 |
您所在的位置：网站首页 › python矩阵运算代码 › 拟牛顿法：python代码实现 |
拟牛顿法：python代码实现
|
拟牛顿法 在梯度类算法原理：最速下降法、牛顿法和拟牛顿法中，介绍了梯度类算法求解优化问题的设计思路，并以最速下降法、牛顿法和拟牛顿法为例，描述了具体的算法迭代过程。 其中，拟牛顿法（Broyden–Fletcher–Goldfarb–Shanno，BFGS）在实际优化场景中被广泛使用，因此本文将自主编写Python代码，实现BFGS的全部过程，并且和现有工具包做对比，从而加深对BFGS的理解。 待优化实例本文使用的待优化实例为10维的Rosenbrock函数： f(x)=\sum_{i=2}^{10} [100(x_i-x_{i-1}^2)^2+(1-x_{i-1})^2] \\ 10维图像难以可视化，此处编写代码绘制2维函数的等高线图，更直观地认识该函数。 import numpy as np import matplotlib.pyplot as plt # 任意维度Rosenbrock函数 def rosen(x): return sum(100 * (x[1:]-x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2) # 绘制等高线图 def plot_contour(f): # x和y区间均为[0, 2] x = np.linspace(0, 2, 100) y = np.linspace(0, 2, 100) # 将原始数据变成网格数据形式 X, Y = np.meshgrid(x, y) # 计算对应目标函数值 Z = np.zeros_like(X) for i in range(Z.shape[0]): for j in range(Z.shape[1]): Z[i, j] = f(np.array([X[i, j], Y[i, j]])) # 设置打开画布大小,长10，宽6 plt.figure(figsize=(10, 6)) # 画等高线 cset = plt.contourf(X, Y, Z, 6, cmap=plt.cm.hot) # 颜色分层，6层 contour = plt.contour(X, Y, Z, 8, colors='k') # 绘制8条等高线 plt.clabel(contour, fontsize=10, colors='k') # 等高线上标明Z值 plt.colorbar(cset) # 显示颜色条 plt.scatter(1, 1, marker='p') # 标明最优解位置 plt.title('2D-Rosenbrock contour') # 增加标题 plt.show() if __name__ == '__main__': plot_contour(rosen)以下为2维Rosenbrock的等高线图。 其大致呈抛物线形，最小值位于抛物线形的山谷中（香蕉型山谷，图中暗红色区域，最小值位置[1,1]）。 这个山谷通常很容易被找到，但由于山谷内的值变化不大，要找到最小值比较困难。  scipy工具包实现BFGSpython的scipy.optimize工具包能够实现BFGS，具体可以参考scipy.optimize优化器的各种使用。 本节仅做简要描述。 主要的求解方案有两种：第一种是只将待优化函数和初值传递给工具包，其余内容全部由工具包自行计算；第二种是将待优化函数的雅可比矩阵也传递给工具包，提升求解的速度。 雅可比矩阵之前没提到过，这里解释一下。 雅可比矩阵是一阶偏导数以一定方式排列成的矩阵。 假设\pmb f=[f_1,f_2,···,f_m]，\pmb x=[x_1,x_2,···,x_n]，那么\pmb f的雅可比矩阵 \pmb J 为 m×n 的矩阵: \pmb J=[\frac{\partial \pmb f}{\partial x_1}···\frac{\partial \pmb f}{\partial x_n}]=\left[ \begin{matrix} \frac{\partial f_1}{\partial x_1} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \vdots & \ddots & \vdots \\ \frac{\partial f_m}{\partial x_1} & \cdots & \frac{\partial f_m}{\partial x_n} \\ \end{matrix} \right]\\ 针对Rosenbrock函数，一阶偏导数为 \begin{eqnarray} \frac{\partial f}{\partial x_j} & =& \sum_{i=2}^{10} [100\frac{\partial (x_i-x_{i-1}^2)^2}{\partial x_j}+\frac{\partial ((1-x_{i-1})^2}{\partial x_j}]\\ & = & \sum_{i=2}^{10} [200(x_i-x_{i-1}^2)(\frac{\partial x_i}{\partial x_j}-2x_{i-1}\frac{\partial x_{i-1}}{\partial x_j})-2(1-x_{i-1})\frac{\partial x_{i-1}}{\partial x_j}] \\ & = & 200(x_j-x_{j-1}^2)-400x_j(x_{j+1}-x_j^2)-2(1-x_j) \end{eqnarray}\\ 当j=1时 \frac{\partial f}{\partial x_1}=-400x_1(x_2-x_1^2)-2(1-x_1) \\ 当j=10时 \frac{\partial f}{\partial x_{10}}=200(x_{10}-x_9^2) \\ 两种求解方案的代码如下： import numpy as np from scipy.optimize import minimize import time def rosen(x): return sum(100 * (x[1:]-x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2) def rosen_der(x): xm = x[1:-1] xm_m1 = x[:-2] xm_p1 = x[2:] der = np.zeros_like(x) der[1:-1] = 200 * (xm - xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1 - xm) der[0] = -400 * x[0] * (x[1] -x[0] ** 2) - 2 * (1 - x[0]) der[-1] = 200 * (x[-1] - x[-2] ** 2) return der if __name__ == '__main__': np.random.seed(0) # 设置种子，保证结果可重复 x0 = np.random.rand(10) # 随机生成初始值 print('------------不提供jac计算--------------') res = minimize(rosen, x0, options={'disp': True}) # 不提供jac print('-------------提供jac计算---------------') res1 = minimize(rosen, x0, jac=rosen_jac, options={'disp': True}) # 提供jac # 不提供jac，统计计算时间 start_time = time.time() for i in range(200): res = minimize(rosen, x0, options={'disp': False}) end_time = time.time() calc_time_non_jac = end_time - start_time # 提供jac，统计计算时间 start_time = time.time() for i in range(200): res = minimize(rosen, x0, jac=rosen_jac, options={'disp': False}) end_time = time.time() time1 = end_time - start_time calc_time_with_jac = end_time - start_time # 对比计算时间 print('-------------评估jac效率提升---------------') print('不提供jac时，计算时间为：{} 秒; 提供jac后，计算时间为：{} 秒，时间降低为原来的：{}'.format( calc_time_non_jac, calc_time_with_jac, calc_time_with_jac / calc_time_non_jac))输出结果如下。 两个方案都可以找到最优解；当不提供雅可比矩阵时，目标函数计算了572次，提供雅可比矩阵后，目标函数的计算次数降至52次；具体到计算效率上，提供雅可比矩阵可以将计算时间降至原先的30%。 ------------不提供jac计算-------------- Optimization terminated successfully. Current function value: 0.000000 Iterations: 41 Function evaluations: 572 Gradient evaluations: 52 -------------提供jac计算--------------- Optimization terminated successfully. Current function value: 0.000000 Iterations: 42 Function evaluations: 52 Gradient evaluations: 52 -------------评估jac效率提升--------------- 不提供jac时，计算时间为：3.727977991104126 秒; 提供jac后，计算时间为：1.1221048831939697 秒，时间降低为原来的：0.3009955761197058 自编Python实现BFGS自行编写代码时，主要关注三项内容： (1) 迭代公式：x_{k+1} = x_k + \alpha * G * g。此处，\alpha是迭代步长，G是替代H^{-1}的函数，g是一阶导数。 (2) \alpha计算：代码中(calc_alpha_by_ArmijoRule函数)使用的是线搜索技术·Armijo准则。该技术此前未提过，此处先给自己挖个坑，后续找时间再补上。 (3) G计算：代码中(update_D_by_BFGS函数)使用的是和梯度类算法原理：最速下降法、牛顿法和拟牛顿法中完全一致的过程。 import numpy as np def rosen(x): return sum(100 * (x[1:] - x[:-1] ** 2) ** 2 + (1 - x[:-1]) ** 2) def rosen_jac(x): xm = x[1:-1] xm_m1 = x[:-2] xm_p1 = x[2:] der = np.zeros_like(x) der[1:-1] = 200 * (xm - xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1 - xm) der[0] = -400 * x[0] * (x[1] - x[0] ** 2) - 2 * (1 - x[0]) der[-1] = 200 * (x[-1] - x[-2] ** 2) return der def calc_alpha_by_ArmijoRule(xCurr, fCurr, gCurr, dCurr, c=1.e-4, v=0.5): i = 0 alpha = v ** i xNext = xCurr + alpha * dCurr fNext = rosen(xNext) while True: if fNext 1e-6: # 迭代方向：GCurr * gCurr dCurr = -np.matmul(GCurr, gCurr) # 迭代步长策略：ArmijoRule alpha = calc_alpha_by_ArmijoRule(xCurr, fCurr, gCurr, dCurr) # 基本迭代公式：x_{k+1} = x_k + alpha * GCurr * gCurr xNext = xCurr + alpha * dCurr fNext = rosen(xNext) gNext = rosen_jac(xNext) # BFGS更新G GNext = update_D_by_BFGS(xCurr, gCurr, xNext, gNext, GCurr) xCurr, fCurr, gCurr, GCurr = xNext, fNext, gNext, GNext iterations += 1 print('iterations: {}, best_f: {}'.format(iterations, fCurr)) if __name__ == '__main__': np.random.seed(0) # 设置种子，保证结果可重复 x0 = np.random.rand(10) # 随机生成初始值 # 自编代码实现 calc_by_self(x0)运行代码后，可以得到如下结果。 显然，自编代码的最优解和使用scipy工具包得到的结果是一致的。 iterations: 207, best_f: [6.79531171e-18] |
【本文地址】
今日新闻 |
推荐新闻 |