SciPy 1.12 中文文档（一）

原文：docs.scipy.org/doc/scipy-1.12.0/index.html

SciPy 用户指南

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/index.html

SciPy 是一组建立在NumPy之上的数学算法和便利函数。它通过提供高级命令和类来操作和可视化数据，显著增强了 Python 的功能。

子包

SciPy 按照不同的科学计算领域划分为多个子包。以下是这些子包的总结表：

Subpackage	描述
cluster	聚类算法
constants	物理和数学常数
fftpack	快速傅里叶变换例程
integrate	积分和常微分方程求解器
interpolate	插值和平滑样条
io	输入和输出
linalg	线性代数
ndimage	N 维图像处理
odr	正交距离回归
optimize	优化和寻根例程
signal	信号处理
sparse	稀疏矩阵及其相关例程
spatial	空间数据结构和算法
special	特殊函数
stats	统计分布和函数

SciPy 子包需要单独导入，例如：

>>> from scipy import linalg, optimize 

下面是按子包组织的完整用户指南。

用户指南

• 特殊函数 (scipy.special)

• 积分 (scipy.integrate)

• 优化 (scipy.optimize)

• 插值 (scipy.interpolate)

• 傅里叶变换 (scipy.fft)

• 信号处理 (scipy.signal)

• 线性代数 (scipy.linalg)

• 稀疏数组 (scipy.sparse)

• 使用 ARPACK 解决稀疏特征值问题

• 压缩稀疏图例程 (scipy.sparse.csgraph)

• 空间数据结构和算法 (scipy.spatial)

• 统计学 (scipy.stats)

• 多维图像处理 (scipy.ndimage)

• 文件 IO (scipy.io)

可执行教程

你还可以在这里找到使用MyST Markdown格式的教程。这些可以通过Jupytext扩展程序打开为 Jupyter 笔记本。

可执行教程

• 插值过渡指南

用户指南

特殊函数（scipy.special）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/special.html

scipy.special包的主要特点是定义了许多数学物理专用函数。可用函数包括阿尔谢尔、椭圆、贝塞尔、伽玛、贝塔、超几何、拋物线圆柱、马修、球形波、斯特鲁维和开尔文函数。还有一些低级别的统计函数，不适合一般用途，因为这些函数的易用接口由stats模块提供。这些函数大多数可以接受数组参数，并返回数组结果，遵循数值 Python 中其他数学函数的广播规则。许多函数还接受复数作为输入。要获取带有一行描述的可用函数的完整列表，请键入>>> help(special). 每个函数还有自己的文档，可通过帮助访问。如果找不到需要的函数，请考虑编写并贡献给该库。您可以使用 C、Fortran 或 Python 编写该函数。在库的源代码中查找这些函数的示例。

实阶贝塞尔函数（jv, jn_zeros)

贝塞尔函数是满足贝塞尔微分方程的解族，其阶数可以是实数或复数α：

[x² \frac{d² y}{dx²} + x \frac{dy}{dx} + (x² - \alpha²)y = 0]

在其他用途中，这些函数出现在波传播问题中，例如薄鼓面的振动模式。这里是一个固定在边缘的圆形鼓面的例子：

>>> from scipy import special
>>> import numpy as np
>>> def drumhead_height(n, k, distance, angle, t):
...    kth_zero = special.jn_zeros(n, k)[-1]
...    return np.cos(t) * np.cos(n*angle) * special.jn(n, distance*kth_zero)
>>> theta = np.r_[0:2*np.pi:50j]
>>> radius = np.r_[0:1:50j]
>>> x = np.array([r * np.cos(theta) for r in radius])
>>> y = np.array([r * np.sin(theta) for r in radius])
>>> z = np.array([drumhead_height(1, 1, r, theta, 0.5) for r in radius]) 

>>> import matplotlib.pyplot as plt
>>> fig = plt.figure()
>>> ax = fig.add_axes(rect=(0, 0.05, 0.95, 0.95), projection='3d')
>>> ax.plot_surface(x, y, z, rstride=1, cstride=1, cmap='RdBu_r', vmin=-0.5, vmax=0.5)
>>> ax.set_xlabel('X')
>>> ax.set_ylabel('Y')
>>> ax.set_xticks(np.arange(-1, 1.1, 0.5))
>>> ax.set_yticks(np.arange(-1, 1.1, 0.5))
>>> ax.set_zlabel('Z')
>>> plt.show() 

特殊函数的 Cython 绑定（scipy.special.cython_special)

SciPy 还为 special 中许多函数提供了标量化、类型化的 Cython 绑定。以下 Cython 代码提供了如何使用这些函数的简单示例：

cimport scipy.special.cython_special as csc

cdef:
    double x = 1
    double complex z = 1 + 1j
    double si, ci, rgam
    double complex cgam

rgam = csc.gamma(x)
print(rgam)
cgam = csc.gamma(z)
print(cgam)
csc.sici(x, &si, &ci)
print(si, ci) 

（参见Cython 文档以获取有关编译 Cython 的帮助。）在这个例子中，函数csc.gamma基本上像其 ufunc 对应物gamma一样工作，尽管它以 C 类型作为参数而不是 NumPy 数组。特别需要注意的是，该函数被重载以支持实数和复数参数；编译时会选择正确的变体。函数csc.sici与sici稍有不同；对于 ufunc，我们可以写成ai, bi = sici(x)，而在 Cython 版本中，多个返回值作为指针传递。可以将其类比为使用输出数组调用 ufunc：sici(x, out=(si, ci))。

使用 Cython 绑定有两个潜在的优势：

• 它们避免 Python 函数开销

• 它们不需要 Python 全局解释器锁（GIL）。

以下部分讨论如何利用这些优势潜在地加快您的代码，当然，首先应该对代码进行分析，确保付出额外的努力是值得的。

避免 Python 函数开销

对于 special 中的 ufuncs，通过向函数传递数组来避免 Python 函数开销，即向量化。通常，这种方法效果很好，但有时在循环内部调用标量输入的特殊函数更方便，例如在实现自己的 ufunc 时。在这种情况下，Python 函数开销可能会显著。考虑以下示例：

import scipy.special as sc
cimport scipy.special.cython_special as csc

def python_tight_loop():
    cdef:
        int n
        double x = 1

    for n in range(100):
        sc.jv(n, x)

def cython_tight_loop():
    cdef:
        int n
        double x = 1

    for n in range(100):
        csc.jv(n, x) 

在一台计算机上，python_tight_loop运行大约需要 131 微秒，而cython_tight_loop运行大约需要 18.2 微秒。显然，这个例子是刻意制造的：可以只调用special.jv(np.arange(100), 1)，就能像在cython_tight_loop中一样快速得到结果。关键是，如果 Python 函数开销在您的代码中变得显著，那么 Cython 绑定可能会有用。

释放 GIL

人们经常需要在许多点评估特殊函数，通常这些评估可以平凡地并行化。由于 Cython 绑定不需要 GIL，因此可以使用 Cython 的prange函数轻松地并行运行它们。例如，假设我们想计算亥姆霍兹方程的基本解：

[\Delta_x G(x, y) + k²G(x, y) = \delta(x - y),]

其中[k]是波数，而[δ]是狄拉克δ函数。已知在二维空间中，唯一的（辐射）解是

[G(x, y) = \frac{i}{4}H_0^{(1)}(k|x - y|),]

其中[H_0^{(1)}]是第一类汉克尔函数，即hankel1函数。以下示例展示了如何并行计算此函数：

from libc.math cimport fabs
cimport cython
from cython.parallel cimport prange

import numpy as np
import scipy.special as sc
cimport scipy.special.cython_special as csc

def serial_G(k, x, y):
    return 0.25j*sc.hankel1(0, k*np.abs(x - y))

@cython.boundscheck(False)
@cython.wraparound(False)
cdef void _parallel_G(double k, double[:,:] x, double[:,:] y,
                      double complex[:,:] out) nogil:
    cdef int i, j

    for i in prange(x.shape[0]):
        for j in range(y.shape[0]):
            out[i,j] = 0.25j*csc.hankel1(0, k*fabs(x[i,j] - y[i,j]))

def parallel_G(k, x, y):
    out = np.empty_like(x, dtype='complex128')
    _parallel_G(k, x, y, out)
    return out 

（如果需要帮助编译 Cython 中的并行代码，请参见这里。）如果上述 Cython 代码在名为 test.pyx 的文件中，那么我们可以编写一个非正式的基准测试，比较该函数的并行和串行版本：

import timeit

import numpy as np

from test import serial_G, parallel_G

def main():
    k = 1
    x, y = np.linspace(-100, 100, 1000), np.linspace(-100, 100, 1000)
    x, y = np.meshgrid(x, y)

    def serial():
        serial_G(k, x, y)

    def parallel():
        parallel_G(k, x, y)

    time_serial = timeit.timeit(serial, number=3)
    time_parallel = timeit.timeit(parallel, number=3)
    print("Serial method took {:.3} seconds".format(time_serial))
    print("Parallel method took {:.3} seconds".format(time_parallel))

if __name__ == "__main__":
    main() 

在一台四核计算机上，串行方法花费了 1.29 秒，而并行方法只花费了 0.29 秒。

不在 scipy.special 中的函数

