SciPy（发音为 saiie）是一个基于 Numpy Python 扩展的数学应用程序包。 借助 SciPy，您的交互式 Python 会话将成为与 MATLAB、IDL、Octave、R-Lab 和 SciLab 相同的完整数据科学和复杂系统原型环境。 今天我想简单讲一下如何使用scipy.optimize包中的一些著名的优化算法。 使用 help() 命令或使用 Shift+Tab 始终可以获得有关使用函数的更详细和最新的帮助。
介绍
为了使您和读者免于搜索和阅读主要来源，方法描述的链接将主要位于维基百科上。 通常，这些信息足以理解一般术语的方法及其应用条件。 要了解数学方法的本质，请点击更权威出版物的链接，这些出版物可以在每篇文章的末尾或您最喜欢的搜索引擎中找到。
因此，scipy.optimize模块包括以下过程的实现：
- 使用各种算法（Nelder-Mead 单纯形、BFGS、牛顿共轭梯度、 и )
- 全局优化（例如： , )
- 最大限度地减少残差 (least_squares) 和使用非线性最小二乘的曲线拟合算法 (curve_fit)
- 最小化一个变量的标量函数 (minim_scalar) 并搜索根 (root_scalar)
- 使用各种算法（混合 Powell、 或大规模方法，例如 ).
在本文中，我们将仅考虑整个列表中的第一项。
多个变量的标量函数的无条件最小化
scipy.optimize 包中的 minim 函数提供了一个通用接口，用于解决多个变量的标量函数的条件和无条件最小化问题。 为了演示它是如何工作的，我们需要一个包含多个变量的合适函数，我们将以不同的方式最小化它。
出于这些目的，N 个变量的 Rosenbrock 函数是完美的，其形式为：
尽管 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 ，让我们看一下如何使用各种 scipy.optimize 过程确定 Rosenbrock 函数的最小值的示例。
Nelder-Mead 单纯形法
设0维空间中有一个初始点x5。 让我们使用算法找到最接近它的 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 方法是解决简单最小化问题的不错选择。 然而，由于它不使用梯度估计，因此可能需要更长的时间才能找到最小值。
鲍威尔法
另一种只计算函数值的优化算法是 。 要使用它，您需要在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-Shanno (BFGS) 算法
为了更快地收敛到解决方案，该过程 使用目标函数的梯度。 梯度可以指定为函数或使用一阶差分来计算。 无论如何，BFGS 方法通常比单纯形方法需要更少的函数调用。
让我们以解析形式求 Rosenbrock 函数的导数：
该表达式对于除第一个和最后一个变量之外的所有变量的导数均有效，其定义为：
让我们看一下计算这个梯度的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梯度计算函数指定为minim函数的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]共轭梯度算法（牛顿）
算法 是一种改进的牛顿法。
牛顿方法基于用二阶多项式逼近局部区域中的函数：
哪里 是二阶导数矩阵（Hessian 矩阵，Hessian）。
如果 Hessian 矩阵是正定的，则可以通过将二次形式的零梯度等于零来找到该函数的局部最小值。 结果将是表达式：
使用共轭梯度法计算逆 Hessian 矩阵。 下面给出了使用此方法最小化 Rosenbrock 函数的示例。 要使用 Newton-CG 方法，您必须指定计算 Hessian 的函数。
Rosenbrock 函数的 Hessian 矩阵的解析形式等于：
哪里 и ，定义矩阵 .
矩阵的其余非零元素等于：
例如，在五维空间 N = 5 中，Rosenbrock 函数的 Hessian 矩阵具有带状形式：
计算该 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 矩阵与另一个任意向量的乘积的向量。 因此，从计算的角度来看，最好立即定义一个返回 Hessian 矩阵与任意向量的乘积结果的函数。
考虑 hess 函数，它将最小化向量作为第一个参数，将任意向量作为第二个参数（以及要最小化的函数的其他参数）。 在我们的例子中，计算 Rosenbrock 函数的 Hessian 矩阵与任意向量的乘积并不是很困难。 如果 p 是任意向量，则乘积 具有以下形式：
计算 Hessian 矩阵和任意向量的乘积的函数作为 hessp 参数的值传递给最小化函数：
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 矩阵的条件不佳和搜索方向不正确可能会导致牛顿共轭梯度算法无效。 在这种情况下，优先考虑 （信任区域）共轭牛顿梯度。
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.]克雷洛夫型方法
与 trust-ncg 方法一样，Krylov 类型的方法非常适合解决大规模问题，因为它们仅使用矩阵向量乘积。 它们的本质是解决由截断的克雷洛夫子空间限制的置信区域中的问题。 对于不确定问题，最好使用此方法，因为与 trust-ncg 方法相比，由于每个子问题的矩阵向量乘积数量较少，因此它使用的非线性迭代次数较少。 此外，二次子问题的解比使用 trust-ncg 方法更准确。
Hessian 矩阵定义的示例：
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 矩阵和任意向量的乘积定义一个函数，因此该算法特别适合处理稀疏（带对角）矩阵。 这可以降低内存成本并节省大量时间。
对于中等规模的问题，存储和分解 Hessian 矩阵的成本并不重要。 这意味着可以通过更少的迭代获得解决方案，几乎精确地解决信任域的子问题。 为此，需要对每个二次子问题迭代求解一些非线性方程。 这样的解决方案通常需要 Hessian 矩阵的 3 或 4 次 Cholesky 分解。 因此，与其他实现的置信区域方法相比，该方法以更少的迭代次数收敛，并且需要更少的目标函数计算。 该算法仅涉及确定完整的Hessian矩阵，不支持使用Hessian与任意向量的乘积函数的能力。
Rosenbrock 函数最小化的示例：
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimization 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 查找方程组的根包裹。
来源：
来源： habr.com