优化 (scipy.optimize)#
scipy.optimize 包提供了几种常用的优化算法。详细列表可参见：
scipy.optimize`（也可通过 ``help(scipy.optimize)` 找到）。
多元标量函数无约束最小化 (minimize)#
minimize 函数为 scipy.optimize 中的多元标量函数提供了无约束和约束最小化算法的通用接口。为了演示最小化函数，考虑最小化具有 \(N\) 个变量的 Rosenbrock 函数的问题：
该函数的最小值为 0，当 \(x_{i}=1.\) 时达到。
注意，Rosenbrock 函数及其导数包含在 scipy.optimize 中。以下各节中的实现展示了如何定义目标函数及其雅可比矩阵和 Hessian 矩阵。scipy.optimize 中的目标函数期望其第一个参数为要优化的 numpy 数组，并返回一个浮点值。确切的调用签名必须是 f(x, *args)，其中 x 表示一个 numpy 数组，args 是一个提供给目标函数的额外参数的元组。
Nelder-Mead 单纯形算法 (method='Nelder-Mead')#
在下面的示例中，minimize 例程与 Nelder-Mead 单纯形算法一起使用（通过 method 参数选择）：
>>> import numpy as np >>> from scipy.optimize import minimize>>> def rosen(x):… “””The Rosenbrock function””” … return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)
>>> x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2]) >>> res = minimize(rosen, x0, method='nelder-mead', ... options={'xatol': 1e-8, 'disp': True}) 优化成功终止。 当前函数值: 0.000000 迭代次数: 339 函数评估次数: 571>>> print(res.x) [1. 1. 1. 1. 1.]
单纯形算法可能是最简单的最小化一个相当良好行为函数的方法。它只需要函数评估，是简单最小化问题的良好选择。然而，由于它不使用任何梯度评估，可能需要更长时间才能找到最小值。
另一个只需要函数调用就能找到最小值的优化算法是 Powell 方法，可以通过在 minimize 中设置 method='powell' 来使用。
为了演示如何向目标函数提供额外参数，让我们用一个额外的缩放因子 a 和一个偏移量 b 来最小化 Rosenbrock 函数：
再次使用 minimize 例程，可以通过以下代码块解决这个问题，示例参数为 a=0.5 和 b=1。
>>> def rosen_with_args(x, a, b):
... """The Rosenbrock function with additional arguments"""
... return sum(a*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) + b
>>> x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
>>> res = minimize(rosen_with_args, x0, method='nelder-mead',
... args=(0.5, 1.), options={'xatol': 1e-8, 'disp': True})
优化成功终止。
当前函数值: 1.000000
迭代次数: 319 # 可能会有所不同
函数评估次数: 525 # 可能会有所不同
>>> print(res.x)
[1. 1. 1. 1. 0.99999999]
作为使用 minimize 的 args 参数的替代方法，只需将目标函数包装在一个新函数中，该函数仅接受 x。当需要将其他参数作为关键字参数传递给目标函数时，这种方法也很有用。
>>> def rosen_with_args(x, a, *, b): # b 是一个仅关键字参数
... return sum(a*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0) + b
>>> def wrapped_rosen_without_args(x):
... return rosen_with_args(x, 0.5, b=1.) # 传入 `a` 和 `b`
>>> x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
>>> res = minimize(wrapped_rosen_without_args, x0, method='nelder-mead',
... options={'xatol': 1e-8,})
>>> print(res.x)
[1. 1. 1. 1. 0.99999999]
另一种替代方法是使用 functools.partial。
>>> from functools import partial
>>> partial_rosen = partial(rosen_with_args, a=0.5, b=1.)
>>> res = minimize(partial_rosen, x0, method='nelder-mead',
... options={'xatol': 1e-8,})
>>> print(res.x)
[1. 1. 1. 1. 0.99999999]
Broyden-Fletcher-Goldfarb-Shanno 算法 (method='BFGS')#
为了更快地收敛到解，此例程使用目标函数的梯度。如果用户没有提供梯度，则使用一阶差分进行估计。Broyden-Fletcher-Goldfarb-Shanno (BFGS) 方法通常在必须估计梯度时需要的函数调用次数比单纯形算法少。
为了演示此算法，再次使用 Rosenbrock 函数。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
此梯度信息通过``jac``参数在:func:`minimize`函数中指定，如下所示。
>>> res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
... options={'disp': True})
优化成功终止。
当前函数值: 0.000000
迭代次数: 25 # 可能会有所不同
函数评估次数: 30
梯度评估次数: 30
>>> res.x
array([1., 1., 1., 1., 1.])
避免冗余计算#
目标函数及其梯度通常会共享计算的一部分。例如，考虑以下问题。
>>> def f(x):
... return -expensive(x[0])**2
>>>
>>> def df(x):
... return -2 * expensive(x[0]) * dexpensive(x[0])
>>>
>>> def expensive(x):
... # 这个函数计算开销很大！
... expensive.count += 1 # 让我们记录它运行了多少次
... return np.sin(x)
>>> expensive.count = 0
>>>
>>> def dexpensive(x):
... return np.cos(x)
>>>
>>> res = minimize(f, 0.5, jac=df)
>>> res.fun
-0.9999999999999174
>>> res.nfev, res.njev
6, 6
>>> expensive.count
12
在这里，expensive 被调用了12次：目标函数中调用了6次，梯度中调用了6次。减少冗余计算的一种方法是创建一个同时返回目标函数和梯度的单一函数。
>>> def f_and_df(x):
... expensive_value = expensive(x[0])
... return (-expensive_value**2, # 目标函数
... -2*expensive_value*dexpensive(x[0])) # 梯度
>>>
>>> expensive.count = 0 # 重置计数器
>>> res = minimize(f_and_df, 0.5, jac=True)
>>> res.fun
-0.9999999999999174
>>> expensive.count
6
当我们调用 minimize 时，我们指定 jac==True 以表明提供的函数返回目标函数及其梯度。虽然方便，但并非所有 scipy.optimize 函数都支持此功能，而且它仅用于在函数及其梯度之间共享计算，而在某些问题中，我们希望与 Hessian（目标函数的二阶导数）和约束共享计算。一种更通用的方法是记忆化计算中的昂贵部分。在简单情况下，这可以通过 functools.lru_cache 包装器来实现。
>>> from functools import lru_cache
>>> expensive.count = 0 # 重置计数器
>>> expensive = lru_cache(expensive)
>>> res = minimize(f, 0.5, jac=df)
>>> res.fun
-0.9999999999999174
>>> expensive.count
6
牛顿共轭梯度算法（method='Newton-CG'）#
牛顿共轭梯度算法是一种修正的牛顿法，它使用共轭梯度算法来（近似）求逆局部海森矩阵 [NW]。牛顿法基于将函数局部拟合为二次形式：
其中 \(\mathbf{H}\left(\mathbf{x}_{0}\right)\) 是二阶导数矩阵（海森矩阵）。如果海森矩阵是正定的，则可以通过将二次形式的梯度设为零来找到该函数的局部最小值，从而得到
海森矩阵的逆通过共轭梯度方法进行求解。以下是使用此方法最小化Rosenbrock函数的示例。为了充分利用牛顿-CG方法，必须提供一个计算海森矩阵的函数。海森矩阵本身不需要构建，只需要一个向量，该向量是海森矩阵与任意向量的乘积，需要提供给最小化例程。因此，用户可以提供一个计算海森矩阵的函数，或者一个计算海森矩阵与任意向量乘积的函数。
完整海森矩阵示例：#
Rosenbrock函数的海森矩阵为
如果 \(i,j\in\left[1,N-2\right]\) 且 \(i,j\in\left[0,N-1\right]\) 定义了 \(N\times N\) 矩阵。矩阵的其他非零项为
例如，当 \(N=5\) 时的 Hessian 矩阵为
计算此 Hessian 矩阵的代码以及使用 Newton-CG 方法最小化函数的代码如下所示：
>>> 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}) 优化成功终止。
当前函数值: 0.000000 迭代次数: 19 # 可能会有所不同 函数评估次数: 22 梯度评估次数: 19 海森矩阵评估次数: 19
>>> res.x array([1., 1., 1., 1., 1.])
海森矩阵乘积示例:#
对于较大的最小化问题，存储整个海森矩阵可能会消耗相当多的时间和内存。牛顿-CG算法只需要海森矩阵与任意向量的乘积。因此，用户可以通过提供一个``hess``函数来计算这个乘积，而不是完整的海森矩阵，该函数将最小化向量作为第一个参数，将任意向量作为第二个参数（以及传递给要最小化的函数的额外参数）。如果可能的话，使用带有海森矩阵乘积选项的牛顿-CG算法可能是最快的方式来最小化函数。
在这种情况下，Rosenbrock函数的海森矩阵与任意向量的乘积并不难计算。如果 \(\mathbf{p}\) 是任意向量，那么 \(\mathbf{H}\left(\mathbf{x}\right)\mathbf{p}\) 的元素为:
以下代码利用这个海森矩阵乘积来使用 minimize 最小化Rosenbrock函数:
>>> 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})
优化成功终止。
当前函数值: 0.000000
迭代次数: 20 # 可能会有所不同
函数评估次数: 23
梯度评估次数: 20
海森评估次数: 44
>>> res.x
array([1., 1., 1., 1., 1.])
根据 [NW] p. 170，当海森矩阵条件数较差时，Newton-CG 算法可能会效率低下，因为在这种情况下，该方法提供的搜索方向质量较差。根据作者的说法，trust-ncg 方法更有效地处理这种问题情况，接下来将对其进行描述。
信赖域牛顿共轭梯度算法 (method='trust-ncg')#
Newton-CG 方法是一种线搜索方法：它找到一个最小化函数二次近似的搜索方向，然后使用线搜索算法在该方向上找到（近乎）最优的步长。另一种方法是，首先固定步长限制 \(\Delta\)，然后在给定的信赖域内通过求解以下二次子问题来找到最优步长 \(\mathbf{p}\)：
然后更新解 \(\mathbf{x}_{k+1} = \mathbf{x}_{k} + \mathbf{p}\)，并根据二次模型与真实函数的吻合程度调整信赖域 \(\Delta\)。这类方法被称为信赖域方法。
trust-ncg 算法是一种信赖域方法，它使用共轭梯度算法来求解信赖域子问题 [NW]。
完整 Hessian 矩阵示例：#
>>> res = minimize(rosen, x0, method='trust-ncg',
... jac=rosen_der, hess=rosen_hess,
... options={'gtol': 1e-8, 'disp': True})
优化成功终止。
当前函数值: 0.000000
迭代次数: 20 # 可能会有所不同
函数评估次数: 21
梯度评估次数: 20
Hessian 评估次数: 19
>>> res.x
array([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})
优化成功终止。
当前函数值: 0.000000
迭代次数: 20 # 可能会有所不同
函数评估次数: 21
梯度评估次数: 20
Hessian 评估次数: 0
>>> res.x
array([1., 1., 1., 1., 1.])
信赖域截断广义 Lanczos / 共轭梯度算法 (method='trust-krylov')#
与 trust-ncg 方法类似，trust-krylov 方法也适用于大规模问题，因为它仅通过矩阵-向量乘积的方式使用 Hessian 矩阵。
它比 trust-ncg 方法更精确地求解二次子问题。
此方法封装了 [TRLIB] 实现的 [GLTR] 方法，该方法精确求解了受限于截断 Krylov 子空间的信赖域子问题。对于不定问题，通常使用此方法会更好，因为它减少了非线性迭代的次数，代价是与 trust-ncg 方法相比，每个子问题求解中多了一些矩阵向量乘积。
完整 Hessian 矩阵示例：#
>>> res = minimize(rosen, x0, method='trust-krylov',
... jac=rosen_der, hess=rosen_hess,
... options={'gtol': 1e-8, 'disp': True})
优化成功终止。
当前函数值: 0.000000
迭代次数: 19 # 可能会有所不同
函数评估次数: 20
梯度评估次数: 20
Hessian 评估次数: 18
>>> res.x
array([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})
优化成功终止。
当前函数值: 0.000000
迭代次数: 19 # 可能会有所不同
函数评估次数: 20
梯度评估次数: 20
Hessian 评估次数: 0
>>> res.x
array([1., 1., 1., 1., 1.])
F. Lenders, C. Kirches, A. Potschka: “trlib: A vector-free implementation of the GLTR method for iterative solution of the trust region problem”, arXiv:1611.04718
N. Gould, S. Lucidi, M. Roma, P. Toint: “Solving the Trust-Region Subproblem using the Lanczos Method”, SIAM J. Optim., 9(2), 504–525, (1999). DOI:10.1137/S1052623497322735
信赖域近似精确算法 (method='trust-exact')#
所有方法 Newton-CG、trust-ncg 和 trust-krylov 都适合处理大规模问题（具有数千个变量的问题）。这是因为共轭梯度算法通过迭代近似解决信赖域子问题（或反转 Hessian 矩阵），而无需显式 Hessian 矩阵的分解。由于只需要 Hessian 矩阵与任意向量的乘积，该算法特别适合处理稀疏 Hessian 矩阵，从而减少存储需求并显著节省稀疏问题的时间。
对于中等规模的问题，其 Hessian 矩阵的存储和分解成本不是关键因素，可以通过几乎精确地解决信赖域子问题，在更少的迭代次数内获得解。为此，对于每个二次子问题，通过迭代求解某个非线性方程 [CGT]。该解通常需要对 Hessian 矩阵进行 3 或 4 次 Cholesky 分解。因此，该方法在较少的迭代次数内收敛，并且比其他实现的信赖域方法需要更少的评估目标函数的次数。该算法不支持 Hessian 乘积选项。以下是使用 Rosenbrock 函数的示例：
>>> res = minimize(rosen, x0, method='trust-exact',
... jac=rosen_der, hess=rosen_hess,
... options={'gtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 13 # 可能会有所不同
Function evaluations: 14
Gradient evaluations: 13
Hessian evaluations: 14
>>> res.x
array([1., 1., 1., 1., 1.])
多元标量函数的有约束最小化（minimize）#
minimize 函数提供了用于有约束最小化的算法，即 'trust-constr'、'SLSQP'、'COBYLA' 和 'COBYQA'。它们要求使用稍微不同的结构来定义约束条件。方法 'trust-constr' 和 'COBYQA' 要求将约束条件定义为一系列对象 LinearConstraint 和 NonlinearConstraint。而方法 'SLSQP' 和 'COBYLA' 则要求将约束条件定义为一系列字典，字典的键为 type、fun 和 jac。
作为一个例子，我们考虑 Rosenbrock 函数的有约束最小化问题：
这个优化问题有唯一解 \([x_0, x_1] = [0.4149,~ 0.1701]\)，其中只有第一个和第四个约束是有效的。
信赖域约束算法（method='trust-constr'）#
信赖域约束方法处理以下形式的有约束最小化问题：
当 \(c^l_j = c^u_j\) 时，该方法将第 \(j\) 个约束视为等式约束并相应处理。除此之外，单边约束
可以通过将上下界设置为 np.inf 并加上适当的符号来指定。
该实现基于 [EQSQP] 用于等式约束问题，基于 [TRIP] 用于不等式约束问题。两者都是适用于大规模问题的信赖域类型算法。
定义边界约束:#
边界约束 \(0 \leq x_0 \leq 1\) 和 \(-0.5 \leq x_1 \leq 2.0\) 使用 Bounds 对象定义。
>>> from scipy.optimize import Bounds
>>> bounds = Bounds([0, -0.5], [1.0, 2.0])
定义线性约束:#
约束 \(x_0 + 2 x_1 \leq 1\) 和 \(2 x_0 + x_1 = 1\) 可以写成线性约束的标准格式:
并使用 LinearConstraint 对象定义。
>>> from scipy.optimize import LinearConstraint
>>> linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1])
定义非线性约束:#
非线性约束:
其雅可比矩阵为:
以及海森矩阵的线性组合:
使用 NonlinearConstraint 对象定义。
>>> def cons_f(x):… return [x[0]**2 + x[1], x[0]**2 - x[1]] >>> def cons_J(x): … return [[2*x[0], 1], [2*x[0], -1]] >>> def cons_H(x, v): … return v[0]*np.array([[2, 0], [0, 0]]) + v[1]*np.array([[2, 0], [0, 0]]) >>> from scipy.optimize import NonlinearConstraint >>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=cons_H)
或者，也可以将 Hessian \(H(x, v)\) 定义为稀疏矩阵，
>>> from scipy.sparse import csc_matrix
>>> def cons_H_sparse(x, v):
... return v[0]*csc_matrix([[2, 0], [0, 0]]) + v[1]*csc_matrix([[2, 0], [0, 0]])
>>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
... jac=cons_J, hess=cons_H_sparse)
或者定义为 LinearOperator 对象。
>>> from scipy.sparse.linalg import LinearOperator
>>> def cons_H_linear_operator(x, v):
... def matvec(p):
... return np.array([p[0]*2*(v[0]+v[1]), 0])
... return LinearOperator((2, 2), matvec=matvec)
>>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1,
... jac=cons_J, hess=cons_H_linear_operator)
当 Hessian \(H(x, v)\) 的计算难以实现或计算上不可行时，可以使用 HessianUpdateStrategy。
目前可用的策略是 BFGS 和 SR1。
>>> from scipy.optimize import BFGS
>>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess=BFGS())
或者，可以使用有限差分来近似 Hessian。
>>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac=cons_J, hess='2-point')
约束的 Jacobian 也可以通过有限差分来近似。然而，在这种情况下，
Hessian 不能通过有限差分计算，需要使用其他方法。
由用户提供或使用 HessianUpdateStrategy 定义。
>>> nonlinear_constraint = NonlinearConstraint(cons_f, -np.inf, 1, jac='2-point', hess=BFGS())
求解优化问题:#
使用以下方法求解优化问题：
>>> x0 = np.array([0.5, 0])
>>> res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess,
... constraints=[linear_constraint, nonlinear_constraint],
... options={'verbose': 1}, bounds=bounds)
# 可能会有所不同
`gtol` 终止条件已满足。
迭代次数: 12, 函数评估次数: 8, CG 迭代次数: 7, 最优性: 2.99e-09, 约束违反: 1.11e-16, 执行时间: 0.016 s.
>>> print(res.x)
[0.41494531 0.17010937]
在需要时，目标函数的 Hessian 矩阵可以使用 LinearOperator 对象定义，
>>> def rosen_hess_linop(x):
... def matvec(p):
... return rosen_hess_p(x, p)
... return LinearOperator((2, 2), matvec=matvec)
>>> res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hess=rosen_hess_linop,
... constraints=[linear_constraint, nonlinear_constraint],
... options={'verbose': 1}, bounds=bounds)
# 可能会有所不同
`gtol` 终止条件已满足。
迭代次数: 12, 函数评估次数: 8, CG 迭代次数: 7, 最优性: 2.99e-09, 约束违反: 1.11e-16, 执行时间: 0.018 s.
>>> print(res.x)
[0.41494531 0.17010937]
或者通过参数 hessp 提供 Hessian-vector 乘积。
>>> res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_p,
... constraints=[linear_constraint, nonlinear_constraint],
... options={'verbose': 1}, bounds=bounds)
# 可能会有所不同
`gtol` 终止条件已满足。
迭代次数: 12, 函数评估次数: 8, CG 迭代次数: 7, 最优性: 2.99e-09, 约束违反: 1.11e-16, 执行时间: 0.018 s.
>>> print(res.x)
[0.41494531 0.17010937]
迭代次数：12，函数评估次数：8，CG迭代次数：7，最优性：2.99e-09，约束违反：1.11e-16，执行时间：0.018秒。
>>> print(res.x)
[0.41494531 0.17010937]
或者，可以近似目标函数的一阶和二阶导数。
例如，可以使用 SR1 拟牛顿近似法来近似 Hessian 矩阵，并使用有限差分来近似梯度。
>>> from scipy.optimize import SR1
>>> res = minimize(rosen, x0, method='trust-constr', jac="2-point", hess=SR1(),
... constraints=[linear_constraint, nonlinear_constraint],
... options={'verbose': 1}, bounds=bounds)
# 可能会有所不同
`gtol` 终止条件已满足。
迭代次数：12，函数评估次数：24，CG迭代次数：7，最优性：4.48e-09，约束违反：0.00e+00，执行时间：0.016秒。
>>> print(res.x)
[0.41494531 0.17010937]
序列最小二乘规划（SLSQP）算法（method='SLSQP'）#
SLSQP 方法处理以下形式的约束最小化问题：
其中 \(\mathcal{E}\) 或 \(\mathcal{I}\) 是包含等式和不等式约束的索引集合。
线性和非线性约束都定义为字典，包含键 type、fun 和 jac。
>>> ineq_cons = {'type': 'ineq',
... 'fun' : lambda x: np.array([1 - x[0] - 2*x[1],
... 1 - x[0]**2 - x[1],
... 1 - x[0]**2 + x[1]]),
... 'jac' : lambda x: np.array([[-1.0, -2.0],
... [-2*x[0], -1.0],
... [-2*x[0], 1.0]])}
>>> eq_cons = {'type': 'eq',
... 'fun' : lambda x: np.array([2*x[0] + x[1] - 1]),
... 'jac' : lambda x: np.array([2.0, 1.0])}
然后使用以下方法求解优化问题：
>>> x0 = np.array([0.5, 0])
>>> res = minimize(rosen, x0, method='SLSQP', jac=rosen_der,
... constraints=[eq_cons, ineq_cons], options={'ftol': 1e-9, 'disp': True},
... bounds=bounds)
# 结果可能有所不同
优化成功终止。 (退出模式 0)
当前函数值: 0.342717574857755
迭代次数: 5
函数评估次数: 6
梯度评估次数: 5
>>> print(res.x)
[0.41494475 0.1701105 ]
大多数适用于方法 'trust-constr' 的选项在 'SLSQP' 中不可用。
全局优化#
全局优化旨在在给定范围内找到函数的全局最小值，即使存在许多局部最小值。通常，全局最小化器高效地搜索参数空间，同时在底层使用局部最小化器（例如，minimize）。SciPy 包含许多优秀的全局优化器。在这里，我们将在相同的优化目标函数上使用这些优化器，即（恰如其名）``eggholder`` 函数:
>>> def eggholder(x):
... return (-(x[1] + 47) * np.sin(np.sqrt(abs(x[0]/2 + (x[1] + 47))))
... -x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47)))))
>>> bounds = [(-512, 512), (-512, 512)]
这个函数看起来像一个蛋盒:
>>> import matplotlib.pyplot as plt
>>> from mpl_toolkits.mplot3d import Axes3D
>>> x = np.arange(-512, 513)
>>> y = np.arange(-512, 513)
>>> xgrid, ygrid = np.meshgrid(x, y)
>>> xy = np.stack([xgrid, ygrid])
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111, projection='3d')
>>> ax.view_init(45, -45)
>>> ax.plot_surface(xgrid, ygrid, eggholder(xy), cmap='terrain')
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('y')
>>> ax.set_zlabel('eggholder(x, y)')
>>> plt.show()
我们现在使用全局优化器来获得最小值和函数在最小值处的值。我们将结果存储在一个字典中，以便稍后比较不同的优化结果。
>>> from scipy import optimize
>>> results = dict()
>>> results['shgo'] = optimize.shgo(eggholder, bounds)
>>> results['shgo']
fun: -935.3379515604197 # 可能会有所不同
funl: array([-935.33795156])
message: 'Optimization terminated successfully.'
nfev: 42
nit: 2
nlfev: 37
nlhev: 0
nljev: 9
success: True
x: array([439.48096952, 453.97740589])
xl: array([[439.48096952, 453.97740589]])
>>> results['DA'] = optimize.dual_annealing(eggholder, bounds)
>>> results['DA']
fun: -956.9182316237413 # 可能会有所不同
message: ['Maximum number of iteration reached']
nfev: 4091
nhev: 0
nit: 1000
njev: 0
x: array([482.35324114, 432.87892901])
所有优化器返回一个 OptimizeResult，除了解之外，还包含有关函数评估次数、优化是否成功等信息。为了简洁起见，我们不会展示其他优化器的完整输出:
>>> results['DE'] = optimize.differential_evolution(eggholder, bounds)
shgo 有第二种方法，它返回所有局部最小值，而不仅仅是它认为的全局最小值:
>>> results['shgo_sobol'] = optimize.shgo(eggholder, bounds, n=200, iters=5,
... sampling_method='sobol')
我们现在将在函数的热图中绘制所有找到的最小值:
>>> fig = plt.figure()
>>> ax = fig.add_subplot(111)
>>> im = ax.imshow(eggholder(xy), interpolation='bilinear', origin='lower',
... cmap='gray')
>>> ax.set_xlabel('x')
>>> ax.set_ylabel('y')
>>>
>>> def plot_point(res, marker='o', color=None):
... ax.plot(512+res.x[0], 512+res.x[1], marker=marker, color=color, ms=10)
>>> plot_point(results['DE'], color='c') # differential_evolution - 青色
>>> plot_point(results['DA'], color='w') # dual_annealing. - 白色
>>> # SHGO 产生多个最小值，绘制它们全部（使用较小的标记大小）
>>> plot_point(results['shgo'], color='r', marker='+')
>>> plot_point(results['shgo_sobol'], color='r', marker='x')
>>> for i in range(results['shgo_sobol'].xl.shape[0]):
... ax.plot(512 + results['shgo_sobol'].xl[i, 0],
... 512 + results['shgo_sobol'].xl[i, 1],
... 'ro', ms=2)
>>> ax.set_xlim([-4, 514*2])
>>> ax.set_ylim([-4, 514*2])
>>> plt.show()
import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize
def eggholder(x):
return (-(x[1] + 47) * np.sin(np.sqrt(abs(x[0]/2 + (x[1] + 47))))
-x[0] * np.sin(np.sqrt(abs(x[0] - (x[1] + 47)))))
bounds = [(-512, 512), (-512, 512)]
x = np.arange(-512, 513)
y = np.arange(-512, 513)
xgrid, ygrid = np.meshgrid(x, y)
xy = np.stack([xgrid, ygrid])
results = dict()
results['shgo'] = optimize.shgo(eggholder, bounds)
results['DA'] = optimize.dual_annealing(eggholder, bounds)
results['DE'] = optimize.differential_evolution(eggholder, bounds)
results['shgo_sobol'] = optimize.shgo(eggholder, bounds, n=256, iters=5,
sampling_method='sobol')
fig = plt.figure(figsize=(4.5, 4.5))
ax = fig.add_subplot(111)
im = ax.imshow(eggholder(xy), interpolation='bilinear', origin='lower',
cmap='gray')
ax.set_xlabel('x')
ax.set_ylabel('y')
def plot_point(res, marker='o', color=None):
ax.plot(512+res.x[0], 512+res.x[1], marker=marker, color=color, ms=10)
plot_point(results['DE'], color='c') # differential_evolution - cyan
plot_point(results['DA'], color='w') # dual_annealing. - white
# SHGO produces multiple minima, plot them all (with a smaller marker size)
plot_point(results['shgo'], color='r', marker='+')
plot_point(results['shgo_sobol'], color='r', marker='x')
for i in range(results['shgo_sobol'].xl.shape[0]):
ax.plot(512 + results['shgo_sobol'].xl[i, 0],
512 + results['shgo_sobol'].xl[i, 1],
'ro', ms=2)
ax.set_xlim([-4, 514*2])
ax.set_ylim([-4, 514*2])
fig.tight_layout()
plt.show()
- alt:
“这个X-Y图是一个热图，Z值用最低点表示为黑色，最高值表示为白色。图像类似于一个旋转了45度的棋盘，但经过了重度平滑。SHGO优化器在网格上的许多最小值处放置了一个红点。SHGO在右上角显示全局最小值为红色X。使用双重退火找到的局部最小值在左上角用白色圆圈标记。使用basinhopping找到的另一个局部最小值在顶部中心用黄色标记。代码将差分进化结果绘制为青色圆圈，但在图中不可见。乍一看，不清楚这些谷底中哪一个才是真正的全局最小值。”
- include-source:
0
最小二乘最小化（least_squares）#
SciPy能够解决鲁棒化的边界约束非线性最小二乘问题：
这里 \(f_i(\mathbf{x})\) 是从 \(\mathbb{R}^n\) 到 \(\mathbb{R}\) 的平滑函数，我们称之为残差。标量值函数 \(\rho(\cdot)\) 的目的是减少异常值残差的影响，并有助于解决方案的鲁棒性，我们称之为损失函数。线性损失函数给出了一个标准的最小二乘问题。此外，允许以某些 \(x_j\) 的下限和上限形式存在的约束。
所有特定于最小二乘最小化的方法都利用了一个称为雅可比矩阵的 \(m \times n\) 偏导数矩阵，定义为 \(J_{ij} = \partial f_i / \partial x_j\)。强烈建议通过解析方法计算此矩阵并将其传递给 least_squares，否则，它将通过有限差分法进行估计，这将花费大量时间。
额外的时间，并且在困难的情况下可能非常不准确。
函数 least_squares 可用于拟合函数
\(\varphi(t; \mathbf{x})\) 到经验数据 \(\{(t_i, y_i), i = 0, \ldots, m-1\}\)。
为此，只需预先计算残差为
\(f_i(\mathbf{x}) = w_i (\varphi(t_i; \mathbf{x}) - y_i)\)，其中 \(w_i\)
是分配给每个观测值的权重。
解决拟合问题的示例#
这里我们考虑一个酶促反应 [1]。定义了11个残差为
其中 \(y_i\) 是测量值，\(u_i\) 是自变量的值。未知参数向量为 \(\mathbf{x} = (x_0, x_1, x_2, x_3)^T\)。如前所述，建议以封闭形式计算雅可比矩阵：
我们将使用 [2] 中定义的“困难”起点。为了找到物理上有意义的解，避免潜在的除以零并确保收敛到全局最小值，我们施加约束 \(0 \leq x_j \leq 100, j = 0, 1, 2, 3\)。
下面的代码实现了 \(\mathbf{x}\) 的最小二乘估计，并最终绘制了原始数据和拟合的模型函数：
- alt:
“此代码绘制了一个X-Y时间序列图。序列从左下角的(0, 0)开始，并迅速上升至最大值0.2，然后趋于平稳。拟合模型以平滑的橙色轨迹显示，并且很好地拟合了数据。”
>>> from scipy.optimize import least_squares
>>> def model(x, u):
... return x[0] * (u ** 2 + x[1] * u) / (u ** 2 + x[2] * u + x[3])
>>> def fun(x, u, y):
... return model(x, u) - y
>>> def jac(x, u, y):
... J = np.empty((u.size, x.size))
... den = u ** 2 + x[2] * u + x[3]
... num = u ** 2 + x[1] * u
... J[:, 0] = num / den
... J[:, 1] = x[0] * u / den
... J[:, 2] = -x[0] * num * u / den ** 2
... J[:, 3] = -x[0] * num / den ** 2
... return J
>>> u = np.array([4.0, 2.0, 1.0, 5.0e-1, 2.5e-1, 1.67e-1, 1.25e-1, 1.0e-1,
... 8.33e-2, 7.14e-2, 6.25e-2])
>>> y = np.array([1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2,
... 4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2])
>>> x0 = np.array([2.5, 3.9, 4.15, 3.9])
>>> res = least_squares(fun, x0, jac=jac, bounds=(0, 100), args=(u, y), verbose=1)
# 可能会有所不同
`ftol` 终止条件已满足。
函数评估次数 130，初始成本 4.4383e+00，最终成本 1.5375e-04，一阶最优性 4.92e-08。
>>> res.x
array([ 0.19280596, 0.19130423, 0.12306063, 0.13607247])
>>> import matplotlib.pyplot as plt
>>> u_test = np.linspace(0, 5)
>>> y_test = model(res.x, u_test)
>>> plt.plot(u, y, 'o', markersize=4, label='数据')
>>> plt.plot(u_test, y_test, label='拟合模型')
>>> plt.xlabel("u")
>>> plt.ylabel("y")
>>> plt.legend(loc='lower right')
>>> plt.show()
进一步的示例#
以下三个交互式示例更详细地展示了 least_squares 的使用方法。
在 scipy 中进行大规模束调整 展示了
least_squares的大规模能力，以及如何高效地计算稀疏雅可比矩阵的有限差分近似。在 scipy 中进行稳健的非线性回归 展示了如何在非线性回归中使用稳健的损失函数处理异常值。
在 scipy 中求解离散边值问题 探讨了如何求解大型方程组，并使用边界条件来实现解的期望性质。
有关实现背后的数学算法的详细信息，请参阅 least_squares 的文档。
单变量函数最小化器 (minimize_scalar)#
通常情况下，只需要找到一个单变量函数（即输入为标量的函数）的最小值。在这些情况下，已经开发了其他可以更快工作的优化技术。这些技术可以通过 minimize_scalar 函数访问，该函数提供了几种算法。
无约束最小化 (method='brent')#
实际上，有两种方法可以用于最小化单变量函数：brent 和 golden，但 golden 仅用于学术目的，很少使用。这些方法可以通过 minimize_scalar 中的 method 参数分别选择。brent 方法使用 Brent 算法来定位最小值。理想情况下，应提供一个包含所需最小值的区间（bracket 参数）。
bracket 是一个三元组 \(\left( a, b, c \right)\)，使得 \(f
\left( a \right) > f \left( b \right) < f \left( c \right)\) 且 \(a <
b < c\)。如果未提供此三元组，则可以选择两个起始点，并使用简单的行进算法从这些点找到一个bracket。如果这两个起始点未提供，将使用 0 和 `1`（这可能不是您函数的选择，并可能导致返回意外的最小值）。
以下是一个示例：
>>> from scipy.optimize import minimize_scalar
>>> f = lambda x: (x - 2) * (x + 1)**2
>>> res = minimize_scalar(f, method='brent')
>>> print(res.x)
1.0
有界最小化（method='bounded'）#
通常，在最小化发生之前，可以对解空间施加约束。minimize_scalar 中的 bounded 方法是一个约束最小化的示例，它为标量函数提供了基本的区间约束。区间约束允许最小化仅在两个固定端点之间进行，这些端点使用强制性的 bounds 参数指定。
例如，要在 \(x=5\) 附近找到 \(J_{1}\left( x \right)\) 的最小值，可以使用区间 \(\left[ 4, 7 \right]\) 作为约束调用 minimize_scalar。结果是 \(x_{\textrm{min}}=5.3314\)：
>>> from scipy.special import j1
>>> res = minimize_scalar(j1, bounds=(4, 7), method='bounded')
>>> res.x
5.33144184241
自定义最小化器#
有时，使用自定义方法作为（多变量或单变量）最小化器可能很有用，例如在使用 minimize 的某些库包装器（如 basinhopping）时。
我们可以通过传递一个可调用对象（函数或实现 __call__ 方法的对象）作为 method 参数来实现这一点，而不是传递方法名称。 让我们考虑一个（坦白说相当虚拟的）需求，即使用一个简单的自定义多变量最小化方法，该方法仅在每个维度上独立地以固定步长搜索邻域:
>>> from scipy.optimize import OptimizeResult
>>> def custmin(fun, x0, args=(), maxfev=None, stepsize=0.1,
... maxiter=100, callback=None, **options):
... bestx = x0
... besty = fun(x0)
... funcalls = 1
... niter = 0
... improved = True
... stop = False
...
... while improved and not stop and niter < maxiter:
... improved = False
... niter += 1
... for dim in range(np.size(x0)):
... for s in [bestx[dim] - stepsize, bestx[dim] + stepsize]:
... testx = np.copy(bestx)
... testx[dim] = s
... testy = fun(testx, *args)
... funcalls += 1
... if testy < besty:
... besty = testy
... bestx = testx
... improved = True
... if callback is not None:
... callback(bestx)
... if maxfev is not None and funcalls >= maxfev:
... stop = True
... break
...
... return OptimizeResult(fun=besty, x=bestx, nit=niter,
... nfev=funcalls, success=(niter > 1))
>>> x0 = [1.35, 0.9, 0.8, 1.1, 1.2]
>>> res = minimize(rosen, x0, method=custmin, options=dict(stepsize=0.05))
>>> res.x
array([1., 1., 1., 1., 1.])
在单变量优化的情况下，这同样适用:
>>> def custmin(fun, bracket, args=(), maxfev=None, stepsize=0.1,
... maxiter=100, callback=None, **options):
... bestx = (bracket[1] + bracket[0]) / 2.0
... besty = fun(bestx)
... funcalls = 1
... niter = 0
...
... improved = True
... stop = False
...
... while improved and not stop and niter < maxiter:
... improved = False
... niter += 1
... for testx in [bestx - stepsize, bestx + stepsize]:
... testy = fun(testx, *args)
... funcalls += 1
... if testy < besty:
... besty = testy
... bestx = testx
... improved = True
... if callback is not None:
... callback(bestx)
... if maxfev is not None and funcalls >= maxfev:
... stop = True
... break
...
... return OptimizeResult(fun=besty, x=bestx, nit=niter,
... nfev=funcalls, success=(niter > 1))
>>> def f(x):
... return (x - 2)**2 * (x + 2)**2
>>> res = minimize_scalar(f, bracket=(-3.5, 0), method=custmin,
... options=dict(stepsize = 0.05))
>>> res.x
-2.0
根查找#
标量函数#
如果有一个单变量方程，可以尝试多种不同的根查找算法。这些算法大多需要一个区间端点，在该区间内预期存在根（因为函数在该区间内改变符号）。通常，brentq 是最佳选择，但其他方法在某些情况下或出于学术目的可能有用。当没有区间可用，但有一个或多个导数可用时，newton`（或 ``halley`, secant）可能适用。这种情况尤其适用于函数定义在复平面的某个子集上，而区间方法无法使用的情况。
固定点求解#
与寻找函数零点密切相关的问题是寻找函数的固定点。函数的固定点是函数求值返回该点的点。
点：\(g\left(x\right)=x.\) 显然，\(g\) 的不动点是 \(f\left(x\right)=g\left(x\right)-x\) 的根。
等价地，\(f\) 的根是 \(g\left(x\right)=f\left(x\right)+x\) 的不动点。例程 fixed_point 提供了一种简单的迭代方法，使用Aitkens序列加速来估计给定起始点的 \(g\) 的不动点。
方程组#
求解一组非线性方程的根可以通过 root 函数实现。有几种方法可供选择，其中包括 hybr``（默认）和 ``lm，它们分别使用Powell的混合方法和MINPACK的Levenberg-Marquardt方法。
以下示例考虑单变量超越方程
其根可以通过以下方式找到:
>>> import numpy as np
>>> from scipy.optimize import root
>>> def func(x):
... return x + 2 * np.cos(x)
>>> sol = root(func, 0.3)
>>> sol.x
array([-1.02986653])
>>> sol.fun
array([ -6.66133815e-16])
现在考虑一组非线性方程
我们定义目标函数，使其也返回雅可比矩阵，并通过将 jac 参数设置为 True 来指示这一点。此外，这里使用Levenberg-Marquardt求解器。
>>> def func2(x):
... f = [x[0] * np.cos(x[1]) - 4,
... x[1]*x[0] - x[1] - 5]
... df = np.array([[np.cos(x[1]), -x[0] * np.sin(x[1])],
... [x[1], x[0] - 1]])
... return f, df
>>> sol = root(func2, [1, 1], jac=True, method='lm')
>>> sol.x
array([ 6.50409711, 0.90841421])
大规模问题的根查找#
root 中的方法 hybr 和 lm 无法处理非常大的问题
变量数量（N），因为它们需要在每个牛顿步中计算和求逆一个稠密的 N x N 雅可比矩阵。当 N 增长时，这变得相当低效。
例如，考虑以下问题：我们需要在正方形 \([0,1]\times[0,1]\) 上求解以下积分微分方程：
边界条件为在上边缘 \(P(x,1) = 1\)，在正方形的其他边界上 \(P=0\)。这可以通过将连续函数 P 近似为其在网格上的值来实现，即 \(P_{n,m}\approx{}P(n h, m h)\)，其中网格间距 h 很小。然后可以近似导数和积分；例如 \(\partial_x^2 P(x,y)\approx{}(P(x+h,y) - 2 P(x,y) + P(x-h,y))/h^2\)。问题等价于找到某个函数 residual(P) 的根，其中 P 是一个长度为 \(N_x N_y\) 的向量。
现在，由于 \(N_x N_y\) 可能很大，root 中的方法 hybr 或 lm 将花费很长时间来解决这个问题。然而，可以使用大规模求解器之一来找到解决方案，例如 krylov、broyden2 或 anderson。这些方法使用所谓的非精确牛顿法，它不是精确计算雅可比矩阵，而是形成其近似值。
我们现在可以如下解决这个问题：
import numpy as np
from scipy.optimize import root
from numpy import cosh, zeros_like, mgrid, zeros
# 参数
nx, ny = 75, 75
hx, hy = 1./(nx-1), 1./(ny-1)
P_left, P_right = 0, 0
P_top, P_bottom = 1, 0
def residual(P):
d2x = zeros_like(P)
d2y = zeros_like(P)
d2x[1:-1] = (P[2:] - 2*P[1:-1] + P[:-2]) / hx/hx
d2x[0] = (P[1] - 2*P[0] + P_left)/hx/hx
d2x[-1] = (P_right - 2*P[-1] + P[-2])/hx/hx
d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
d2y[:,0] = (P[:,1] - 2*P[:,0] + P_bottom)/hy/hy
d2y[:,-1] = (P_top - 2*P[:,-1] + P[:,-2])/hy/hy
return d2x + d2y + 5*cosh(P).mean()**2
# 求解
guess = zeros((nx, ny), float)
sol = root(residual, guess, method='krylov', options={'disp': True})
#sol = root(residual, guess, method='broyden2', options={'disp': True, 'max_rank': 50})
#sol = root(residual, guess, method='anderson', options={'disp': True, 'M': 10})
print('Residual: %g' % abs(residual(sol.x)).max())
# 可视化
import matplotlib.pyplot as plt
x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
plt.pcolormesh(x, y, sol.x, shading='gouraud')
plt.colorbar()
plt.show()
仍然太慢？预处理。#
在寻找函数 \(f_i({\bf x}) = 0\) 的零点时，
i = 1, 2, …, N，krylov 求解器大部分时间都花在
求解雅可比矩阵的逆上，
如果你有一个近似于逆矩阵的近似矩阵 \(M\approx{}J^{-1}\)，你可以用它来*预处理* 线性求解问题。其思想是，与其求解 \(J{\bf s}={\bf y}\)，不如求解 \(MJ{\bf s}=M{\bf y}\)：由于矩阵 \(MJ\) 比 \(J\) 更接近单位矩阵， 因此方程应该更容易被 Krylov 方法处理。
矩阵 M 可以作为选项 options['jac_options']['inner_M']
传递给使用 krylov 方法的 root。它可以是一个
（稀疏）矩阵或 scipy.sparse.linalg.LinearOperator 实例。
对于上一节中的问题，我们注意到要 求解由两部分组成：第一部分是拉普拉斯算子的应用，\([\partial_x^2 + \partial_y^2] P\)，第二部分是积分。我们实际上可以很容易地计算出对应于拉普拉斯算子部分的雅可比矩阵：我们知道在一维情况下
因此，整个二维算子由以下矩阵表示
对应于积分的雅可比矩阵 \(J_2\) 更难计算，而且由于其所有元素都非零，将很难求逆。另一方面，\(J_1\) 是一个相对简单的矩阵，可以通过 scipy.sparse.linalg.splu 求逆（或者可以通过 scipy.sparse.linalg.spilu 近似求逆）。因此，我们满足于取 \(M\approx{}J_1^{-1}\) 并希望得到最好的结果。
在下面的示例中，我们使用预条件子 \(M=J_1^{-1}\)。
from scipy.optimize import root
from scipy.sparse import spdiags, kron
from scipy.sparse.linalg import spilu, LinearOperator
from numpy import cosh, zeros_like, mgrid, zeros, eye
# parameters
nx, ny = 75, 75
hx, hy = 1./(nx-1), 1./(ny-1)
P_left, P_right = 0, 0
P_top, P_bottom = 1, 0
def get_preconditioner():
"""Compute the preconditioner M"""
diags_x = zeros((3, nx))
diags_x[0,:] = 1/hx/hx
diags_x[1,:] = -2/hx/hx
diags_x[2,:] = 1/hx/hx
Lx = spdiags(diags_x, [-1,0,1], nx, nx)
diags_y = zeros((3, ny))
diags_y[0,:] = 1/hy/hy
diags_y[1,:] = -2/hy/hy
diags_y[2,:] = 1/hy/hy
Ly = spdiags(diags_y, [-1,0,1], ny, ny)
J1 = kron(Lx, eye(ny)) + kron(eye(nx), Ly)
# Now we have the matrix `J_1`. We need to find its inverse `M` --
# however, since an approximate inverse is enough, we can use
# the *incomplete LU* decomposition
J1_ilu = spilu(J1)
# This returns an object with a method .solve() that evaluates
# the corresponding matrix-vector product. We need to wrap it into
# a LinearOperator before it can be passed to the Krylov methods:
M = LinearOperator(shape=(nx*ny, nx*ny), matvec=J1_ilu.solve)
return M
def solve(preconditioning=True):
"""Compute the solution"""
count = [0]
def residual(P):
count[0] += 1
d2x = zeros_like(P)
d2y = zeros_like(P)
d2x[1:-1] = (P[2:] - 2*P[1:-1] + P[:-2])/hx/hx
d2x[0] = (P[1] - 2*P[0] + P_left)/hx/hx
d2x[-1] = (P_right - 2*P[-1] + P[-2])/hx/hx
d2y[:,1:-1] = (P[:,2:] - 2*P[:,1:-1] + P[:,:-2])/hy/hy
d2y[:,0] = (P[:,1] - 2*P[:,0] + P_bottom)/hy/hy
d2y[:,-1] = (P_top - 2*P[:,-1] + P[:,-2])/hy/hy
return d2x + d2y + 5*cosh(P).mean()**2
# preconditioner
if preconditioning:
M = get_preconditioner()
else:
M = None
# solve
guess = zeros((nx, ny), float)
sol = root(residual, guess, method='krylov',
options={'disp': True,
'jac_options': {'inner_M': M}})
print('Residual', abs(residual(sol.x)).max())
print('Evaluations', count[0])
return sol.x
def main():
sol = solve(preconditioning=True)
# visualize
import matplotlib.pyplot as plt
x, y = mgrid[0:1:(nx*1j), 0:1:(ny*1j)]
plt.clf()
plt.pcolor(x, y, sol)
plt.clim(0, 1)
plt.colorbar()
plt.show()
if __name__ == "__main__":
main()
运行结果如下，首先是未使用预条件子的情况:
0: |F(x)| = 803.614; step 1; tol 0.000257947
1: |F(x)| = 345.912; step 1; tol 0.166755
2: |F(x)| = 139.159; step 1; tol 0.145657
3: |F(x)| = 27.3682; step 1; tol 0.0348109
4: |F(x)| = 1.03303; step 1; tol 0.00128227
5: |F(x)| = 0.0406634; step 1; tol 0.00139451
6: |F(x)| = 0.00344341; step 1; tol 0.00645373
7: |F(x)| = 0.000153671; step 1; tol 0.00179246
8: |F(x)| = 6.7424e-06; step 1; tol 0.00173256
Residual 3.57078908664e-07
Evaluations 317
然后是使用预条件子的情况:
0: |F(x)| = 136.993; step 1; tol 7.49599e-06
1: |F(x)| = 4.80983; step 1; tol 0.00110945
2: |F(x)| = 0.195942; step 1; tol 0.00149362
3: |F(x)| = 0.000563597; step 1; tol 7.44604e-06
4: |F(x)| = 1.00698e-09; 步数 1; 容差 2.87308e-12
残差 9.29603061195e-11
评估次数 77
使用预条件子将 residual 函数的评估次数减少了 4 倍。对于残差计算成本高的问题，良好的预条件子至关重要——它甚至可以决定问题在实践中是否可解。
预条件子是一门艺术、科学和行业。在这里，我们幸运地选择了一个简单且效果不错的选项，但这个主题的深度远不止于此。
线性规划 (linprog)#
函数 linprog 可以在满足线性等式和不等式约束的条件下最小化线性目标函数。这类问题被称为线性规划。线性规划解决以下形式的问题：
其中 \(x\) 是决策变量的向量；\(c\), \(b_{ub}\), \(b_{eq}\), \(l\), 和 \(u\) 是向量；而 \(A_{ub}\) 和 \(A_{eq}\) 是矩阵。
在本教程中，我们将尝试使用 linprog 解决一个典型的线性规划问题。
线性规划示例#
考虑以下简单的线性规划问题：
我们需要一些数学操作来将目标问题转换为 linprog 接受的格式。
首先，考虑目标函数。
我们希望最大化目标函数，但 linprog 只能接受最小化问题。这可以通过将最大化 \(29x_1 + 45x_2\) 转换为最小化 \(-29x_1 -45x_2\) 来轻松解决。此外，\(x_3, x_4\) 未在目标函数中显示。这意味着与 \(x_3, x_4\) 对应的权重为零。因此，目标函数可以转换为：
如果我们定义决策变量向量 \(x = [x_1, x_2, x_3, x_4]^T\)，则在此问题中 linprog 的目标权重向量 \(c\) 应为
接下来，让我们考虑两个不等式约束。第一个是“小于”不等式，因此它已经是 linprog 接受的形式。第二个是“大于”不等式，因此我们需要将两边乘以 \(-1\) 以将其转换为“小于”不等式。显式显示零系数，我们有：
这些方程可以转换为矩阵形式：
其中
接下来，让我们考虑两个等式约束。显式显示零权重，这些是：
这些方程可以转换为矩阵形式：
其中
最后，让我们考虑对单个决策变量的单独不等式约束，这些约束被称为“箱式约束”或“简单边界”。这些约束可以通过 linprog 的 bounds 参数来应用。如 linprog 文档中所述，bounds 的默认值是 (0, None)，这意味着每个决策变量的下界为 0，上界为无穷大：所有决策变量都是非负的。我们的边界不同，因此我们需要将每个决策变量的上下界指定为一个元组，并将这些元组组合成一个列表。
最后，我们可以使用 linprog 来求解转换后的问题。
>>> import numpy as np
>>> from scipy.optimize import linprog
>>> c = np.array([-29.0, -45.0, 0.0, 0.0])
>>> A_ub = np.array([[1.0, -1.0, -3.0, 0.0],
... [-2.0, 3.0, 7.0, -3.0]])
>>> b_ub = np.array([5.0, -10.0])
>>> A_eq = np.array([[2.0, 8.0, 1.0, 0.0],
... [4.0, 4.0, 0.0, 1.0]])
>>> b_eq = np.array([60.0, 60.0])
>>> x0_bounds = (0, None)
>>> x1_bounds = (0, 5.0)
>>> x2_bounds = (-np.inf, 0.5) # +/- np.inf 可以用来代替 None
>>> x3_bounds = (-3.0, None)
>>> bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
>>> result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
>>> print(result.message)
The problem is infeasible. (HiGHS Status 8: model_status is Infeasible; primal_status is At lower/fixed bound)
结果表明我们的问题是不可行的，这意味着没有满足所有约束的解向量。这并不一定意味着我们做错了什么；有些问题确实是不可行的。
假设我们决定将 \(x_1\) 的边界约束过于严格，并可以放宽到 \(0 \leq x_1 \leq 6\)。在调整代码 x1_bounds = (0, 6) 以反映更改并再次执行后：
>>> x1_bounds = (0, 6)
>>> bounds = [x0_bounds, x1_bounds, x2_bounds, x3_bounds]
>>> result = linprog(c, A_ub=A_ub, b_ub=b_ub, A_eq=A_eq, b_eq=b_eq, bounds=bounds)
>>> print(result.message)
Optimization terminated successfully. (HiGHS Status 7: Optimal)
结果显示优化成功。
我们可以检查目标值 (result.fun) 是否与 \(c^Tx\) 相同：
>>> x = np.array(result.x)
>>> obj = result.fun
>>> print(c @ x)
-505.97435889013434 # 可能会有所不同
>>> print(obj)
-505.97435889013434 # 可能会有所不同
我们还可以检查所有约束是否在合理容差范围内满足：
>>> print(b_ub - (A_ub @ x).flatten()) # 这等同于 result.slack
[ 6.52747190e-10, -2.26730279e-09] # 可能会有所不同
>>> print(b_eq - (A_eq @ x).flatten()) # 这等同于 result.con
[ 9.78840831e-09, 1.04662945e-08]] # 可能会有所不同
>>> print([0 <= result.x[0], 0 <= result.x[1] <= 6.0, result.x[2] <= 0.5, -3.0 <= result.x[3]])
[True, True, True, True]
分配问题#
线性求和分配问题示例#
考虑为游泳混合接力队选拔学生的问题。 我们有一张表格显示了五名学生在每种游泳风格的时间：
我们需要为四种游泳方式各选择一名学生，以使总接力时间最小化。
这是一个典型的线性求和分配问题。我们可以使用 linear_sum_assignment 来解决它。
线性求和分配问题是其中最著名的组合优化问题之一。 给定一个“成本矩阵” \(C\)，问题是选择
每行中恰好一个元素
且不从任何列中选择超过一个元素
使得所选元素的总和最小化
换句话说，我们需要将每一行分配给一列，使得相应条目的总和最小化。
形式上，设 \(X\) 是一个布尔矩阵，其中 \(X[i,j] = 1\) 当且仅当行 \(i\) 被分配给列 \(j\)。 那么最优分配的成本为
第一步是定义成本矩阵。
在这个例子中，我们希望将每种游泳方式分配给一名学生。
linear_sum_assignment 能够将成本矩阵的每一行分配给一列。
因此，为了形成成本矩阵，需要将上表转置，使得行对应游泳方式，列对应学生：
>>> import numpy as np
>>> cost = np.array([[43.5, 45.5, 43.4, 46.5, 46.3],
... [47.1, 42.1, 39.1, 44.1, 47.8],
... [48.4, 49.6, 42.1, 44.5, 50.4],
... [38.2, 36.8, 43.2, 41.2, 37.2]])
我们可以使用 linear_sum_assignment 来解决分配问题：
>>> from scipy.optimize import linear_sum_assignment
>>> row_ind, col_ind = linear_sum_assignment(cost)
row_ind 和 col_ind 是成本矩阵的最优分配矩阵索引：
>>> row_ind
array([0, 1, 2, 3])
>>> col_ind
array([0, 2, 3, 1])
最优分配为：
>>> styles = np.array(["backstroke", "breaststroke", "butterfly", "freestyle"])[row_ind]
>>> students = np.array(["A", "B", "C", "D", "E"])[col_ind]
>>> dict(zip(styles, students))
{'backstroke': 'A', 'breaststroke': 'C', 'butterfly': 'D', 'freestyle': 'B'}
最优的总混合泳时间为：
>>> cost[row_ind, col_ind].sum()
163.89999999999998
注意，这个结果与每种泳姿的最小时间总和不同：
>>> np.min(cost, axis=1).sum()
161.39999999999998
因为学生“C”在“蛙泳”和“蝶泳”两种泳姿中都是最好的游泳者。 我们不能将学生“C”分配到两种泳姿中，所以我们将学生C分配到“蛙泳”泳姿， 并将D分配到“蝶泳”泳姿以最小化总时间。
参考文献
一些进一步阅读和相关的软件，如Newton-Krylov [KK]， PETSc [PP]，和PyAMG [AMG]：
D.A. Knoll 和 D.E. Keyes，“Jacobian-free Newton-Krylov methods”， J. Comp. Phys. 193, 357 (2004). DOI:10.1016/j.jcp.2003.08.010
PETSc https://www.mcs.anl.gov/petsc/ 及其 Python 绑定 https://bitbucket.org/petsc/petsc4py/
PyAMG（代数多重网格预处理器/求解器） pyamg/pyamg#issues
混合整数线性规划#
背包问题示例#
背包问题是众所周知的组合优化问题。 给定一组物品，每个物品都有一个大小和一个价值，问题是选择 那些物品以最大化总价值，同时总大小低于某个阈值。
形式上，设
\(x_i\) 是一个布尔变量，表示物品 \(i\) 是否包含在背包中，
\(n\) 是物品的总数，
\(v_i\) 是物品 \(i\) 的价值，
\(s_i\) 表示物品 \(i\) 的大小，
\(C\) 表示背包的容量。
那么问题可以表述为：
尽管目标函数和不平等约束在*决策变量* \(x_i\) 中是线性的，但这与典型的线性规划问题不同，因为决策变量只能取整数值。具体来说，我们的决策变量只能是 \(0\) 或 \(1\)，因此这被称为*二进制整数线性规划*（BILP）。此类问题属于更大的*混合整数线性规划*（MILPs）类别，我们可以使用 milp 来解决。
在我们的示例中，有8个物品可供选择，每个物品的大小和价值如下所示。
>>> import numpy as np
>>> from scipy import optimize
>>> sizes = np.array([21, 11, 15, 9, 34, 25, 41, 52])
>>> values = np.array([22, 12, 16, 10, 35, 26, 42, 53])
我们需要将八个决策变量约束为二进制。我们通过添加 Bounds 约束来确保它们位于 \(0\) 和 \(1\) 之间，并应用“整数性”约束以确保它们*要么*是 \(0\) *要么*是 \(1\)。
>>> bounds = optimize.Bounds(0, 1) # 0 <= x_i <= 1
>>> integrality = np.full_like(values, True) # x_i 是整数
背包容量约束使用 LinearConstraint 指定。
>>> capacity = 100
>>> constraints = optimize.LinearConstraint(A=sizes, lb=0, ub=capacity)
如果我们遵循线性代数的常规规则，输入 A 应该是一个二维矩阵，而上下界 lb 和 ub 应该是一维向量，但只要输入可以广播到一致的形状，LinearConstraint 是宽容的。
使用上面定义的变量，我们可以使用以下方法解决背包问题：
milp。注意，milp 最小化目标函数，但我们希望最大化总价值，因此我们将 c 设置为价值的负数。
>>> from scipy.optimize import milp
>>> res = milp(c=-values, constraints=constraints,
... integrality=integrality, bounds=bounds)
让我们检查一下结果：
>>> res.success
True
>>> res.x
array([1., 1., 0., 1., 1., 1., 0., 0.])
这意味着我们应该选择物品 1、2、4、5、6 以在尺寸约束下优化总价值。请注意，这与我们解决 *线性规划松弛*（没有整数约束）并尝试对决策变量进行四舍五入所得到的结果不同。
>>> from scipy.optimize import milp
>>> res = milp(c=-values, constraints=constraints,
... integrality=False, bounds=bounds)
>>> res.x
array([1. , 1. , 1. , 1. ,
0.55882353, 1. , 0. , 0. ])
如果我们对这个解进行向上取整，得到 array([1., 1., 1., 1., 1., 1., 0., 0.])，我们的背包将超过容量约束，而如果我们向下取整，得到 array([1., 1., 1., 1., 0., 1., 0., 0.])，我们将得到一个次优解。
有关更多 MILP 教程，请参阅 SciPy Cookbooks 上的 Jupyter 笔记本：