有些函数未包含在 scipy.special 中，因为它们可以利用 NumPy 和 SciPy 中现有的函数直接实现。为了避免重复造轮子，本节提供了几个这样的函数的实现示例，希望能说明如何处理类似的函数。在所有示例中，NumPy 被导入为 np，而 special 被导入为 sc。

二进熵函数：

def binary_entropy(x):
    return -(sc.xlogy(x, x) + sc.xlog1py(1 - x, -x))/np.log(2) 

[0, 1] 上的矩形阶跃函数：

def step(x):
    return 0.5*(np.sign(x) + np.sign(1 - x)) 

可以使用平移和缩放来得到任意阶跃函数。

阶梯函数：

def ramp(x):
    return np.maximum(0, x) 

积分（scipy.integrate）

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/integrate.html

子包scipy.integrate 提供了几种积分技术，包括普通微分方程积分器。您可以通过帮助命令了解模块的概述：

>>> help(integrate)
 Methods for Integrating Functions given function object.

 quad          -- General purpose integration.
 dblquad       -- General purpose double integration.
 tplquad       -- General purpose triple integration.
 fixed_quad    -- Integrate func(x) using Gaussian quadrature of order n.
 quadrature    -- Integrate with given tolerance using Gaussian quadrature.
 romberg       -- Integrate func using Romberg integration.

 Methods for Integrating Functions given fixed samples.

 trapezoid            -- Use trapezoidal rule to compute integral.
 cumulative_trapezoid -- Use trapezoidal rule to cumulatively compute integral.
 simpson              -- Use Simpson's rule to compute integral from samples.
 romb                 -- Use Romberg Integration to compute integral from
 -- (2**k + 1) evenly-spaced samples.

 See the special module's orthogonal polynomials (special) for Gaussian
 quadrature roots and weights for other weighting factors and regions.

 Interface to numerical integrators of ODE systems.

 odeint        -- General integration of ordinary differential equations.
 ode           -- Integrate ODE using VODE and ZVODE routines. 

通用积分（quad)

函数quad 提供了一种计算单变量函数在两点之间积分的方法。这些点可以是(\pm\infty)（(\pm) inf），表示无限限制。例如，假设您希望计算贝塞尔函数 jv(2.5, x) 在区间 ([0, 4.5]) 上的积分。

[I=\int_{0}^{4.5}J_{2.5}\left(x\right), dx.]

这可以使用 quad 计算：

>>> import scipy.integrate as integrate
>>> import scipy.special as special
>>> result = integrate.quad(lambda x: special.jv(2.5,x), 0, 4.5)
>>> result
(1.1178179380783249, 7.8663172481899801e-09) 

>>> from numpy import sqrt, sin, cos, pi
>>> I = sqrt(2/pi)*(18.0/27*sqrt(2)*cos(4.5) - 4.0/27*sqrt(2)*sin(4.5) +
...                 sqrt(2*pi) * special.fresnel(3/sqrt(pi))[0])
>>> I
1.117817938088701 

>>> print(abs(result[0]-I))
1.03761443881e-11 

quad 的第一个参数是一个“可调用”的 Python 对象（即函数、方法或类实例）。请注意在此情况下使用 lambda- 函数作为参数。接下来的两个参数是积分的上下限。返回值是一个元组，第一个元素是积分估计值，第二个元素是绝对积分误差的估计值。请注意，在这种情况下，这个积分的真实值是

[I=\sqrt{\frac{2}{\pi}}\left(\frac{18}{27}\sqrt{2}\cos\left(4.5\right)-\frac{4}{27}\sqrt{2}\sin\left(4.5\right)+\sqrt{2\pi}\textrm{Si}\left(\frac{3}{\sqrt{\pi}}\right)\right),]

其中

[\textrm{Si}\left(x\right)=\int_{0}^{x}\sin\left(\frac{\pi}{2}t^{2}\right), dt.]

是 Fresnel 正弦积分。请注意，数值计算的积分结果比精确结果高出 (1.04\times10^{-11}) — 远低于报告的误差估计。

如果要积分的函数需要额外的参数，可以在 args 参数中提供。假设要计算以下积分：

[I(a,b)=\int_{0}^{1} ax²+b , dx.]

这个积分可以通过以下代码计算：

>>> from scipy.integrate import quad
>>> def integrand(x, a, b):
...     return a*x**2 + b
...
>>> a = 2
>>> b = 1
>>> I = quad(integrand, 0, 1, args=(a,b))
>>> I
(1.6666666666666667, 1.8503717077085944e-14) 

quad 还允许使用 (\pm) inf 作为参数之一进行无限输入。例如，假设要计算指数积分的数值值：

[E_{n}\left(x\right)=\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}, dt.]

是所需的（并且忘记了可以将这个积分计算为special.expn(n,x)的事实）。函数special.expn的功能可以通过基于quad例程定义新函数vec_expint来复制：

>>> from scipy.integrate import quad
>>> import numpy as np
>>> def integrand(t, n, x):
...     return np.exp(-x*t) / t**n
... 

>>> def expint(n, x):
...     return quad(integrand, 1, np.inf, args=(n, x))[0]
... 

>>> vec_expint = np.vectorize(expint) 

>>> vec_expint(3, np.arange(1.0, 4.0, 0.5))
array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049])
>>> import scipy.special as special
>>> special.expn(3, np.arange(1.0,4.0,0.5))
array([ 0.1097,  0.0567,  0.0301,  0.0163,  0.0089,  0.0049]) 

被积函数甚至可以使用quad参数（尽管误差界限可能会由于使用quad中的积分函数而低估误差）。在这种情况下，积分是

[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}, dt, dx=\frac{1}{n}.]

>>> result = quad(lambda x: expint(3, x), 0, np.inf)
>>> print(result)
(0.33333333324560266, 2.8548934485373678e-09) 

>>> I3 = 1.0/3.0
>>> print(I3)
0.333333333333 

>>> print(I3 - result[0])
8.77306560731e-11 

最后一个例子显示，可以使用重复调用quad来处理多重积分。

警告

数值积分算法在有限数量的点上采样被积函数。因此，它们不能保证对任意被积函数和积分限的准确结果（或准确性估计）。例如，考虑高斯积分：

>>> def gaussian(x):
...     return np.exp(-x**2)
>>> res = integrate.quad(gaussian, -np.inf, np.inf)
>>> res
(1.7724538509055159, 1.4202636756659625e-08)
>>> np.allclose(res[0], np.sqrt(np.pi))  # compare against theoretical result
True 

由于被积函数除了在原点附近几乎为零，我们预期大但有限的积分限会得到相同的结果。然而：

>>> integrate.quad(gaussian, -10000, 10000)
(1.975190562208035e-203, 0.0) 

这是因为在quad中实现的自适应积分例程虽然按设计工作，但没有注意到函数在如此大的有限区间内的小而重要部分。为了获得最佳结果，请考虑使用紧密环绕被积函数重要部分的积分限。

>>> integrate.quad(gaussian, -15, 15)
(1.772453850905516, 8.476526631214648e-11) 

必要时可以将具有几个重要区域的被积函数分成若干部分。

一般的多重积分（dblquad, tplquad, nquad)

双重和三重积分的机制已经封装到dblquad和tplquad函数中。这些函数分别接受要积分的函数及四个或六个参数。所有内积分的限制必须定义为函数。

下面展示了使用双重积分计算几个(I_{n})值的示例：

>>> from scipy.integrate import quad, dblquad
>>> def I(n):
...     return dblquad(lambda t, x: np.exp(-x*t)/t**n, 0, np.inf, lambda x: 1, lambda x: np.inf)
... 

>>> print(I(4))
(0.2500000000043577, 1.29830334693681e-08)
>>> print(I(3))
(0.33333333325010883, 1.3888461883425516e-08)
>>> print(I(2))
(0.4999999999985751, 1.3894083651858995e-08) 

对于非常数限制的积分的一个例子

[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y , dx, dy=\frac{1}{96}.]

可以使用下面的表达式计算这个积分（请注意使用非常数 lambda 函数作为内积分上限）：

>>> from scipy.integrate import dblquad
>>> area = dblquad(lambda x, y: x*y, 0, 0.5, lambda x: 0, lambda x: 1-2*x)
>>> area
(0.010416666666666668, 1.1564823173178715e-16) 

对于 n 重积分，scipy 提供了函数nquad。积分边界是一个可迭代对象：要么是常数边界的列表，要么是非常数积分边界的函数列表。积分的顺序（因此也是边界）是从最内层的积分到最外层的。

上述积分

[I_{n}=\int_{0}^{\infty}\int_{1}^{\infty}\frac{e^{-xt}}{t^{n}}, dt, dx=\frac{1}{n}]

可以计算为

>>> from scipy import integrate
>>> N = 5
>>> def f(t, x):
...    return np.exp(-x*t) / t**N
...
>>> integrate.nquad(f, [[1, np.inf],[0, np.inf]])
(0.20000000000002294, 1.2239614263187945e-08) 

注意f的参数顺序必须与积分边界的顺序匹配；即，对于(t)的内积分区间为([1, \infty])，对于(x)的外积分区间为([0, \infty])。

非常数积分边界可以以类似的方式处理；如上例所示。

[I=\int_{y=0}^{1/2}\int_{x=0}^{1-2y} x y , dx, dy=\frac{1}{96}.]

可通过以下方式进行评估

