拟牛顿法

在梯度类算法原理：最速下降法、牛顿法和拟牛顿法中，介绍了梯度类算法求解优化问题的设计思路，并以最速下降法、牛顿法和拟牛顿法为例，描述了具体的算法迭代过程。其中，拟牛顿法（Broyden–Fletcher–Goldfarb–Shanno，BFGS）在实际优化场景中被广泛使用，因此本文将自主编写Python代码，实现BFGS的全部过程，并且和现有工具包做对比，从而加深对BFGS的理解。

待优化实例

本文使用的待优化实例为10维的Rosenbrock函数：
f(x)=∑i=210[100(xi−xi−12)2+(1−xi−1)2]f(x)=\\sum_{i=2}^{10} [100(x_i-x_{i-1}^2)^2+(1-x_{i-1})^2]f(x)=i=2∑10​[100(xi​−xi−12​)2+(1−xi−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，宽6plt.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工具包实现BFGS

python的scipy.optimize工具包能够实现BFGS，具体可以参考scipy.optimize优化器的各种使用。本节仅做简要描述。

主要的求解方案有两种：第一种是只将待优化函数和初值传递给工具包，其余内容全部由工具包自行计算；第二种是将待优化函数的雅可比矩阵也传递给工具包，提升求解的速度。

雅可比矩阵之前没提到过，这里解释一下。雅可比矩阵是一阶偏导数以一定方式排列成的矩阵。假设f=[f1,f2,⋅⋅⋅,fm]\\pmb f=[f_1,f_2,···,f_m]f=[f1​,f2​,⋅⋅⋅,fm​]，x=[x1,x2,⋅⋅⋅,xn]\\pmb x=[x_1,x_2,···,x_n]x=[x1​,x2​,⋅⋅⋅,xn​]，那么$f$的雅可比矩阵 $J$ 为 $m\times n$ 的矩阵:
J=[∂f∂x1⋅⋅⋅∂f∂xn]=[∂f1∂x1⋯∂f1∂xn⋮⋱⋮∂fm∂x1⋯∂fm∂xn]\\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]J=[∂x1​∂f​⋅⋅⋅∂xn​∂f​]=​∂x1​∂f1​​⋮∂x1​∂fm​​​⋯⋱⋯​∂xn​∂f1​​⋮∂xn​∂fm​​​​

针对Rosenbrock函数，一阶偏导数为

∂f∂xj=200(xj−xj−12)−400xj(xj+1−xj2)−2(1−xj)\\frac{\\partial f}{\\partial x_j} =200(x_j-x_{j-1}^2)-400x_j(x_{j+1}-x_j^2)-2(1-x_j)∂xj​∂f​=200(xj​−xj−12​)−400xj​(xj+1​−xj2​)−2(1−xj​)
当$j=1$时
∂f∂x1=−400x1(x2−x12)−2(1−x1)\\frac{\\partial f}{\\partial x_1}=-400x_1(x_2-x_1^2)-2(1-x_1)∂x1​∂f​=−400x1​(x2​−x12​)−2(1−x1​)
当$j=10$时
∂f∂x10=200(x10−x92)\\frac{\\partial f}{\\partial x_{10}}=200(x_{10}-x_9^2)∂x10​∂f​=200(x10​−x92​)

两种求解方案的代码如下：

import numpy as np
from scipy.optimize import minimize
import timedef 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 derif __name__ == '__main__':np.random.seed(0)  # 设置种子，保证结果可重复x0 = np.random.rand(10)  # 随机生成初始值print('------------不提供jac计算--------------')res = minimize(rosen, x0, options={'disp': True})  # 不提供jacprint('-------------提供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_timecalc_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.000000Iterations: 41Function evaluations: 572Gradient evaluations: 52
-------------提供jac计算---------------
Optimization terminated successfully.Current function value: 0.000000Iterations: 42Function evaluations: 52Gradient evaluations: 52
-------------评估jac效率提升---------------
不提供jac时，计算时间为：3.727977991104126 秒; 提供jac后，计算时间为：1.1221048831939697 秒，时间降低为原来的：0.3009955761197058

自编Python实现BFGS

自行编写代码时，主要关注三项内容：

(1) 迭代公式：xk+1=xk+α∗G∗gx_{k+1} = x_k + \\alpha * G * gxk+1​=xk​+α∗G∗g。此处，$\alpha$是迭代步长，$G$是替代H−1H^{-1}H−1的函数，$g$是一阶导数。

(2) $\alpha$计算：代码中(calc_alpha_by_ArmijoRule函数)使用的是线搜索技术·Armijo准则。该技术此前未提过，此处先给自己挖个坑，后续找时间再补上。

(3) $G$计算：代码中(update_D_by_BFGS函数)使用的是和梯度类算法原理：最速下降法、牛顿法和拟牛顿法中完全一致的过程。

import numpy as npdef 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 derdef calc_alpha_by_ArmijoRule(xCurr, fCurr, gCurr, dCurr, c=1.e-4, v=0.5):i = 0alpha = v ** ixNext = xCurr + alpha * dCurrfNext = rosen(xNext)while True:if fNext <= fCurr + c * alpha * np.matmul(dCurr.T, gCurr):breaki += 1alpha = v ** ixNext = xCurr + alpha * dCurrfNext = rosen(xNext)return alphadef update_D_by_BFGS(xCurr, gCurr, xNext, gNext, GCurr):sk = xNext - xCurryk = gNext - gCurrterm1 = GCurr @ yk @ yk.T @ GCurr.T / (yk.T @ GCurr.T @ yk)term2 = sk @ sk.T / (sk.T @ yk)Dk = GCurr - term1 + term2return Dkdef calc_by_self(initial_x):# x：变量, f：目标函数, g：目标函数一阶导数, D：H^(-1)的替代表达式xCurr = initial_x.reshape((len(initial_x), 1))fCurr = rosen(xCurr)gCurr = rosen_jac(xCurr).reshape((len(xCurr), 1))GCurr = np.identity(len(xCurr))iterations = 0# 一阶导数不为0，则持续迭代while np.linalg.norm(gCurr) > 1e-6:# 迭代方向：GCurr * gCurrdCurr = -np.matmul(GCurr, gCurr)# 迭代步长策略：ArmijoRulealpha = calc_alpha_by_ArmijoRule(xCurr, fCurr, gCurr, dCurr)# 基本迭代公式：x_{k+1} = x_k + alpha * GCurr * gCurrxNext = xCurr + alpha * dCurrfNext = rosen(xNext)gNext = rosen_jac(xNext)# BFGS更新GGNext = update_D_by_BFGS(xCurr, gCurr, xNext, gNext, GCurr)xCurr, fCurr, gCurr, GCurr = xNext, fNext, gNext, GNextiterations += 1print('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]