👱🏿 🙎 🤲🏽 科学，优化 👰🏻 🕊️ 🏩

SciPy（发音为sai pie）是基于Numpy Python扩展的数学应用程序包。借助SciPy，交互式Python会话可为复杂系统（如MATLAB，IDL，Octave，R-Lab和SciLab）转变为相同的完整数据处理和原型开发环境。今天，我想简单地谈谈如何在scipy.optimize包中使用一些著名的优化算法。始终可以使用help（）命令或使用Shift + Tab获得有关使用功能的更详细和最新的帮助。

引言

为了使您自己和读者免于搜索和阅读源代码，指向方法描述的链接将主要链接到Wikipedia。通常，此信息足以理解常规方法及其应用条件。为了理解数学方法的本质，我们单击了指向更多权威出版物的链接，这些出版物可以在每篇文章的末尾或您喜欢的搜索引擎中找到。

因此，scipy.optimize模块包括以下过程的实现：

使用各种算法（Nelder-Mead单形，BFGS，共轭牛顿梯度， COBYLA和SLSQP ）对几个变量（最小值）的标量函数进行有条件和无条件的最小化
全局优化（例如： Basinhopping ， diff_evolution ）
最小二乘残差的最小化（least_squares）和将曲线拟合到非线性最小二乘的算法（curve_fit）
最小化一个变量的标量函数（minim_scalar）并找到根（root_scalar）
使用各种算法（混合鲍威尔，莱文贝格-马夸特或大规模方法，例如牛顿-克里洛夫）的方程组（根）的多维求解器。

在本文中，我们将只考虑整个列表中的第一项。

无条件最小化多个变量的标量函数

scipy.optimize包中的minim函数提供了一个通用接口，用于解决几个变量的标量函数的有条件和无条件最小化问题。为了演示其工作，我们将需要一个适合多个变量的函数，并将以不同的方式将其最小化。

为此，N个变量的Rosenbrock函数是完美的，其形式为：

左 （ 右 ） 左 （ 右 ） 左 （ 右 ）

$f \左（\ mathbf {x} \右）= \ sum_ {i = 1} ^ {N-1} [100 \左（x_ {i + 1} -x_ {i} ^ {2} \右） ^ {2} + \左（1-x_ {i} \右）^ {2}]$

尽管Rosenbrock函数及其Jacobi和Hessian矩阵（分别是一阶和二阶导数）已经在scipy.optimize包中定义了，但我们还是自己定义了它。

import numpy as np def rosen(x): """The Rosenbrock function""" return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)

为了清楚起见，我们在3D中绘制了两个变量的Rosenbrock函数的值。

渲染代码

 from mpl_toolkits.mplot3d import Axes3D import matplotlib.pyplot as plt from matplotlib import cm from matplotlib.ticker import LinearLocator, FormatStrFormatter #  3D  fig = plt.figure(figsize=[15, 10]) ax = fig.gca(projection='3d') #    ax.view_init(45, 30) #     X = np.arange(-2, 2, 0.1) Y = np.arange(-1, 3, 0.1) X, Y = np.meshgrid(X, Y) Z = rosen(np.array([X,Y])) #   surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm) plt.show()

预先知道最小值为0 $x_i = 1$ 请考虑使用各种scipy.optimize过程确定Rosenbrock函数的最小值的示例。

Nelder-Mead单纯形法（Nelder-Mead）