>>> from scipy import integrate
>>> def f(x, y):
...     return x*y
...
>>> def bounds_y():
...     return [0, 0.5]
...
>>> def bounds_x(y):
...     return [0, 1-2*y]
...
>>> integrate.nquad(f, [bounds_x, bounds_y])
(0.010416666666666668, 4.101620128472366e-16) 

这与之前的结果相同。

高斯积分

还提供了一些函数，以便在固定区间上执行简单的高斯积分。第一个是fixed_quad，执行固定阶数的高斯积分。第二个函数是quadrature，执行多阶高斯积分，直到积分估计的差异低于用户提供的某个容差。这些函数都使用了模块scipy.special.orthogonal，该模块可以计算多种正交多项式的根和积分权重（这些多项式本身作为特殊函数返回多项式类的实例，例如，special.legendre）。

罗姆伯格积分法

罗姆伯格方法[WPR]是另一种用于数值积分的方法。请参阅romberg的帮助函数以获取更多细节。

使用样本进行积分

如果样本等间距，并且可用样本数为 (2^{k}+1)（其中 (k) 是整数），那么可以使用 Romberg romb 积分来获得高精度的积分估计。Romberg 积分使用梯形规则在与二的幂相关的步长上，并对这些估计进行理查逊外推，以更高精度地近似积分。

在任意间隔样本的情况下，有两个函数trapezoid和simpson可用。它们分别使用牛顿-科特斯一阶和二阶公式进行积分。梯形规则将函数近似为相邻点之间的直线，而辛普森规则将函数在三个相邻点之间近似为抛物线。

对于样本数为奇数且等间距的情况，如果函数是三阶或更低阶的多项式，则辛普森规则是精确的。如果样本不是等间距的，则结果只有在函数是二阶或更低阶的多项式时才是精确的。

>>> import numpy as np
>>> def f1(x):
...    return x**2
...
>>> def f2(x):
...    return x**3
...
>>> x = np.array([1,3,4])
>>> y1 = f1(x)
>>> from scipy import integrate
>>> I1 = integrate.simpson(y1, x)
>>> print(I1)
21.0 

这恰好对应于

[\int_{1}^{4} x² , dx = 21,]

而积分第二个函数

>>> y2 = f2(x)
>>> I2 = integrate.simpson(y2, x)
>>> print(I2)
61.5 

不对应于

[\int_{1}^{4} x³ , dx = 63.75]

因为 f2 中的多项式阶数大于二阶。

使用低级回调函数进行更快的积分

如果用户希望减少集成时间，可以通过scipy.LowLevelCallable将 C 函数指针传递给quad、dblquad、tplquad或nquad，它将在 Python 中进行集成并返回结果。这里的性能提升来自两个因素。主要改进是函数本身的编译提供的更快的函数评估。此外，在quad中，通过消除 C 和 Python 之间的函数调用，我们还提供了加速。对于像正弦这样的简单函数，这种方法可能提供大约 2 倍的速度改进，但对于更复杂的函数，可能会产生更明显的改进（10 倍以上）。因此，此功能专为希望通过写一些 C 来显著减少计算时间的用户而设计。

例如，可以通过ctypes在几个简单的步骤中使用该方法：

1.) 在 C 中编写一个带有函数签名double f(int n, double *x, void *user_data)的积分函数，其中x是包含函数 f 评估点的数组，user_data是您想要提供的任意附加数据。

/* testlib.c */
double  f(int  n,  double  *x,  void  *user_data)  {
  double  c  =  *(double  *)user_data;
  return  c  +  x[0]  -  x[1]  *  x[2];  /* corresponds to c + x - y * z */
} 

2.) 现在将此文件编译为共享/动态库（快速搜索将帮助解决这个问题，因为它依赖于操作系统）。用户必须链接任何使用的数学库等。在 Linux 上，看起来像这样：

$ gcc -shared -fPIC -o testlib.so testlib.c 

输出库将被称为testlib.so，但它可能具有不同的文件扩展名。现在已经创建了一个库，可以使用ctypes加载到 Python 中。

3.) 使用ctypes将共享库加载到 Python 中，并设置restypes和argtypes - 这使得 SciPy 能够正确解释函数：

import os, ctypes
from scipy import integrate, LowLevelCallable

lib = ctypes.CDLL(os.path.abspath('testlib.so'))
lib.f.restype = ctypes.c_double
lib.f.argtypes = (ctypes.c_int, ctypes.POINTER(ctypes.c_double), ctypes.c_void_p)

c = ctypes.c_double(1.0)
user_data = ctypes.cast(ctypes.pointer(c), ctypes.c_void_p)

func = LowLevelCallable(lib.f, user_data) 

函数中最后的void *user_data是可选的，如果不需要，可以省略（在 C 函数和 ctypes argtypes 中都是如此）。注意，坐标被传递为双精度数组，而不是单独的参数。

4.) 现在像往常一样集成库函数，这里使用nquad：

>>> integrate.nquad(func, [[0, 10], [-10, 0], [-1, 1]])
(1200.0, 1.1102230246251565e-11) 

返回的 Python 元组在减少的时间内正常返回。此方法可以使用所有可选参数，包括指定奇点、无限界等。

普通微分方程 (solve_ivp)

对一组普通微分方程（ODEs）进行积分，并给出初始条件是另一个有用的例子。SciPy 中提供了函数 solve_ivp 用于积分第一阶向量微分方程：

[\frac{d\mathbf{y}}{dt}=\mathbf{f}\left(\mathbf{y},t\right),]

给定初始条件 (\mathbf{y}\left(0\right)=y_{0})，其中 (\mathbf{y}) 是长度为 (N) 的向量，(\mathbf{f}) 是从 (\mathcal{R}^{N}) 到 (\mathcal{R}^{N}) 的映射。通过将中间导数引入 (\mathbf{y}) 向量，任何高阶常微分方程总可以通过这种类型的微分方程来减少。

例如，假设要找到以下二阶微分方程的解：

[\frac{d^{2}w}{dz^{2}}-zw(z)=0]

初始条件为 (w\left(0\right)=\frac{1}{\sqrt[3]{3^{2}}\Gamma\left(\frac{2}{3}\right)}) 和 (\left.\frac{dw}{dz}\right|_{z=0}=-\frac{1}{\sqrt[3]{3}\Gamma\left(\frac{1}{3}\right)}.) 已知带有这些边界条件的解为艾里函数

[w=\textrm{Ai}\left(z\right),]

这提供了使用 special.airy 来检查积分器的方法。

首先，通过设定 (\mathbf{y}=\left[\frac{dw}{dz},w\right]) 和 (t=z) 将这个 ODE 转换为标准形式。因此，微分方程变为

[\begin{split}\frac{d\mathbf{y}}{dt}=\left[\begin{array}{c} ty_{1}\ y_{0}\end{array}\right]=\left[\begin{array}{cc} 0 & t\ 1 & 0\end{array}\right]\left[\begin{array}{c} y_{0}\ y_{1}\end{array}\right]=\left[\begin{array}{cc} 0 & t\ 1 & 0\end{array}\right]\mathbf{y}.\end{split}]

换句话说，

[\mathbf{f}\left(\mathbf{y},t\right)=\mathbf{A}\left(t\right)\mathbf{y}.]

作为一个有趣的提醒，如果 (\mathbf{A}\left(t\right)) 与 (\int_{0}^{t}\mathbf{A}\left(\tau\right), d\tau) 在矩阵乘法下交换，则此线性微分方程在使用矩阵指数的精确解时有解：

[\mathbf{y}\left(t\right)=\exp\left(\int_{0}^{t}\mathbf{A}\left(\tau\right)d\tau\right)\mathbf{y}\left(0\right),]

但在这种情况下，(\mathbf{A}\left(t\right)) 及其积分不对易。

可以使用函数solve_ivp来解决这个微分方程。它需要导数fprime，时间跨度*[t_start, t_end]和初始条件向量y0作为输入参数，并返回一个对象，其y*字段是连续解值的数组，作为列。因此，初始条件给出在第一个输出列中。

>>> from scipy.integrate import solve_ivp
>>> from scipy.special import gamma, airy
>>> y1_0 = +1 / 3**(2/3) / gamma(2/3)
>>> y0_0 = -1 / 3**(1/3) / gamma(1/3)
>>> y0 = [y0_0, y1_0]
>>> def func(t, y):
...     return [t*y[1],y[0]]
...
>>> t_span = [0, 4]
>>> sol1 = solve_ivp(func, t_span, y0)
>>> print("sol1.t: {}".format(sol1.t))
sol1.t:    [0\.         0.10097672 1.04643602 1.91060117 2.49872472 3.08684827
 3.62692846 4\.        ] 

正如可以看到的，如果未另有指定，solve_ivp会自动确定其时间步长。为了比较solve_ivp的解与airy函数，由solve_ivp创建的时间向量被传递给airy函数。

>>> print("sol1.y[1]: {}".format(sol1.y[1]))
sol1.y[1]: [0.35502805 0.328952   0.12801343 0.04008508 0.01601291 0.00623879
 0.00356316 0.00405982]
>>> print("airy(sol.t)[0]: {}".format(airy(sol1.t)[0]))
airy(sol.t)[0]: [0.35502805 0.328952   0.12804768 0.03995804 0.01575943 0.00562799
 0.00201689 0.00095156] 

具有其标准参数的solve_ivp的解显示出与 airy 函数的显著偏差。为了减小这种偏差，可以使用相对和绝对容差。

>>> rtol, atol = (1e-8, 1e-8)
>>> sol2 = solve_ivp(func, t_span, y0, rtol=rtol, atol=atol)
>>> print("sol2.y[1][::6]: {}".format(sol2.y[1][0::6]))
sol2.y[1][::6]: [0.35502805 0.19145234 0.06368989 0.0205917  0.00554734 0.00106409]
>>> print("airy(sol2.t)[0][::6]: {}".format(airy(sol2.t)[0][::6]))
airy(sol2.t)[0][::6]: [0.35502805 0.19145234 0.06368989 0.0205917  0.00554733 0.00106406] 

要为solve_ivp的解指定用户定义的时间点，solve_ivp提供了两种可能性，也可以互补使用。通过将t_eval选项传递给函数调用solve_ivp返回在其输出中这些时间点的解。

>>> import numpy as np
>>> t = np.linspace(0, 4, 100)
>>> sol3 = solve_ivp(func, t_span, y0, t_eval=t) 

如果函数的雅可比矩阵已知，则可以将其传递给solve_ivp以获得更好的结果。但请注意，默认的积分方法RK45不支持雅可比矩阵，因此必须选择另一种积分方法。支持雅可比矩阵的积分方法之一是例如以下示例中的Radau方法。

>>> def gradient(t, y):
...     return [[0,t], [1,0]]
>>> sol4 = solve_ivp(func, t_span, y0, method='Radau', jac=gradient) 

解决具有带状雅可比矩阵的系统

odeint可以告知雅可比矩阵为banded。对于已知是僵硬的大型微分方程系统，这可以显著提高性能。

例如，我们将使用线性方法解 1-D Gray-Scott 偏微分方程[MOL]。在区间(x \in [0, L])上，函数(u(x, t))和(v(x, t))的 Gray-Scott 方程为

[\begin{split}\begin{split} \frac{\partial u}{\partial t} = D_u \frac{\partial² u}{\partial x²} - uv² + f(1-u) \ \frac{\partial v}{\partial t} = D_v \frac{\partial² v}{\partial x²} + uv² - (f + k)v \ \end{split}\end{split}]

其中(D_u)和(D_v)分别是分量(u)和(v)的扩散系数，(f)和(k)是常数。（有关系统的更多信息，请参见groups.csail.mit.edu/mac/projects/amorphous/GrayScott/）

我们假设诺依曼（即“无通量”）边界条件：

[\frac{\partial u}{\partial x}(0,t) = 0, \quad \frac{\partial v}{\partial x}(0,t) = 0, \quad \frac{\partial u}{\partial x}(L,t) = 0, \quad \frac{\partial v}{\partial x}(L,t) = 0]

为了应用线性方法，我们通过定义均匀间隔的(N)个点的网格(\left{x_0, x_1, \ldots, x_{N-1}\right})来离散化(x)变量，其中(x_0 = 0)，(x_{N-1} = L)。我们定义(u_j(t) \equiv u(x_k, t))和(v_j(t) \equiv v(x_k, t))，并用有限差分替换(x)导数。即，

[\frac{\partial² u}{\partial x²}(x_j, t) \rightarrow \frac{u_{j-1}(t) - 2 u_{j}(t) + u_{j+1}(t)}{(\Delta x)²}]

然后我们得到一个由(2N)个常微分方程组成的系统：

(1)[\begin{split} \begin{split} \frac{du_j}{dt} = \frac{D_u}{(\Delta x)²} \left(u_{j-1} - 2 u_{j} + u_{j+1}\right) -u_jv_j² + f(1 - u_j) \ \frac{dv_j}{dt} = \frac{D_v}{(\Delta x)²} \left(v_{j-1} - 2 v_{j} + v_{j+1}\right) + u_jv_j² - (f + k)v_j \end{split}\end{split}]

为方便起见，已省略了((t))参数。

为了强制边界条件，我们引入“虚拟”点(x_{-1})和(x_N)，定义(u_{-1}(t) \equiv u_1(t))，(u_N(t) \equiv u_{N-2}(t))；(v_{-1}(t))和(v_N(t))类似地定义。

然后

(2)[\begin{split} \begin{split} \frac{du_0}{dt} = \frac{D_u}{(\Delta x)²} \left(2u_{1} - 2 u_{0}\right) -u_0v_0² + f(1 - u_0) \ \frac{dv_0}{dt} = \frac{D_v}{(\Delta x)²} \left(2v_{1} - 2 v_{0}\right) + u_0v_0² - (f + k)v_0 \end{split}\end{split}]

和

(3)[\begin{split} \begin{split} \frac{du_{N-1}}{dt} = \frac{D_u}{(\Delta x)²} \left(2u_{N-2} - 2 u_{N-1}\right) -u_{N-1}v_{N-1}² + f(1 - u_{N-1}) \ \frac{dv_{N-1}}{dt} = \frac{D_v}{(\Delta x)²} \left(2v_{N-2} - 2 v_{N-1}\right) + u_{N-1}v_{N-1}² - (f + k)v_{N-1} \end{split}\end{split}]

我们的完整的(2N)个常微分方程系统是(1)对于(k = 1, 2, \ldots, N-2)，以及(2)和(3)。

我们现在可以开始在代码中实现这个系统。我们必须将 ({u_k}) 和 ({v_k}) 合并成长度为 (2N) 的单一向量。两个明显的选择是 ({u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}}) 和 ({u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}})。数学上讲，选择没有影响，但是选择会影响odeint如何高效地解决系统。原因在于变量顺序如何影响雅可比矩阵非零元素的模式。

当变量按 ({u_0, u_1, \ldots, u_{N-1}, v_0, v_1, \ldots, v_{N-1}}) 排序时，雅可比矩阵的非零元素模式是

[\begin{split}\begin{smallmatrix} * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \ * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 \ 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 \ 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 \ 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 \ 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & * & 0 \ 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 & 0 & * \ * & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 \ 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 & 0 \ 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 & 0 \ 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & 0 \ 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 \ 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & 0 & * & * & * \ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & 0 & 0 & 0 & * & * & * \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & 0 & 0 & 0 & 0 & 0 \ \end{smallmatrix}\end{split}]

使用变量交错排列的雅可比模式为 ({u_0, v_0, u_1, v_1, \ldots, u_{N-1}, v_{N-1}}) 是

[\begin{split}\begin{smallmatrix} * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * & 0 & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * & * & 0 \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * & 0 & * \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & 0 & * & * \ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 & * & * & * \ \end{smallmatrix}\end{split}]

在这两种情况下，只有五条非平凡对角线，但当变量交错时，带宽要小得多。也就是说，主对角线和主对角线上下两侧的两条对角线是非零的。这很重要，因为odeint 的输入 mu 和 ml 是雅可比矩阵的上下带宽。当变量交错时，mu 和 ml 都是 2。当变量以 ({v_k}) 排列跟随 ({u_k}) 时，上下带宽是 (N)。

决定好后，我们可以编写实现微分方程系统的函数。

首先，我们定义系统的源项和反应项的函数：

def G(u, v, f, k):
    return f * (1 - u) - u*v**2

def H(u, v, f, k):
    return -(f + k) * v + u*v**2 

接下来，我们定义计算微分方程组右端的函数：

def grayscott1d(y, t, f, k, Du, Dv, dx):
  """
 Differential equations for the 1-D Gray-Scott equations.

 The ODEs are derived using the method of lines.
 """
    # The vectors u and v are interleaved in y.  We define
    # views of u and v by slicing y.
    u = y[::2]
    v = y[1::2]

    # dydt is the return value of this function.
    dydt = np.empty_like(y)

    # Just like u and v are views of the interleaved vectors
    # in y, dudt and dvdt are views of the interleaved output
    # vectors in dydt.
    dudt = dydt[::2]
    dvdt = dydt[1::2]

    # Compute du/dt and dv/dt.  The end points and the interior points
    # are handled separately.
    dudt[0]    = G(u[0],    v[0],    f, k) + Du * (-2.0*u[0] + 2.0*u[1]) / dx**2
    dudt[1:-1] = G(u[1:-1], v[1:-1], f, k) + Du * np.diff(u,2) / dx**2
    dudt[-1]   = G(u[-1],   v[-1],   f, k) + Du * (- 2.0*u[-1] + 2.0*u[-2]) / dx**2
    dvdt[0]    = H(u[0],    v[0],    f, k) + Dv * (-2.0*v[0] + 2.0*v[1]) / dx**2
    dvdt[1:-1] = H(u[1:-1], v[1:-1], f, k) + Dv * np.diff(v,2) / dx**2
    dvdt[-1]   = H(u[-1],   v[-1],   f, k) + Dv * (-2.0*v[-1] + 2.0*v[-2]) / dx**2

    return dydt 

我们不会实现一个计算雅可比矩阵的函数，但我们会告诉odeint 雅可比矩阵是带状的。这使得底层求解器（LSODA）可以避免计算已知为零的值。对于大型系统，这显著提高了性能，正如以下 ipython 会话中所示。

首先，我们定义所需的输入：

In [30]: rng = np.random.default_rng()

In [31]: y0 = rng.standard_normal(5000)

In [32]: t = np.linspace(0, 50, 11)