假设在5维空间中有一个初始点x0。使用Nelder-Mead单纯形算法（该算法指定为method参数的值）找到最接近Rosenbrock函数的最小点：

 from scipy.optimize import minimize x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) res = minimize(rosen, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 339 Function evaluations: 571 [1. 1. 1. 1. 1.]

单纯形法是最小化清晰定义且相当平滑的函数的最简单方法。它不需要计算函数的导数；仅指定其值就足够了。对于简单的最小化问题，Nelder-Mead方法是一个不错的选择。但是，由于它不使用梯度估计，因此可能需要更长的时间才能找到最小值。

鲍威尔法

仅计算函数值的另一种优化算法是Powell方法。要使用它，您需要在minim函数中设置method ='powell'。

 x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) res = minimize(rosen, x0, method='powell', options={'xtol': 1e-8, 'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 19 Function evaluations: 1622 [1. 1. 1. 1. 1.]

Broyden-Fletcher-Goldfarb-Channo算法（BFGS）

为了更快地收敛到解， BFGS过程使用目标函数的梯度。可以将梯度指定为函数或使用一阶差分来计算。无论如何，通常BFGS方法比单纯形方法需要更少的函数调用。

我们以解析形式找到Rosenbrock函数的导数：

部 分 部 分 和 限 制 （ ） （ ， ， ） （ ） ，

$\ frac {\部分f} {\部分x_j} = \和\限制_ {i = 1} ^ N 200（x_i-x_ {i-1} ^ 2）（\ delta_ {i，j}-2x_ {i -1，j}）-2（1-x_ {i-1}）\ delta_ {i-1，j} =$

（ ） （ ） （ ）

$= 200（x_j-x_ {j-1} ^ 2）-400x_j（x_ {j + 1}-x_j ^ 2）-2（1-x_j）$

该表达式对除第一个和最后一个定义为的所有变量的导数有效：

部 分 部 分 （ ） （ ） ，

$\ frac {\部分f} {\部分x_0} = -400 x_0（x_1-x_0 ^ 2）-2（1-x_0），$

$\ frac {\部分f} {\部分x_ {N-1}} = 200（x_ {N-1}-x_ {N-2} ^ 2）。$

让我们看一下计算此梯度的Python函数：

 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

梯度计算函数被指定为最小值函数的jac参数的值，如下所示。

 res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 25 Function evaluations: 30 Gradient evaluations: 30 [1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]

共轭梯度算法（牛顿）

牛顿共轭梯度算法是一种改进的牛顿方法。
牛顿法基于第二级多项式对局部区域中的函数进行逼近：

$f \左（\ mathbf {x} \右）\约f \左（\ mathbf {x} _ {0} \右）+ \ nabla f \左（\ mathbf {x} _ {0} \右） \ cdot \左（\ mathbf {x}-\ mathbf {x} _ {0} \右）+ \ frac {1} {2} \左（\ mathbf {x}-\ mathbf {x} _ {0} \ right）^ {T} \ mathbf {H} \ left（\ mathbf {x} _ {0} \ right）\ left（\ mathbf {x}-\ mathbf {x} _ {0} \ right）$

在哪里 $\ mathbf {H} \左（\ mathbf {x} _ {0} \右）$ 是二阶导数的矩阵（Hessian矩阵，Hessian）。
如果Hessian是正定的，则可以通过将二次形式的零梯度等于零来找到该函数的局部最小值。结果是一个表达式：

$\ mathbf {x} _ {\ textrm {opt}} = \ mathbf {x} _ {0}-\ mathbf {H} ^ {-1} \ nabla f$

使用共轭梯度法计算反Hessian。下面给出了使用此方法最小化Rosenbrock函数的示例。要使用Newton-CG方法，必须定义一个评估Hessian的函数。
分析形式的Rosenbrock函数的Hessian等于：

$H_ {i，j} = \ frac {\部分^ 2 f} {\部分x_i x_j} = 200（\ delta_ {i，j}-2x_ {i-1} \ delta {i-1，j}- 400x_i（\ delta_ {i + 1，j}-2x_i \ delta {i，j}）-400 \ delta_ {i，j}（x_ {i + 1}-x_i ^ 2）+ 2 \ delta_ {i，j } =$

$=（202 + 1200x_i ^ 2-400x_ {i + 1}）\ delta_ {i，j}-400x_i \ delta_ {i + 1，j}-400x_ {i-1} \ delta_ {i-1，j}$

在哪里 $i，j \ in \左[1，N-2 \右]$ 和 $i，j \ in \左[0，N-1 \右]$ 确定矩阵 $N \倍N$ 。

矩阵的其余非零元素等于：

$\ frac {\部分^ 2 f} {\部分x_0 ^ 2} = 1200x_0 ^ 2-400x_1 +2$

$\ frac {\部分^ 2 f} {\部分x_0 x_1} = \ frac {\部分^ 2 f} {\部分x_1 x_0} = -400x_0$

$\ frac {\部分^ 2 f} {\部分x_ {N-1} x_ {N-2}} = \ frac {\部分^ 2 f} {\部分x_ {N-2} x_ {N-1 }} = -400x_ {N-2}$

$\ frac {\部分^ 2 f} {\部分x_ {N-1} ^ 2} = 200x$

例如，在五维空间N = 5中，Rosenbrock函数的Hessian矩阵具有带形式：

$\ tiny \ mathbf {H} = \开始{bmatrix} 1200x_ {0} ^ {2} -400x_ {1} +2＆-400x_ {0}＆0＆0＆0 \\ -400x_ {0}＆202 + 1200x_ {1} ^ {2} -400x_ {2}＆-400x_ {1}＆0＆0 \\ 0＆-400x_ {1}＆202 + 1200x_ {2} ^ {2} -400x_ {3}＆ -400x_ {2}＆0 \\ 0＆＆-400x_ {2}＆202 + 1200x_ {3} ^ {2} -400x_ {4}＆-400x_ {3} \\ 0＆0＆0＆-400x_ { 3}和200 \结尾{bmatrix}$

使用共轭梯度方法（牛顿）计算此Hessian的代码以及使Rosenbrock函数最小化的代码：

 def rosen_hess(x): x = np.asarray(x) H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1) diagonal = np.zeros_like(x) diagonal[0] = 1200*x[0]**2-400*x[1]+2 diagonal[-1] = 200 diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:] H = H + np.diag(diagonal) return H res = minimize(rosen, x0, method='Newton-CG', jac=rosen_der, hess=rosen_hess, options={'xtol': 1e-8, 'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 24 Function evaluations: 33 Gradient evaluations: 56 Hessian evaluations: 24 [1. 1. 1. 0.99999999 0.99999999]

定义Hessian与任意向量的乘积函数的示例

在实际问题中，计算和存储整个黑森州矩阵可能需要大量的时间和内存资源。此外，实际上，由于Hessian矩阵本身没有必要指定，因为最小化过程只需要一个等于Hessian与另一个任意向量的乘积的向量。因此，从计算的角度来看，最好立即确定返回Hessian与任意向量的乘积结果的函数。

考虑hess函数，该函数将最小化向量作为第一个参数，并将任意向量作为第二个参数（以及最小化函数的其他参数）。在我们的情况下，计算Rosenbrock函数的Hessian与任意向量的乘积并不是很困难。如果p是任意向量，则乘积 $H（x）\ cdot p$ 具有以下形式：

$\ mathbf {H} \左（\ mathbf {x} \右）\ mathbf {p} = \开始{bmatrix} \左（1200x_ {0} ^ {2} -400x_ {1} +2 \ right）p_ {0} -400x_ {0} p_ {1} \\ \ vdots \\ -400x_ {i-1} p_ {i-1} + \左（202 + 1200x_ {i} ^ {2} -400x_ {i + 1} \对）p_ {i} -400x_ {i} p_ {i + 1} \\ \ vdots \\ -400x_ {N-2} p_ {N-2} + 200p_ {N-1} \ end {bmatrix }。$

计算hessian和任意向量乘积的函数将作为hessp参数的值传递给minimum函数：

 def rosen_hess_p(x, p): x = np.asarray(x) Hp = np.zeros_like(x) Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1] Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \ -400*x[1:-1]*p[2:] Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1] return Hp res = minimize(rosen, x0, method='Newton-CG', jac=rosen_der, hessp=rosen_hess_p, options={'xtol': 1e-8, 'disp': True})

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 24 Function evaluations: 33 Gradient evaluations: 56 Hessian evaluations: 66