In [33]: f = 0.024

In [34]: k = 0.055

In [35]: Du = 0.01

In [36]: Dv = 0.005

In [37]: dx = 0.025 

不利用雅可比矩阵的带状结构来计算时间：

In [38]: %timeit sola = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx))
1 loop, best of 3: 25.2 s per loop 

现在设置 ml=2 和 mu=2，这样odeint 知道雅可比矩阵是带状的：

In [39]: %timeit solb = odeint(grayscott1d, y0, t, args=(f, k, Du, Dv, dx), ml=2, mu=2)
10 loops, best of 3: 191 ms per loop 

这样快了不少！

让我们确保它们计算出了相同的结果：

In [41]: np.allclose(sola, solb)
Out[41]: True 

参考资料

[WPR]

en.wikipedia.org/wiki/Romberg’s_method

[MOL]

en.wikipedia.org/wiki/Method_of_lines

优化 (scipy.optimize)

原文：docs.scipy.org/doc/scipy-1.12.0/tutorial/optimize.html

目录

• 优化 (scipy.optimize)

  • 多元标量函数的无约束最小化 (minimize)

    • 尼尔德-米德单纯形算法 (method='Nelder-Mead')

    • 布洛伊登-弗莱彻-戈尔德法-沙诺算法 (method='BFGS')

      • 避免冗余计算

    • 牛顿共轭梯度算法 (method='Newton-CG')

      • 完整海森矩阵示例：

      • 海森矩阵乘积示例：

    • 信赖区域牛顿共轭梯度算法 (method='trust-ncg')

      • 完整海森矩阵示例：

      • 海森矩阵乘积示例：

    • 信赖区域截断广义兰兹-共轭梯度算法 (method='trust-krylov')

      • 完整海森矩阵示例：

      • 海森矩阵乘积示例：

    • 信赖区域近乎精确算法 (method='trust-exact')

  • 多元标量函数的约束最小化 (minimize)

    • 信赖区域约束算法 (method='trust-constr')

      • 定义边界约束：

      • 定义线性约束：

      • 定义非线性约束：

      • 解决优化问题：

    • 顺序最小二乘编程 (SLSQP) 算法 (method='SLSQP')

  • 全局优化

  • 最小二乘最小化 (least_squares)

    • 解决拟合问题的示例

    • 更多示例

  • 一元函数最小化器 (minimize_scalar)

    • 无约束最小化 (method='brent')

    • 有界最小化 (method='bounded')

  • 自定义最小化器

  • 根查找

    • 标量函数

    • 固定点求解

    • 方程组

    • 大问题的根查找

    • 仍然太慢？预调节。

  • 线性规划 (linprog)

    • 线性规划示例

  • 分配问题

    • 线性和分配问题示例

  • 混合整数线性规划

    • 背包问题示例

scipy.optimize 包提供了几种常用的优化算法。详细清单可参阅：scipy.optimize （也可通过help(scipy.optimize)查找）。

无约束多变量标量函数最小化 (minimize)

minimize 函数为多变量标量函数的无约束和约束最小化算法提供了通用接口，在scipy.optimize中使用。为了演示最小化函数，考虑如下问题，即最小化具有(N)个变量的 Rosenbrock 函数：

[f\left(\mathbf{x}\right)=\sum_{i=1}^{N-1}100\left(x_{i+1}-x_{i}^{2}\right)^{2}+\left(1-x_{i}\right)^{2}.]

此函数的最小值为 0，当(x_{i}=1.)时达到。

注意，Rosenbrock 函数及其导数包含在scipy.optimize中。以下部分中展示的实现示例说明了如何定义一个目标函数以及其雅可比和 Hessian 函数。scipy.optimize中的目标函数期望其第一个参数是 numpy 数组，用于优化，并且必须返回一个浮点值。确切的调用签名必须是f(x, *args)，其中x表示 numpy 数组，args是传递给目标函数的额外参数的元组。

Nelder-Mead Simplex 算法 (method='Nelder-Mead')

在下面的示例中，使用Nelder-Mead simplex 算法（通过method参数选择）调用了minimize例程：

>>> 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})
Optimization terminated successfully.
 Current function value: 0.000000
 Iterations: 339
 Function evaluations: 571 

>>> print(res.x)
[1\. 1\. 1\. 1\. 1.] 

Simplex 算法可能是最简单的最小化相当良好函数的方法。它只需要函数评估，并且对于简单的最小化问题是一个很好的选择。然而，由于它不使用任何梯度评估，可能需要更长时间找到最小值。

另一个仅需函数调用来找到最小值的优化算法是Powell方法，通过在minimize中设置method='powell'来选择。

为了演示如何向目标函数提供额外参数，让我们最小化 Rosenbrock 函数，其中包括一个额外的缩放因子a和一个偏移量b：

[f\left(\mathbf{x}, a, b\right)=\sum_{i=1}^{N-1}a\left(x_{i+1}-x_{i}^{2}\right)^{2}+\left(1-x_{i}\right)^{2} + b.]

再次使用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})
Optimization terminated successfully.
 Current function value: 1.000000
 Iterations: 319 # may vary
 Function evaluations: 525 # may vary 

>>> print(res.x)
[1\.         1\.         1\.         1\.         0.99999999] 

作为使用minimize的args参数的替代方法，只需将目标函数包装在一个新函数中，该函数仅接受x作为参数即可。当需要将额外参数作为关键字参数传递给目标函数时，这种方法也很有用。

>>> def rosen_with_args(x, a, *, b):  # b is a keyword-only argument
...     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.)  # pass in `a` and `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 函数的梯度是以下向量：

\begin{eqnarray*} \frac{\partial f}{\partial x_{j}} & = & \sum_{i=1}^{N}200\left(x_{i}-x_{i-1}^{2}\right)\left(\delta_{i,j}-2x_{i-1}\delta_{i-1,j}\right)-2\left(1-x_{i-1}\right)\delta_{i-1,j}.\ & = & 200\left(x_{j}-x_{j-1}^{2}\right)-400x_{j}\left(x_{j+1}-x_{j}^{2}\right)-2\left(1-x_{j}\right).\end{eqnarray*}

这个表达式对内部导数是有效的。特殊情况包括

\begin{eqnarray*} \frac{\partial f}{\partial x_{0}} & = & -400x_{0}\left(x_{1}-x_{0}^{2}\right)-2\left(1-x_{0}\right),\ \frac{\partial f}{\partial x_{N-1}} & = & 200\left(x_{N-1}-x_{N-2}^{2}\right).\end{eqnarray*}

用以下代码段构建一个计算此梯度的 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 

这些梯度信息在minimize函数中通过jac参数指定，如下所示。

>>> res = minimize(rosen, x0, method='BFGS', jac=rosen_der,
...                options={'disp': True})
Optimization terminated successfully.
 Current function value: 0.000000
 Iterations: 25                     # may vary
 Function evaluations: 30
 Gradient evaluations: 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):
...     # this function is computationally expensive!
...     expensive.count += 1  # let's keep track of how many times it runs
...     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 次：在目标函数中六次，在梯度中六次。减少冗余计算的一种方法是创建一个返回目标函数和梯度的单一函数。

>>> def f_and_df(x):
...     expensive_value = expensive(x[0])
...     return (-expensive_value**2,  # objective function
...             -2*expensive_value*dexpensive(x[0]))  # gradient
>>>
>>> expensive.count = 0  # reset the counter
>>> 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  # reset the counter
>>> expensive = lru_cache(expensive)
>>> res = minimize(f, 0.5, jac=df)
>>> res.fun
-0.9999999999999174
>>> expensive.count
6 

Newton-Conjugate-Gradient 算法（method='Newton-CG'）

Newton-Conjugate Gradient 算法是改进的 Newton 方法，使用共轭梯度算法（近似）求解局部 Hessian 矩阵的逆[NW]。Newton 方法基于将函数局部拟合为二次型：

[f\left(\mathbf{x}\right)\approx f\left(\mathbf{x}{0}\right)+\nabla f\left(\mathbf{x}{0}\right)\cdot\left(\mathbf{x}-\mathbf{x}{0}\right)+\frac{1}{2}\left(\mathbf{x}-\mathbf{x}{0}\right)^{T}\mathbf{H}\left(\mathbf{x}{0}\right)\left(\mathbf{x}-\mathbf{x}{0}\right).]

where (\mathbf{H}\left(\mathbf{x}_{0}\right)) is a matrix of second-derivatives (the Hessian). If the Hessian is positive definite then the local minimum of this function can be found by setting the gradient of the quadratic form to zero, resulting in

[\mathbf{x}{\textrm{opt}}=\mathbf{x}{0}-\mathbf{H}^{-1}\nabla f.]

The inverse of the Hessian is evaluated using the conjugate-gradient method. An example of employing this method to minimizing the Rosenbrock function is given below. To take full advantage of the Newton-CG method, a function which computes the Hessian must be provided. The Hessian matrix itself does not need to be constructed, only a vector which is the product of the Hessian with an arbitrary vector needs to be available to the minimization routine. As a result, the user can provide either a function to compute the Hessian matrix, or a function to compute the product of the Hessian with an arbitrary vector.

Full Hessian example:

The Hessian of the Rosenbrock function is

\begin{eqnarray*} H_{ij}=\frac{\partial^{2}f}{\partial x_{i}\partial x_{j}} & = & 200\left(\delta_{i,j}-2x_{i-1}\delta_{i-1,j}\right)-400x_{i}\left(\delta_{i+1,j}-2x_{i}\delta_{i,j}\right)-400\delta_{i,j}\left(x_{i+1}-x_{i}^{2}\right)+2\delta_{i,j},\ & = & \left(202+1200x_{i}^{2}-400x_{i+1}\right)\delta_{i,j}-400x_{i}\delta_{i+1,j}-400x_{i-1}\delta_{i-1,j},\end{eqnarray*}

if (i,j\in\left[1,N-2\right]) with (i,j\in\left[0,N-1\right]) defining the (N\times N) matrix. Other non-zero entries of the matrix are

\begin{eqnarray*} \frac{\partial^{2}f}{\partial x_{0}^{2}} & = & 1200x_{0}^{2}-400x_{1}+2,\ \frac{\partial^{2}f}{\partial x_{0}\partial x_{1}}=\frac{\partial^{2}f}{\partial x_{1}\partial x_{0}} & = & -400x_{0},\ \frac{\partial^{2}f}{\partial x_{N-1}\partial x_{N-2}}=\frac{\partial^{2}f}{\partial x_{N-2}\partial x_{N-1}} & = & -400x_{N-2},\ \frac{\partial^{2}f}{\partial x_{N-1}^{2}} & = & 200.\end{eqnarray*}

For example, the Hessian when (N=5) is

[\begin{split}\mathbf{H}=\begin{bmatrix} 1200x_{0}^{2}+2\mkern-2em\&1200x_{1}^{2}+202\mkern-2em\&&1200x_{1}^{2}+202\mkern-2em\&&&1200x_{3}^{2}+202\mkern-1em\&&&&200\end{bmatrix}-400\begin{bmatrix} x_1 & x_0 \ x_0 & x_2 & x_1 \ & x_1 & x_3 & x_2\ & & x_2 & x_4 & x_3 \ & & & x_3 & 0\end{bmatrix}.\end{split}]

The code which computes this Hessian along with the code to minimize the function using Newton-CG method is shown in the following example:

>>> 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})
Optimization terminated successfully.
 Current function value: 0.000000
 Iterations: 19                       # may vary
 Function evaluations: 22
 Gradient evaluations: 19
 Hessian evaluations: 19
>>> res.x
array([1.,  1.,  1.,  1.,  1.]) 

Hessian product example:

对于更大的最小化问题，存储整个 Hessian 矩阵可能会消耗大量时间和内存。Newton-CG 算法仅需 Hessian 与任意向量的乘积。因此，用户可以提供计算该乘积的代码，而不是完整的 Hessian 矩阵，方法是通过提供一个 hess 函数，该函数将最小化向量作为第一个参数，任意向量作为第二个参数（以及传递给要最小化函数的额外参数）。如果可能，使用具有 Hessian 乘积选项的 Newton-CG 可能是最快的最小化函数的方式。

在这种情况下，用 Rosenbrock Hessian 乘以任意向量的乘积并不难计算。如果 (\mathbf{p}) 是任意向量，则 (\mathbf{H}\left(\mathbf{x}\right)\mathbf{p}) 的元素为：

[\begin{split}\mathbf{H}\left(\mathbf{x}\right)\mathbf{p}=\begin{bmatrix} \left(1200x_{0}^{2}-400x_{1}+2\right)p_{0}-400x_{0}p_{1}\ \vdots\ -400x_{i-1}p_{i-1}+\left(202+1200x_{i}^{2}-400x_{i+1}\right)p_{i}-400x_{i}p_{i+1}\ \vdots\ -400x_{N-2}p_{N-2}+200p_{N-1}\end{bmatrix}.\end{split}]

使用 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})
Optimization terminated successfully.
 Current function value: 0.000000
 Iterations: 20                    # may vary
 Function evaluations: 23
 Gradient evaluations: 20
 Hessian evaluations: 44
>>> res.x
array([1., 1., 1., 1., 1.]) 

根据 [NW] 第 170 页，当 Hessian 矩阵条件差时，Newton-CG 算法可能效率低下，因为在这些情况下方法提供的搜索方向质量较差。作者表示，trust-ncg 方法更有效地处理这种问题情况，并将在下文中详细描述。

Trust-Region Newton-Conjugate-Gradient Algorithm (method='trust-ncg')

Newton-CG 方法是一种线搜索方法：它找到一个搜索方向，最小化函数的二次近似，然后使用一种线搜索算法在该方向上找到（几乎）最优的步长。另一种方法是先固定步长限制 (\Delta)，然后通过解以下二次子问题找到给定信赖半径内的最优步 (\mathbf{p})：

\begin{eqnarray*} \min_{\mathbf{p}} f\left(\mathbf{x}{k}\right)+\nabla f\left(\mathbf{x}{k}\right)\cdot\mathbf{p}+\frac{1}{2}\mathbf{p}^{T}\mathbf{H}\left(\mathbf{x}_{k}\right)\mathbf{p};&\ \text{subject to: } |\mathbf{p}|\le \Delta.& \end{eqnarray*}

解被更新为 (\mathbf{x}{k+1} = \mathbf{x}{k} + \mathbf{p})，并且信赖半径 (\Delta) 根据二次模型与实际函数的一致程度进行调整。这类方法被称为信赖域方法。trust-ncg 算法是一种信赖域方法，它使用共轭梯度算法来解决信赖域子问题 [NW]。

Full Hessian example:

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

Hessian product example:

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

Trust-Region Truncated Generalized Lanczos / Conjugate Gradient Algorithm (method='trust-krylov')

与trust-ncg方法类似，trust-krylov方法适用于大规模问题，因为它只使用海森矩阵作为线性操作符，通过矩阵-向量乘积来解决二次子问题。它比trust-ncg方法更准确地解决了二次子问题。

\begin{eqnarray*} \min_{\mathbf{p}} f\left(\mathbf{x}{k}\right)+\nabla f\left(\mathbf{x}{k}\right)\cdot\mathbf{p}+\frac{1}{2}\mathbf{p}^{T}\mathbf{H}\left(\mathbf{x}_{k}\right)\mathbf{p};&\ \text{subject to: } |\mathbf{p}|\le \Delta.& \end{eqnarray*}

该方法使用[TRLIB]实现了[GLTR]方法，精确求解了限制在截断 Krylov 子空间中的信任域子问题。对于不定问题，使用该方法通常更好，因为它减少了非线性迭代的次数，而在每个子问题求解中增加了少量的矩阵-向量乘积。

Full Hessian example:

>>> 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                    # may vary
 Function evaluations: 20
 Gradient evaluations: 20
 Hessian evaluations: 18
>>> res.x
array([1., 1., 1., 1., 1.]) 

Hessian product example:

>>> 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                    # may vary
 Function evaluations: 20
 Gradient evaluations: 20
 Hessian evaluations: 0
>>> res.x
array([1., 1., 1., 1., 1.]) 

[TRLIB]

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

[GLTR]

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

Trust-Region Nearly Exact Algorithm (method='trust-exact')

所有方法Newton-CG，trust-ncg和trust-krylov都适用于处理大规模问题（具有数千个变量的问题）。这是因为共轭梯度算法通过迭代近似地解决信任域子问题（或反转海森矩阵），而不需要显式的海森矩阵分解。由于只需海森矩阵与任意向量的乘积，因此该算法特别适用于处理稀疏海森矩阵，从而实现低存储需求和在这些稀疏问题中显著的时间节省。

对于中等大小的问题，其中海森矩阵的存储和分解成本并不关键，通过几乎精确地求解信任域子问题，可以在较少的迭代次数内获得解决方案。为了实现这一点，对每个二次子问题进行了某些非线性方程的迭代求解 [CGT]。这种解决方案通常需要对海森矩阵进行 3 到 4 次乔列斯基分解。因此，该方法收敛的迭代次数较少，并且比其他实现的信任域方法需要更少的目标函数评估。该算法不支持海森乘积选项。以下是使用 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                    # may vary
 Function evaluations: 14
 Gradient evaluations: 13
 Hessian evaluations: 14
>>> res.x
array([1., 1., 1., 1., 1.]) 

[NW] (1,2,3)

J. Nocedal, S.J. Wright：“Numerical optimization.” 第 2 版. Springer Science (2006).

[CGT]

康恩, A. R., 戈尔德, N. I., & 琼特, P. L. “信赖域方法”. Siam. (2000). pp. 169-200.

多变量标量函数的约束最小化 (minimize)

minimize 函数提供了约束最小化的算法，即 'trust-constr' 、 'SLSQP' 和 'COBYLA' 。它们要求约束使用稍微不同的结构定义。 'trust-constr' 方法要求约束以 LinearConstraint 和 NonlinearConstraint 对象序列的形式给出。另一方面， 'SLSQP' 和 'COBYLA' 方法要求约束以字典序列的形式给出，其中包括 type 、 fun 和 jac 键。

例如，让我们考虑对 Rosenbrock 函数的约束最小化：