共轭梯度的信赖域算法（牛顿）

Hessian矩阵的条件性差和错误的搜索方向会导致共轭牛顿梯度算法效率低下。在这种情况下，优先考虑共轭牛顿梯度的信赖域方法。

黑森州矩阵定义示例：

 res = minimize(rosen, x0, method='trust-ncg', jac=rosen_der, hess=rosen_hess, options={'gtol': 1e-8, 'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 20 Function evaluations: 21 Gradient evaluations: 20 Hessian evaluations: 19 [1. 1. 1. 1. 1.]

一个具有Hessian乘积函数和任意向量的示例：

 res = minimize(rosen, x0, method='trust-ncg', jac=rosen_der, hessp=rosen_hess_p, options={'gtol': 1e-8, 'disp': True}) print(res.x)

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 20 Function evaluations: 21 Gradient evaluations: 20 Hessian evaluations: 0 [1. 1. 1. 1. 1.]

Krylovsky类型方法

像trust-ncg方法一样，Krylovsky类型的方法仅使用矩阵向量乘积，因此非常适合解决大规模问题。它们的本质是解决受截断的Krylov子空间限制的机密领域中的问题。对于不确定的任务，最好使用此方法，因为与trust-ncg方法相比，由于每个子任务的矩阵向量乘积较少，因此使用的非线性迭代次数较少。此外，二次子任务的解决方案比trust-ncg方法更准确。
黑森州矩阵定义示例：

 res = minimize(rosen, x0, method='trust-krylov', jac=rosen_der, hess=rosen_hess, options={'gtol': 1e-8, 'disp': True}) Optimization terminated successfully. Current function value: 0.000000 Iterations: 19 Function evaluations: 20 Gradient evaluations: 20 Hessian evaluations: 18 print(res.x) [1. 1. 1. 1. 1.]

一个具有Hessian乘积函数和任意向量的示例：

 res = minimize(rosen, x0, method='trust-krylov', jac=rosen_der, hessp=rosen_hess_p, options={'gtol': 1e-8, 'disp': True}) Optimization terminated successfully. Current function value: 0.000000 Iterations: 19 Function evaluations: 20 Gradient evaluations: 20 Hessian evaluations: 0 print(res.x) [1. 1. 1. 1. 1.]

基于置信度的近似算法

所有方法（Newton-CG，trust-ncg和trust-krylov）都非常适合解决大规模任务（具有数千个变量）。这是由于以下事实：共轭梯度的基本算法意味着对逆黑森州矩阵的近似确定。该解决方案是迭代的，而没有Hessian的明确分解。由于仅需要确定Hessian和任意矢量的乘积的函数，因此该算法特别适用于处理稀疏（磁带对角线）矩阵。这样可以降低存储器成本并节省大量时间。

在中等大小的问题中，存储和分解Hessian的成本并不重要。这意味着可以通过更少的迭代获得解决方案，从而几乎完全解决了信任区域的子任务。为此，为每个二次子任务迭代求解了一些非线性方程。这种解决方案通常需要Holets Hessian矩阵的3或4个分解。结果，与置信域的其他实现方法相比，该方法收敛于较少的迭代，并且需要较少的目标函数计算。该算法仅意味着确定完整的Hessian矩阵，不支持使用Hessian乘积函数和任意向量的能力。

最小化Rosenbrock函数的示例：

 res = minimize(rosen, x0, method='trust-exact', jac=rosen_der, hess=rosen_hess, options={'gtol': 1e-8, 'disp': True}) res.x

 Optimization terminated successfully. Current function value: 0.000000 Iterations: 13 Function evaluations: 14 Gradient evaluations: 13 Hessian evaluations: 14 array([1., 1., 1., 1., 1.])

这也许住。在下一篇文章中，我将尝试介绍有关条件最小化的最有趣的东西，即最小化在解决逼近问题，最小化一个变量的功能，任意最小化器以及使用scipy.optimize程序包找到方程组根的应用。

资料来源： https : //docs.scipy.org/doc/scipy/reference/

科学，优化

引言