\begin{eqnarray*} \min_{x_0, x_1} & ~~100\left(x_{1}-x_{0}^{2}\right)^{2}+\left(1-x_{0}\right)^{2} &\ \text{subject to: } & x_0 + 2 x_1 \leq 1 & \ & x_0² + x_1 \leq 1 & \ & x_0² - x_1 \leq 1 & \ & 2 x_0 + x_1 = 1 & \ & 0 \leq x_0 \leq 1 & \ & -0.5 \leq x_1 \leq 2.0. & \end{eqnarray*}

该优化问题有唯一解 ([x_0, x_1] = [0.4149,~ 0.1701])，其中仅第一和第四个约束是活动约束。

信赖域约束算法 (method='trust-constr')

信赖域约束方法处理以下形式的约束最小化问题：

\begin{eqnarray*} \min_x & f(x) & \ \text{subject to: } & ~~~ c^l \leq c(x) \leq c^u, &\ & x^l \leq x \leq x^u. & \end{eqnarray*}

当 (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) 可以用线性约束标准格式写成：

\begin{equation*} \begin{bmatrix}-\infty \1\end{bmatrix} \leq \begin{bmatrix} 1& 2 \ 2& 1\end{bmatrix} \begin{bmatrix} x_0 \x_1\end{bmatrix} \leq \begin{bmatrix} 1 \ 1\end{bmatrix},\end{equation*}

并使用 LinearConstraint 对象定义。

>>> from scipy.optimize import LinearConstraint
>>> linear_constraint = LinearConstraint([[1, 2], [2, 1]], [-np.inf, 1], [1, 1]) 

定义非线性约束:

非线性约束：

\begin{equation*} c(x) = \begin{bmatrix} x_0² + x_1 \ x_0² - x_1\end{bmatrix} \leq \begin{bmatrix} 1 \ 1\end{bmatrix}, \end{equation*}

其中 Jacobian 矩阵为：

\begin{equation*} J(x) = \begin{bmatrix} 2x_0 & 1 \ 2x_0 & -1\end{bmatrix},\end{equation*}

和 Hessian 的线性组合：

\begin{equation*} H(x, v) = \sum_{i=0}¹ v_i \nabla² c_i(x) = v_0\begin{bmatrix} 2 & 0 \ 0 & 0\end{bmatrix} + v_1\begin{bmatrix} 2 & 0 \ 0 & 0\end{bmatrix}, \end{equation*}

使用 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)
# may vary
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.016 s.
>>> print(res.x)
[0.41494531 0.17010937] 

当需要时，可以使用 LinearOperator 对象定义目标函数的 Hessian，

>>> 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)
# may vary
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.018 s.
>>> print(res.x)
[0.41494531 0.17010937] 

或者通过参数 hessp 实现 Hessian-向量乘积。

>>> res = minimize(rosen, x0, method='trust-constr', jac=rosen_der, hessp=rosen_hess_p,
...                constraints=[linear_constraint, nonlinear_constraint],
...                options={'verbose': 1}, bounds=bounds)
# may vary
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 8, CG iterations: 7, optimality: 2.99e-09, constraint violation: 1.11e-16, execution time: 0.018 s.
>>> 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)
# may vary
`gtol` termination condition is satisfied.
Number of iterations: 12, function evaluations: 24, CG iterations: 7, optimality: 4.48e-09, constraint violation: 0.00e+00, execution time: 0.016 s.
>>> print(res.x)
[0.41494531 0.17010937] 

[旅行]

Byrd, Richard H., Mary E. Hribar, and Jorge Nocedal. 1999. An interior point algorithm for large-scale nonlinear programming. SIAM Journal on Optimization 9.4: 877-900.

[EQSQP]

Lalee, Marucha, Jorge Nocedal, and Todd Plantega. 1998. On the implementation of an algorithm for large-scale equality constrained optimization. SIAM Journal on Optimization 8.3: 682-706.

顺序最小二乘规划（SLSQP）算法（method='SLSQP'）

SLSQP 方法处理形如以下约束极小化问题：

\begin{eqnarray*} \min_x & f(x) \ \text{subject to: } & c_j(x) = 0 , &j \in \mathcal{E}\ & c_j(x) \geq 0 , &j \in \mathcal{I}\ & \text{lb}_i \leq x_i \leq \text{ub}_i , &i = 1,...,N. \end{eqnarray*}

其中(\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)
# may vary
Optimization terminated successfully.    (Exit mode 0)
 Current function value: 0.342717574857755
 Iterations: 5
 Function evaluations: 6
 Gradient evaluations: 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  # may vary
 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  # may vary
 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 - 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])
>>> plt.show() 

最小二乘最小化（least_squares)

SciPy 能够解决鲁棒的有界非线性最小二乘问题：

\begin{align} &\min_\mathbf{x} \frac{1}{2} \sum_{i = 1}^m \rho\left(f_i(\mathbf{x})²\right) \ &\text{subject to }\mathbf{lb} \leq \mathbf{x} \leq \mathbf{ub} \end{align}

这里 ( 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 个残差为

[ f_i(x) = \frac{x_0 (u_i² + u_i x_1)}{u_i² + u_i x_2 + x_3} - y_i, \quad i = 0, \ldots, 10, ]

其中(y_i)是测量值，(u_i)是自变量值。未知参数向量为(\mathbf{x} = (x_0, x_1, x_2, x_3)^T)。如前所述，建议以闭合形式计算雅可比矩阵：

\begin{align} &J_{i0} = \frac{\partial f_i}{\partial x_0} = \frac{u_i² + u_i x_1}{u_i² + u_i x_2 + x_3} \ &J_{i1} = \frac{\partial f_i}{\partial x_1} = \frac{u_i x_0}{u_i² + u_i x_2 + x_3} \ &J_{i2} = \frac{\partial f_i}{\partial x_2} = -\frac{x_0 (u_i² + u_i x_1) u_i}{(u_i² + u_i x_2 + x_3)²} \ &J_{i3} = \frac{\partial f_i}{\partial x_3} = -\frac{x_0 (u_i² + u_i x_1)}{(u_i² + u_i x_2 + x_3)²} \end{align}

我们将使用在[2]中定义的“难点”起始点。为了找到物理上有意义的解，避免可能的除零，并确保收敛到全局最小值，我们施加约束条件(0 \leq x_j \leq 100, j = 0, 1, 2, 3)。

下面的代码实现了对(\mathbf{x})的最小二乘估计，并最终绘制了原始数据和拟合的模型函数：

>>> 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)
# may vary
`ftol` termination condition is satisfied.
Function evaluations 130, initial cost 4.4383e+00, final cost 1.5375e-04, first-order optimality 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='data')
>>> plt.plot(u_test, y_test, label='fitted model')
>>> plt.xlabel("u")
>>> plt.ylabel("y")
>>> plt.legend(loc='lower right')
>>> plt.show() 

更多示例

下面的三个互动示例详细说明了如何使用least_squares。

1. Scipy 中的大规模束调整展示了least_squares的大规模能力，以及如何高效地计算稀疏雅可比矩阵的有限差分近似。

2. Scipy 中的鲁棒非线性回归展示了如何在非线性回归中使用鲁棒损失函数处理异常值。

3. Scipy 中解离散边值问题介绍了如何解决大型方程系统，并使用边界来实现解的期望性质。

有关实现背后数学算法的详细信息，请参阅least_squares的文档。

单变量函数最小化器 (minimize_scalar)

经常只需要对单变量函数（即以标量作为输入的函数）进行最小化。在这些情况下，已经开发出了其他可以更快工作的优化技术。这些技术可以从minimize_scalar函数中访问，该函数提供了几种算法。

无约束最小化（method='brent'）

实际上，有两种方法可以用来最小化单变量函数：brent和golden，但golden仅用于学术目的，应该很少使用。这些可以通过minimize_scalar中的method参数分别选择。brent方法使用 Brent 算法来寻找最小值。理想情况下，应该提供一个包含所需最小值的区间（bracket参数），这是一个三元组 (\left( a, b, c \right))，使得 (f \left( a \right) > f \left( b \right) < f \left( c \right)) 并且 (a < b < c)。如果未提供这个参数，则可以选择两个起始点，并且使用简单的逐步算法从这些点中找到一个区间。如果没有提供这两个起始点，则会使用 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 函数来找到一组非线性方程的根。有多种方法可用，其中包括使用 Powell 的混合方法和 MINPACK 中的 Levenberg-Marquardt 方法的 hybr（默认）和 lm。

以下示例考虑单变量超越方程

[x+2\cos\left(x\right)=0,]

可以如下找到一个根：

>>> 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]) 

现在考虑一组非线性方程：

\begin{eqnarray*} x_{0}\cos\left(x_{1}\right) & = & 4,\ x_{0}x_{1}-x_{1} & = & 5. \end{eqnarray*}

我们定义目标函数，使其返回雅可比矩阵，并通过设置 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]) 

解决大问题的根查找

方法 hybr 和 lm 在 root 中不能处理非常大数量的变量 (N)，因为它们需要在每次牛顿步骤中计算并求逆密集的 N x N 雅可比矩阵，当 N 增大时效率变得相当低下。

例如，考虑在正方形 ([0,1]\times[0,1]) 上解决以下积分微分方程：

[(\partial_x² + \partial_y²) P + 5 \left(\int_0¹\int_0¹\cosh(P),dx,dy\right)² = 0]

在边界条件 (P(x,1) = 1) 上以及方形边界上的其他地方 (P=0)。可以通过在网格上近似连续函数 P 的值 (P_{n,m}\approx{}P(n h, m h))，其中网格间距 h 很小，来完成此操作。然后可以近似导数和积分；例如 (\partial_x² P(x,y)\approx{}(P(x+h,y) - 2 P(x,y) + P(x-h,y))/h²)。然后，问题等效于找到某些函数 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

# 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 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

# solve
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())

# visualize
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 求解器大部分时间用于求解雅可比矩阵的逆。

[J_{ij} = \frac{\partial f_i}{\partial x_j} .]

如果你对逆矩阵 (M\approx{}J^{-1}) 有一个近似值，可以用它来预处理线性反演问题。其思想是，不是解决 (J{\bf s}={\bf y})，而是解决 (MJ{\bf s}=M{\bf y})：因为矩阵 (MJ) 比 (J) 更接近单位矩阵，所以对于 Krylov 方法来说，这个方程应该更容易处理。

矩阵M可以作为选项options['jac_options']['inner_M']传递给root的krylov方法。它可以是一个（稀疏）矩阵或scipy.sparse.linalg.LinearOperator实例。

对于前一节中的问题，我们注意到要解决的函数由两部分组成：第一部分是拉普拉斯算子的应用，([\partial_x² + \partial_y²] P)，第二部分是积分。实际上，我们可以很容易地计算与拉普拉斯算子部分对应的雅可比矩阵：我们知道在一维中

[\begin{split}\partial_x² \approx \frac{1}{h_x²} \begin{pmatrix} -2 & 1 & 0 & 0 \cdots \ 1 & -2 & 1 & 0 \cdots \ 0 & 1 & -2 & 1 \cdots \ \ldots \end{pmatrix} = h_x^{-2} L\end{split}]

使得整个二维算子表示为

[J_1 = \partial_x² + \partial_y² \simeq h_x^{-2} L \otimes I + h_y^{-2} I \otimes L]

与积分对应的雅可比矩阵(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; step 1; tol 2.87308e-12
Residual 9.29603061195e-11
Evaluations 77 

使用预处理器将residual函数的评估次数减少了4倍。对于计算成本高昂的残差的问题，良好的预处理至关重要 —— 它甚至可以决定问题在实践中是否可解。

预处理是一门艺术、科学和工业。在这里，我们很幸运地做出了一个简单的选择，效果还不错，但这个主题比这里展示的要深入得多。

线性规划（linprog)

函数linprog可以最小化一个线性目标函数，同时满足线性等式和不等式约束。这种问题被称为线性规划。线性规划解决以下形式的问题：

[\begin{split}\min_x \ & c^T x \ \mbox{such that} \ & A_{ub} x \leq b_{ub},\ & A_{eq} x = b_{eq},\ & l \leq x \leq u ,\end{split}]

其中(x)是决策变量向量；(c)、(b_{ub})、(b_{eq})、(l)和(u)是向量；(A_{ub})和(A_{eq})是矩阵。

在本教程中，我们将尝试使用linprog解决典型的线性规划问题。

线性规划示例

考虑以下简单的线性规划问题：

[\begin{split}\max_{x_1, x_2, x_3, x_4} \ & 29x_1 + 45x_2 \ \mbox{such that} \ & x_1 -x_2 -3x_3 \leq 5\ & 2x_1 -3x_2 -7x_3 + 3x_4 \geq 10\ & 2x_1 + 8x_2 + x_3 = 60\ & 4x_1 + 4x_2 + x_4 = 60\ & 0 \leq x_0\ & 0 \leq x_1 \leq 5\ & x_2 \leq 0.5\ & -3 \leq x_3\\end{split}]

我们需要一些数学操作来将目标问题转换为linprog接受的形式。

首先，让我们考虑目标函数。我们想要最大化目标函数，但linprog只能接受最小化问题。这很容易通过将最大化(29x_1 + 45x_2)转换为最小化(-29x_1 -45x_2)来修正。另外，(x_3, x_4)在目标函数中没有显示。这意味着与(x_3, x_4)对应的权重为零。因此，目标函数可以转换为：

[\min_{x_1, x_2, x_3, x_4} \ -29x_1 -45x_2 + 0x_3 + 0x_4]

如果我们定义决策变量向量(x = [x_1, x_2, x_3, x_4]^T)，那么linprog在这个问题中的目标权重向量(c)应为

[c = [-29, -45, 0, 0]^T]

接下来，让我们考虑这两个不等式约束。第一个是“小于”不等式，因此已经是linprog接受的形式。第二个是“大于”不等式，因此我们需要将两边乘以(-1)将其转换为“小于”不等式。明确显示零系数，我们有：

[\begin{split}x_1 -x_2 -3x_3 + 0x_4 &\leq 5\ -2x_1 + 3x_2 + 7x_3 - 3x_4 &\leq -10\\end{split}]

这些方程可以转换为矩阵形式：

[\begin{split}A_{ub} x \leq b_{ub}\\end{split}]

其中

\begin{equation*} A_{ub} = \begin{bmatrix} 1 & -1 & -3 & 0 \ -2 & 3 & 7 & -3 \end{bmatrix} \end{equation*}\begin{equation*} b_{ub} = \begin{bmatrix} 5 \ -10 \end{bmatrix} \end{equation*}

接下来，让我们考虑两个等式约束。明确显示零权重，它们是：

[\begin{split}2x_1 + 8x_2 + 1x_3 + 0x_4 &= 60\ 4x_1 + 4x_2 + 0x_3 + 1x_4 &= 60\\end{split}]

这些方程可以转换为矩阵形式：

[\begin{split}A_{eq} x = b_{eq}\\end{split}]

其中

\begin{equation*} A_{eq} = \begin{bmatrix} 2 & 8 & 1 & 0 \ 4 & 4 & 0 & 1 \end{bmatrix} \end{equation*}\begin{equation*} b_{eq} = \begin{bmatrix} 60 \ 60 \end{bmatrix} \end{equation*}

最后，让我们考虑对单独决策变量的分离不等式约束，这些约束被称为“箱约束”或“简单边界”。这些约束可以使用 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 can be used instead of 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  # may vary
>>> print(obj)
-505.97435889013434  # may vary 

我们还可以检查所有约束是否在合理的容差范围内得到满足：

>>> print(b_ub - (A_ub @ x).flatten())  # this is equivalent to result.slack
[ 6.52747190e-10, -2.26730279e-09]  # may vary
>>> print(b_eq - (A_eq @ x).flatten())  # this is equivalent to result.con
[ 9.78840831e-09, 1.04662945e-08]]  # may vary
>>> 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] 

分配问题

线性和分配问题示例

考虑选择一个游泳混合接力队的学生问题。我们有一个表格显示了五名学生每种游泳风格的时间：

Student	backstroke	breaststroke	butterfly	freestyle
A	43.5	47.1	48.4	38.2
B	45.5	42.1	49.6	36.8
C	43.4	39.1	42.1	43.2
D	46.5	44.1	44.5	41.2
E	46.3	47.8	50.4	37.2

我们需要为每种四种游泳风格选择一个学生，以使接力总时间最小化。这是一个典型的线性和分配问题。我们可以使用 linear_sum_assignment 来解决它。

线性和分配问题是最著名的组合优化问题之一。给定一个“成本矩阵” (C)，问题是选择每行中的一个元素

• 从而确保每列中不选择超过一个元素

• 而不选择任何列中超过一个元素

• 以使所选元素的和最小化

换句话说，我们需要将每一行分配给一个列，使得对应条目的总和最小化。

正式地说，设 (X) 是一个布尔矩阵，其中 (X[i,j] = 1) 当且仅当第 (i) 行被分配给第 (j) 列。那么最优的分配成本为

[\min \sum_i \sum_j C_{i,j} X_{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 分配到“蝶泳”项目以最小化总时间。

参考文献

一些进一步阅读和相关软件，例如牛顿-克里洛夫 [KK]，PETSc [PP]，和 PyAMG [AMG]：

[KK]

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

[PP]

PETSc www.mcs.anl.gov/petsc/ 和其 Python 绑定 bitbucket.org/petsc/petsc4py/

[AMG]

PyAMG（代数多重网格预处理器/求解器）github.com/pyamg/pyamg/issues

混合整数线性规划

背包问题示例

背包问题是一个著名的组合优化问题。给定一组物品，每个物品有一个大小和一个价值，问题是在总大小不超过一定阈值的条件下选择物品以最大化总价值。

正式地说，设

• (x_i) 是一个布尔变量，表示是否将物品 (i) 放入背包，

• (n) 表示物品的总数，

• (v_i) 是物品 (i) 的价值，

• (s_i) 表示物品 (i) 的大小，以及

• (C) 表示背包的容量。

然后问题是:

[\max \sum_i^n v_{i} x_{i}][\text{subject to} \sum_i^n s_{i} x_{i} \leq C, x_{i} \in {0, 1}]

尽管目标函数和不等式约束在决策变量 (x_i) 中是线性的，但这与典型的线性规划问题不同，因为决策变量只能取整数值。具体来说，我们的决策变量只能是 (0) 或 (1)，因此这被称为二进制整数线性规划（BILP）。这种问题属于更大的混合整数线性规划（MILP）类别，我们可以使用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 are integers 

使用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 Cookbook 上的 Jupyter 笔记本：

• 压缩感知 l1 程序

• 压缩感知 l0 程序