SciPy-1-12-中文文档-一-

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

"This code generates a 3-D representation of the vibrational modes on a drum head viewed at a three-quarter angle. A circular region on the X-Y plane is defined with a Z value of 0 around the edge. Within the circle a single smooth valley exists on the -X side and a smooth peak exists on the +X side. The image resembles a yin-yang at this angle."

特殊函数的 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\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}\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} 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}\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} 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}\right){2}+\left(1-x_{i}\right).]

此函数的最小值为 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}\right){2}+\left(1-x_{i}\right) + 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}\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}\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}\right)\cdot\left(\mathbf{x}-\mathbf{x}{0}\right)+\frac{1}{2}\left(\mathbf{x}-\mathbf{x}\right)^{T}\mathbf{H}\left(\mathbf{x}{0}\right)\left(\mathbf{x}-\mathbf{x}\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}-\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}+202\mkern-2em\&&1200x_{1}{2}+202\mkern-2em\&&&1200x_{3}+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}\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} + \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}\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)+\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()

"这个 X-Y 图是一个热图，Z 值用最低点为黑色，最高值为白色表示。图像类似于旋转了 45 度但严重平滑的棋盘格。一个红色的点位于由 SHGO 优化器产生的许多最小值格点上。SHGO 在右上角显示全局最小值为红色 X。用双退火找到的局部最小值为左上角的白色圆形标记。用盆地跳跃找到的另一个局部最小值为顶部中心的黄色标记。代码正在绘制差分进化结果作为青色圆圈，但在图上看不见。乍一看不清楚哪个山谷是真正的全局最小值。"

最小二乘最小化（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()

"此代码绘制了一个 X-Y 时间序列。该序列从左下角的(0, 0)开始，迅速趋向最大值 0.2，然后趋于平缓。拟合模型显示为平滑的橙色曲线，与数据非常匹配。"

更多示例

下面的三个互动示例详细说明了如何使用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()

"This code generates a 2-D heatmap with Z values from 0 to 1. The graph resembles a smooth, dark blue-green, U shape, with an open yellow top. The right, bottom, and left edges have a value near zero and the top has a value close to 1. The center of the solution space has a value close to 0.8."

还是太慢？预处理。

当寻找函数 (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 笔记本：

插值（scipy.interpolate）

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

SciPy 提供了几种用于 1、2 和更高维数据的插值和平滑的通用工具。选择特定的插值程序取决于数据：是一维数据、是在结构化网格上给定，还是非结构化的。另一个因素是插值器的期望平滑度。简而言之，推荐用于插值的程序可以总结如下：

种类 例程 连续性 备注
1D 线性 numpy.interp 分段连续 来自 numpy
三次样条 CubicSpline 2 阶导数
单调三次样条 PchipInterpolator 1 阶导数 不越界
非立方样条 make_interp_spline (k-1)阶导数 k=3 等同于 CubicSpline
最近邻 interp1d 种类='nearest'、'previous'、'next'
N-D 曲线 最近邻、线性、样条 make_interp_spline (k-1)阶导数 使用 N 维 y 数组
N-D 规则（矩形）网格 最近邻 RegularGridInterpolator 方法='nearest'
线性 方法='linear'
样条 2 阶导数 方法='cubic'、'quintic'
单调样条 1 阶导数 方法='pchip'
N-D 散点 最近邻 NearestNDInterpolator 别名：griddata
线性 LinearNDInterpolator
立方体（仅限 2D） CloughTocher2DInterpolator 1st derivatives
径向基函数 RBFInterpolator

对于数据平滑，提供了使用三次样条的基于 FORTRAN 库 FITPACK 的 1-D 和 2-D 数据的功能。

此外，还提供了使用径向基函数进行插值/平滑的例程，其中包括几种核函数。

更多细节请查看以下链接。

  • 1-D 插值

    • 分段线性插值

    • 三次样条

    • 单调插值器

    • B 样条插值

    • 参数样条曲线

    • 1-D 插值的旧接口 (interp1d)

    • 缺失数据

  • 分段多项式和样条

    • 操纵 PPoly 对象

    • B 样条：节点和系数

      • B 样条基函数

      • B 样条基设计矩阵

  • 平滑样条

    • 1-D 中的样条平滑

      • 程序化 (splrep)

      • 面向对象的 (UnivariateSpline)

    • 2-D 平滑样条

      • 散点数据的双变量样条拟合

      • 在网格上数据的双变量样条拟合

      • 球坐标中的数据的双变量样条拟合

  • 在规则网格上的多变量数据插值 (RegularGridInterpolator)

    • 均匀间隔数据

  • 散点数据插值 (griddata)

  • 使用径向基函数进行平滑/插值](interpolate/ND_unstructured.html#using-radial-basis-functions-for-smoothing-interpolation)

    • 1-D 示例

    • 2-D 示例

  • 外推技巧和窍门](interpolate/extrapolation_examples.html)

    • interp1d：复制numpy.interp的左右填充值

    • CubicSpline 扩展边界条件

    • 手动实现渐近值](interpolate/extrapolation_examples.html#manually-implement-the-asymptotics)

      • 设置

      • 使用已知的渐近值](interpolate/extrapolation_examples.html#use-the-known-asymptotics)

    • 在D > 1中的外推

傅里叶变换（scipy.fft）

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

内容

  • 傅里叶变换（scipy.fft)

    • 快速傅里叶变换

      • 1D 离散傅里叶变换

      • 2D 和 ND 离散傅里叶变换

    • 离散余弦变换

      • Type I DCT

      • Type II DCT

      • Type III DCT

      • Type IV DCT

      • DCT 和 IDCT

      • 示例

    • 离散正弦变换

      • Type I DST

      • Type II DST

      • Type III DST

      • Type IV DST

      • DST 和 IDST

    • 快速汉克尔变换

    • 参考文献

Fourier 分析是将函数表达为周期分量之和并从这些分量中恢复信号的方法。当函数及其傅里叶变换被替换为离散化的对应物时，称为离散傅里叶变换（DFT）。由于一种名为快速傅里叶变换（FFT）的非常快速计算算法的存在，DFT 已成为数值计算的重要工具之一，该算法由高斯（1805 年）知晓，并由库利和图基在现代形式中首次揭示[CT65]。Press 等人提供了傅里叶分析及其应用的易于理解的介绍[NR07]。

快速傅里叶变换

1D 离散傅里叶变换

FFT y[k] 长度为(N)的序列 x[n] 的傅里叶变换定义为

[y[k] = \sum_{n=0}^{N-1} e^{-2 \pi j \frac{k n}{N} } x[n] , ,]

反变换定义如下

[x[n] = \frac{1}{N} \sum_{k=0}^{N-1} e^{2 \pi j \frac{k n}{N} } y[k] , .]

这些变换可以通过fft和ifft来计算，如下例所示。

>>> from scipy.fft import fft, ifft
>>> import numpy as np
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
>>> y = fft(x)
>>> y
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
 -1.83155948+1.60822041j, -1.83155948-1.60822041j,
 2.08155948+1.65109876j])
>>> yinv = ifft(y)
>>> yinv
array([ 1.0+0.j,  2.0+0.j,  1.0+0.j, -1.0+0.j,  1.5+0.j])

从 FFT 的定义可以看出

[y[0] = \sum_{n=0}^{N-1} x[n] , .]

在例子中

>>> np.sum(x)
4.5

对应于(y[0])。对于偶数 N，元素(y[1]...y[N/2-1])包含正频率项，元素(y[N/2]...y[N-1])包含负频率项，按递减负频率顺序排列。对于奇数 N，元素(y[1]...y[(N-1)/2])包含正频率项，元素(y[(N+1)/2]...y[N-1])包含负频率项，按递减负频率顺序排列。

如果序列 x 是实值的，则正频率下的(y[n])值是负频率下(y[n])值的共轭（因为频谱是对称的）。通常只绘制与正频率对应的 FFT。

该示例绘制了两个正弦波的 FFT。

>>> from scipy.fft import fft, fftfreq
>>> import numpy as np
>>> # Number of sample points
>>> N = 600
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
>>> yf = fft(y)
>>> xf = fftfreq(N, T)[:N//2]
>>> import matplotlib.pyplot as plt
>>> plt.plot(xf, 2.0/N * np.abs(yf[0:N//2]))
>>> plt.grid()
>>> plt.show()

"这段代码生成一个振幅在 Y 轴上，频率在 X 轴上的 X-Y 图。单个蓝色迹象在整个图中除了两个峰值外，其余位置振幅为零。第一个更高的峰值在 50 Hz 处，第二个峰值在 80 Hz 处。"

FFT 输入信号本质上是被截断的。这种截断可以被建模为无限信号与矩形窗口函数的乘积。在频谱域中，这种乘积变成了信号频谱与窗口函数频谱的卷积，形式为(\sin(x)/x)。这种卷积引起了称为频谱泄漏的效应（参见[WPW]）。使用专用窗口函数对信号进行窗函数处理有助于减轻频谱泄漏。下面的示例使用了 scipy.signal 中的 Blackman 窗口，并展示了窗函数效果（FFT 的零分量已被截断，仅用于说明目的）。

>>> from scipy.fft import fft, fftfreq
>>> import numpy as np
>>> # Number of sample points
>>> N = 600
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.sin(50.0 * 2.0*np.pi*x) + 0.5*np.sin(80.0 * 2.0*np.pi*x)
>>> yf = fft(y)
>>> from scipy.signal.windows import blackman
>>> w = blackman(N)
>>> ywf = fft(y*w)
>>> xf = fftfreq(N, T)[:N//2]
>>> import matplotlib.pyplot as plt
>>> plt.semilogy(xf[1:N//2], 2.0/N * np.abs(yf[1:N//2]), '-b')
>>> plt.semilogy(xf[1:N//2], 2.0/N * np.abs(ywf[1:N//2]), '-r')
>>> plt.legend(['FFT', 'FFT w. window'])
>>> plt.grid()
>>> plt.show()

"这段代码生成一个振幅在 Y 轴上，频率在 X 轴上的 X-Y 对数线性图。第一个迹象是 FFT，在 50 和 80 Hz 有两个峰值，噪声底噪约为 1e-2。第二个迹象是窗口 FFT，具有相同的两个峰，但由于窗口函数，噪声底噪较低，约为 1e-7。"

如果序列 x 是复值的，则频谱不再对称。为了简化使用 FFT 函数的工作，scipy 提供了以下两个辅助函数。

函数fftfreq返回 FFT 样本频率点。

>>> from scipy.fft import fftfreq
>>> freq = fftfreq(8, 0.125)
>>> freq
array([ 0., 1., 2., 3., -4., -3., -2., -1.])

类似地，函数fftshift允许交换向量的下半部分和上半部分，以便于显示。

>>> from scipy.fft import fftshift
>>> x = np.arange(8)
>>> fftshift(x)
array([4, 5, 6, 7, 0, 1, 2, 3])

下面的示例绘制了两个复指数的 FFT；请注意其不对称的频谱。

>>> from scipy.fft import fft, fftfreq, fftshift
>>> import numpy as np
>>> # number of signal points
>>> N = 400
>>> # sample spacing
>>> T = 1.0 / 800.0
>>> x = np.linspace(0.0, N*T, N, endpoint=False)
>>> y = np.exp(50.0 * 1.j * 2.0*np.pi*x) + 0.5*np.exp(-80.0 * 1.j * 2.0*np.pi*x)
>>> yf = fft(y)
>>> xf = fftfreq(N, T)
>>> xf = fftshift(xf)
>>> yplot = fftshift(yf)
>>> import matplotlib.pyplot as plt
>>> plt.plot(xf, 1.0/N * np.abs(yplot))
>>> plt.grid()
>>> plt.show()

"这段代码生成一个振幅在 Y 轴上，频率在 X 轴上的 X-Y 图。在整个图中，迹象除了在-80 和 50 Hz 处有两个尖峰外，其他位置值为零。右侧的 50 Hz 峰值是左侧峰值的两倍。"

函数rfft计算实序列的 FFT，并仅输出半频率范围的复 FFT 系数(y[n])。剩余的负频率分量由 FFT 的 Hermitian 对称性隐含处理，用于实输入的情况（y[n] = conj(y[-n])）。对于 N 为偶数：([Re(y[0]) + 0j, y[1], ..., Re(y[N/2]) + 0j])；对于 N 为奇数：([Re(y[0]) + 0j, y[1], ..., y[N/2]]。显示为(Re(y[k]) + 0j)的项受到限制，因为根据 Hermitian 属性，它们是它们自己的复共轭。

函数irfft用特殊顺序计算 FFT 系数的 IFFT。

>>> from scipy.fft import fft, rfft, irfft
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5, 1.0])
>>> fft(x)
array([ 5.5 +0.j        ,  2.25-0.4330127j , -2.75-1.29903811j,
 1.5 +0.j        , -2.75+1.29903811j,  2.25+0.4330127j ])
>>> yr = rfft(x)
>>> yr
array([ 5.5 +0.j        ,  2.25-0.4330127j , -2.75-1.29903811j,
 1.5 +0.j        ])
>>> irfft(yr)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5,  1\. ])
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])
>>> fft(x)
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
 -1.83155948+1.60822041j, -1.83155948-1.60822041j,
 2.08155948+1.65109876j])
>>> yr = rfft(x)
>>> yr
array([ 4.5       +0.j        ,  2.08155948-1.65109876j,
 -1.83155948+1.60822041j])

注意rfft对于奇偶长度的信号具有相同的形状。默认情况下，irfft假定输出信号应为偶数长度。因此，对于奇数信号，它将给出错误的结果：

>>> irfft(yr)
array([ 1.70788987,  2.40843925, -0.37366961,  0.75734049])

要恢复原始的奇数长度信号，必须通过n参数传递输出形状。

>>> irfft(yr, n=len(x))
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

2- and N-D 离散傅里叶变换

函数fft2和ifft2分别提供 2-D FFT 和 IFFT。类似地，fftn和ifftn提供 N-D FFT 和 IFFT。

对于实输入信号，类似于rfft，我们有函数rfft2和irfft2用于 2-D 实变换；rfftn和irfftn用于 N-D 实变换。

下面的示例演示了 2-D IFFT 并绘制了结果的（2-D）时域信号。

>>> from scipy.fft import ifftn
>>> import matplotlib.pyplot as plt
>>> import matplotlib.cm as cm
>>> import numpy as np
>>> N = 30
>>> f, ((ax1, ax2, ax3), (ax4, ax5, ax6)) = plt.subplots(2, 3, sharex='col', sharey='row')
>>> xf = np.zeros((N,N))
>>> xf[0, 5] = 1
>>> xf[0, N-5] = 1
>>> Z = ifftn(xf)
>>> ax1.imshow(xf, cmap=cm.Reds)
>>> ax4.imshow(np.real(Z), cmap=cm.gray)
>>> xf = np.zeros((N, N))
>>> xf[5, 0] = 1
>>> xf[N-5, 0] = 1
>>> Z = ifftn(xf)
>>> ax2.imshow(xf, cmap=cm.Reds)
>>> ax5.imshow(np.real(Z), cmap=cm.gray)
>>> xf = np.zeros((N, N))
>>> xf[5, 10] = 1
>>> xf[N-5, N-10] = 1
>>> Z = ifftn(xf)
>>> ax3.imshow(xf, cmap=cm.Reds)
>>> ax6.imshow(np.real(Z), cmap=cm.gray)
>>> plt.show()

"这段代码生成了一个 2x3 网格中排列的六个热图。顶行显示大多是空白画布，除了每幅图像上有两个微小的红色峰。底行显示了每个上面图像的逆 FFT 的实部。第一列有两个水平排列的点在上面的图像中，在下面的图像中有一个平滑的灰度图，显示了五条黑色垂直条纹，代表 2D 时域信号。第二列有两个垂直排列的点在上面的图像中，在下面的图像中有一个平滑的灰度图，显示了五条水平黑色条纹，代表 2D 时域信号。在最后一列中，顶部图像有两个对角线位置的点；对应的下面的图像有大概 20 条 60 度角的黑色条纹。"

离散余弦变换

SciPy 提供了函数 dct 和相应的 IDCT 函数 idct 进行离散余弦变换。总共有 8 种 DCT 类型[WPC], [Mak]; 然而，SciPy 只实现了前 4 种。通常“DCT” 指的是 DCT 类型 2，“反向 DCT” 指的是 DCT 类型 3. 此外，DCT 系数可以有不同的标准化方式（对于大多数类型，SciPy 提供 None 和 ortho）。dct/idct 函数调用的两个参数允许设置 DCT 类型和系数标准化。

对于一维数组 x，使用 norm='ortho' 的 dct(x) 等同于 MATLAB 的 dct(x)。

第一类离散余弦变换

SciPy 使用以下未归一化 DCT-I 的定义 (norm=None)：

[y[k] = x_0 + (-1)^k x_{N-1} + 2\sum_{n=1}^{N-2} x[n] \cos\left(\frac{\pi nk}{N-1}\right), \qquad 0 \le k < N.]

注意，DCT-I 仅支持输入大小 > 1。

第二类离散余弦变换

SciPy 使用以下未归一化 DCT-II 的定义 (norm=None)：

[y[k] = 2 \sum_{n=0}^{N-1} x[n] \cos \left({\pi(2n+1)k \over 2N} \right) \qquad 0 \le k < N.]

在归一化 DCT 的情况下 (norm='ortho')，DCT 系数 (y[k]) 乘以一个缩放因子 f：

[\begin{split}f = \begin{cases} \sqrt{1/(4N)}, & \text{如果 \(k = 0\)} \ \sqrt{1/(2N)}, & \text{否则} \end{cases} , .\end{split}]

在这种情况下，DCT 的“基函数” (\phi_k[n] = 2 f \cos \left({\pi(2n+1)k \over 2N} \right)) 变为正交的：

[\sum_{n=0}^{N-1} \phi_k[n] \phi_l[n] = \delta_{lk}.]

第三类离散余弦变换

SciPy 使用以下未归一化 DCT-III 的定义 (norm=None)：

[y[k] = x_0 + 2 \sum_{n=1}^{N-1} x[n] \cos\left({\pi n(2k+1) \over 2N}\right) \qquad 0 \le k < N,]

or, for norm='ortho':

[y[k] = {x_0\over\sqrt{N}} + {2\over\sqrt{N}} \sum_{n=1}^{N-1} x[n] \cos\left({\pi n(2k+1) \over 2N}\right) \qquad 0 \le k < N.]

第四类离散余弦变换

SciPy 使用以下未归一化 DCT-IV 的定义 (norm=None)：

[y[k] = 2 \sum_{n=0}^{N-1} x[n] \cos\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,]

或者，对于 norm='ortho'：

[y[k] = \sqrt{2\over N}\sum_{n=0}^{N-1} x[n] \cos\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N]

DCT 和 IDCT

（未标准化的）DCT-III 是（未标准化的）DCT-II 的反变换，只差一个因子 2N。标准化的 DCT-III 则正好是标准化的 DCT-II 的反变换。函数 idct 执行了 DCT 和 IDCT 类型之间的映射，以及正确的归一化。

以下示例展示了不同类型和归一化的 DCT 和 IDCT 之间的关系。

>>> from scipy.fft import dct, idct
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

DCT-II 和 DCT-III 是彼此的反变换，因此对于正交变换，我们可以恢复到原始信号。

>>> dct(dct(x, type=2, norm='ortho'), type=3, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

在默认归一化下进行相同操作时，我们得到一个额外的缩放因子 (2N=10)，因为正向变换是未标准化的。

>>> dct(dct(x, type=2), type=3)
array([ 10.,  20.,  10., -10.,  15.])

因此，我们应该使用函数 idct 为两者使用相同类型，以获得正确归一化的结果。

>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=2), type=2)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

对于 DCT-I，类似的结果可以看到它本身的反变换，只差一个因子 (2(N-1))。

>>> dct(dct(x, type=1, norm='ortho'), type=1, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])
>>> # Unnormalized round-trip via DCT-I: scaling factor 2*(N-1) = 8
>>> dct(dct(x, type=1), type=1)
array([ 8\. ,  16.,  8\. , -8\. ,  12.])
>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=1), type=1)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

对于 DCT-IV，也是其自身的反变换，只差一个因子 (2N)。

>>> dct(dct(x, type=4, norm='ortho'), type=4, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])
>>> # Unnormalized round-trip via DCT-IV: scaling factor 2*N = 10
>>> dct(dct(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>> # Normalized inverse: no scaling factor
>>> idct(dct(x, type=4), type=4)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

示例

DCT 表现出“能量压缩特性”，意味着对于许多信号，只有前几个 DCT 系数具有显著的幅度。将其他系数置零会导致小的重构误差，这一事实在损失信号压缩（例如 JPEG 压缩）中得到利用。

下面的示例展示了信号 x 及其从 DCT 系数重构的两个重构（(x_{20}) 和 (x_{15})）。从使用 20 个系数重构的信号 (x_{20}) 中可以看出，相对误差仍然非常小（约 0.1%），但提供了五倍的压缩率。

>>> from scipy.fft import dct, idct
>>> import matplotlib.pyplot as plt
>>> N = 100
>>> t = np.linspace(0,20,N, endpoint=False)
>>> x = np.exp(-t/3)*np.cos(2*t)
>>> y = dct(x, norm='ortho')
>>> window = np.zeros(N)
>>> window[:20] = 1
>>> yr = idct(y*window, norm='ortho')
>>> sum(abs(x-yr)**2) / sum(abs(x)**2)
0.0009872817275276098
>>> plt.plot(t, x, '-bx')
>>> plt.plot(t, yr, 'ro')
>>> window = np.zeros(N)
>>> window[:15] = 1
>>> yr = idct(y*window, norm='ortho')
>>> sum(abs(x-yr)**2) / sum(abs(x)**2)
0.06196643004256714
>>> plt.plot(t, yr, 'g+')
>>> plt.legend(['x', '$x_{20}$', '$x_{15}$'])
>>> plt.grid()
>>> plt.show()

"此代码生成一个 X-Y 图，显示 Y 轴上的振幅和 X 轴上的时间。第一个蓝色轨迹是原始信号，从振幅 1 开始，在绘图持续时间内振幅下降至 0，类似频率啁啾。第二个红色轨迹是使用 DCT 生成的 x_20 重构，高振幅区域紧随原始信号，但在绘图右侧不够清晰。第三个绿色轨迹是使用 DCT 生成的 x_15 重构，比 x_20 重构不够精确，但仍类似于 x。"

离散正弦变换

SciPy 提供了使用函数 dst 进行 DST（Mak）计算，并使用函数 idst 进行对应的 IDST 计算。

理论上，有 8 种 DST 类型，适用于不同的偶数/奇数边界条件和边界偏移组合[WPS]，但 scipy 只实现了前 4 种类型。

类型 I DST

DST-I 假设输入在 n=-1 和 n=N 周围为奇数。SciPy 使用以下未归一化定义 DST-I (norm=None):

[y[k] = 2\sum_{n=0}^{N-1} x[n] \sin\left( \pi {(n+1) (k+1)}\over{N+1} \right), \qquad 0 \le k < N.]

还要注意，DST-I 仅支持输入大小 > 1。DST-I（未归一化）是其自身的逆变换，除了一个 2(N+1) 的因子。

类型 II DST

DST-II 假设输入在 n=-1/2 周围为奇数，在 n=N 周围为偶数。SciPy 使用以下未归一化定义 DST-II (norm=None):

[y[k] = 2 \sum_{n=0}^{N-1} x[n] \sin\left( {\pi (n+1/2)(k+1)} \over N \right), \qquad 0 \le k < N.]

类型 III DST

DST-III 假设输入在 n=-1 周围为奇数，在 n=N-1 周围为偶数。SciPy 使用以下未归一化定义 DST-III (norm=None):

[y[k] = (-1)^k x[N-1] + 2 \sum_{n=0}^{N-2} x[n] \sin \left( {\pi (n+1)(k+1/2)} \over N \right), \qquad 0 \le k < N.]

类型 IV DST

SciPy 使用以下未归一化定义 DST-IV (norm=None):

[y[k] = 2 \sum_{n=0}^{N-1} x[n] \sin\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,]

或者，对于 norm='ortho'：

[y[k] = \sqrt{2\over N}\sum_{n=0}^{N-1} x[n] \sin\left({\pi (2n+1)(2k+1) \over 4N}\right) \qquad 0 \le k < N,]

DST 和 IDST

以下示例展示了不同类型和标准化下 DST 与 IDST 之间的关系。

>>> from scipy.fft import dst, idst
>>> x = np.array([1.0, 2.0, 1.0, -1.0, 1.5])

DST-II 和 DST-III 是彼此的逆变换，因此对于正交变换，我们可以返回到原始信号。

>>> dst(dst(x, type=2, norm='ortho'), type=3, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

在默认标准化下进行相同处理时，由于正向变换未归一化，我们会得到一个额外的缩放因子 (2N=10)。

>>> dst(dst(x, type=2), type=3)
array([ 10.,  20.,  10., -10.,  15.])

因此，我们应该使用函数 idst 为两者使用相同的类型，从而得到正确标准化的结果。

>>> idst(dst(x, type=2), type=2)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

对于 DST-I，其本身是其逆变换，仅相差一个因子 (2(N-1))。

>>> dst(dst(x, type=1, norm='ortho'), type=1, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])
>>>  # scaling factor 2*(N+1) = 12
>>> dst(dst(x, type=1), type=1)
array([ 12.,  24.,  12., -12.,  18.])
>>>  # no scaling factor
>>> idst(dst(x, type=1), type=1)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

而对于 DST-IV，其本身也是其逆变换，除了一个因子 (2N)。

>>> dst(dst(x, type=4, norm='ortho'), type=4, norm='ortho')
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])
>>>  # scaling factor 2*N = 10
>>> dst(dst(x, type=4), type=4)
array([ 10.,  20.,  10., -10.,  15.])
>>>  # no scaling factor
>>> idst(dst(x, type=4), type=4)
array([ 1\. ,  2\. ,  1\. , -1\. ,  1.5])

快速汉克尔变换

SciPy 提供了函数 fht 和 ifht 来对对数间隔输入数组执行快速汉克尔变换（FHT）及其逆变换（IFHT）。

福特变换（FHT）是由[Ham00]定义的连续汉克尔变换的离散版本。

[A(k) = \int_{0}^{\infty} ! a(r) , J_{\mu}(kr) , k , dr ;,]

其中 (J_{\mu}) 是阶数为 (\mu) 的贝塞尔函数。在变量变换 (r \to \log r), (k \to \log k) 下，这变为

[A(e^{\log k}) = \int_{0}^{\infty} ! a(e^{\log r}) , J_{\mu}(e^{\log k + \log r}) , e^{\log k + \log r} , d{\log r}]

这是对数空间中的卷积。FHT 算法使用 FFT 对离散输入数据执行此卷积。

要注意由于 FFT 卷积的循环性质，必须小心减少数值环绕现象。为确保低环绕条件[Ham00]成立，可以通过使用fhtoffset函数计算的偏移量稍微偏移输出数组。

References

[CT65]

Cooley, James W., and John W. Tukey, 1965, “An algorithm for the machine calculation of complex Fourier series,” Math. Comput. 19: 297-301.

[NR07]

Press, W., Teukolsky, S., Vetterline, W.T., and Flannery, B.P., 2007, Numerical Recipes: The Art of Scientific Computing, ch. 12-13. Cambridge Univ. Press, Cambridge, UK.

[Mak] (1,2)

J. Makhoul, 1980, ‘A Fast Cosine Transform in One and Two Dimensions’, IEEE Transactions on acoustics, speech and signal processing vol. 28(1), pp. 27-34, DOI:10.1109/TASSP.1980.1163351

[Ham00] (1,2)

A. J. S. Hamilton, 2000, “Uncorrelated modes of the non-linear power spectrum”, MNRAS, 312, 257. DOI:10.1046/j.1365-8711.2000.03071.x

[WPW]

en.wikipedia.org/wiki/Window_function

[WPC]

en.wikipedia.org/wiki/Discrete_cosine_transform

[WPS]

en.wikipedia.org/wiki/Discrete_sine_transform

信号处理（scipy.signal）

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

信号处理工具箱目前包含一些滤波函数，有限的滤波器设计工具集，以及一些用于一维和二维数据的 B-样条插值算法。虽然 B-样条算法技术上可以归类为插值类别，但它们包含在此处是因为它们仅适用于等间距数据，并且广泛利用滤波理论和传递函数形式学来提供快速的 B-样条变换。要理解本节，您需要理解 SciPy 中信号是实数或复数数组。

B-样条

B-样条是在有限区间内以 B-样条系数和结点点的形式对连续函数的近似。如果结点点等间距，间距为 (\Delta x)，那么对一维函数的 B-样条近似是有限基函数展开。

[y\left(x\right)\approx\sum_{j}c_{j}\beta^{o}\left(\frac{x}{\Delta x}-j\right).]

在结点间距 (\Delta x) 和 (\Delta y) 的二维情况下，函数表示为

[z\left(x,y\right)\approx\sum_{j}\sum_{k}c_{jk}\beta^{o}\left(\frac{x}{\Delta x}-j\right)\beta^{o}\left(\frac{y}{\Delta y}-k\right).]

在这些表达式中，(\beta^{o}\left(\cdot\right)) 是阶数为 (o) 的空间限制 B-样条基函数。等间距结点和等间距数据点的要求，允许从样本值 (y_{n}) 中相对轻松地计算系数 (c_{j}) 的快速（逆滤波）算法。与一般样条插值算法不同，这些算法可以快速找到大图像的样条系数。

通过 B-样条基函数表示一组样本的优势在于，可以相对轻松地从样条系数中计算连续域操作（导数、重采样、积分等）。例如，样条的二阶导数是

[y{}^{\prime\prime}\left(x\right)=\frac{1}{\Delta x{2}}\sum_{j}c_{j}\beta\left(\frac{x}{\Delta x}-j\right).]

利用 B-样条的性质，

[\frac{d{2}\beta\left(w\right)}{dw{2}}=\beta\left(w+1\right)-2\beta{o-2}\left(w\right)+\beta\left(w-1\right),]

可以看出

[y^{\prime\prime}\left(x\right)=\frac{1}{\Delta x{2}}\sum_{j}c_{j}\left[\beta\left(\frac{x}{\Delta x}-j+1\right)-2\beta^{o-2}\left(\frac{x}{\Delta x}-j\right)+\beta^{o-2}\left(\frac{x}{\Delta x}-j-1\right)\right].]

如果 (o=3)，则在采样点：

\begin{eqnarray} \Delta x{2}\left.y\left(x\right)\right|{x=n\Delta x} & = & \sumc_{j}\delta_{n-j+1}-2c_{j}\delta_{n-j}+c_{j}\delta_{n-j-1},\ & = & c_{n+1}-2c_{n}+c_{n-1}.\end{eqnarray}

因此，可以轻松地从样条拟合中计算出二阶导数信号。如果需要，可以找到平滑样条，使得二阶导数对随机误差不太敏感。

精明的读者可能已经注意到，数据样本通过卷积算子与结节系数相关联，因此，简单地使用采样的 B 样条函数进行卷积可以从样条系数中恢复原始数据。卷积的输出可以根据边界处理方式而改变（随着数据集维度的增加，这一点变得越来越重要）。信号处理子包中关于 B 样条的算法假定镜像对称边界条件。因此，基于该假设计算样条系数，并且通过假设它们也是镜像对称的，可以精确地从样条系数中恢复数据样本。

目前，该包提供了在一维和二维等间距样本中确定二阶和三阶立方样条系数的函数（qspline1d, qspline2d, cspline1d, cspline2d)。对于较大的 (o)，B 样条基函数可以很好地近似为标准偏差为 (\sigma_{o}=\left(o+1\right)/12) 的零均值高斯函数：

[\beta{o}\left(x\right)\approx\frac{1}{\sqrt{2\pi\sigma_{o}{2}}}\exp\left(-\frac{x^{2}}{2\sigma_{o}}\right).]

可用于任意 (x) 和 (o) 计算该高斯函数的函数为 gauss_spline 。以下代码和图像使用样条滤波计算浣熊面部的边缘图像（平滑样条的二阶导数），浣熊面部是通过命令 scipy.datasets.face 返回的数组。命令 sepfir2d 用于应用具有镜像对称边界条件的可分离 2-D FIR 滤波器到样条系数。该函数非常适合从样条系数中重建样本，并且比 convolve2d 更快，后者可以卷积任意 2-D 滤波器并允许选择镜像对称边界条件。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt 

>>> image = datasets.face(gray=True).astype(np.float32)
>>> derfilt = np.array([1.0, -2, 1.0], dtype=np.float32)
>>> ck = signal.cspline2d(image, 8.0)
>>> deriv = (signal.sepfir2d(ck, derfilt, [1]) +
...          signal.sepfir2d(ck, [1], derfilt))

或者，我们可以这样做：

laplacian = np.array([[0,1,0], [1,-4,1], [0,1,0]], dtype=np.float32)
deriv2 = signal.convolve2d(ck,laplacian,mode='same',boundary='symm') 

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

"This code displays two plots. The first plot is a normal grayscale photo of a raccoon climbing on a palm plant. The second plot has the 2-D spline filter applied to the photo and is completely grey except the edges of the photo have been emphasized, especially on the raccoon fur and palm fronds."

>>> plt.figure()
>>> plt.imshow(deriv)
>>> plt.gray()
>>> plt.title('Output of spline edge filter')
>>> plt.show()

"This code displays two plots. The first plot is a normal grayscale photo of a raccoon climbing on a palm plant. The second plot has the 2-D spline filter applied to the photo and is completely grey except the edges of the photo have been emphasized, especially on the raccoon fur and palm fronds."

过滤

过滤是修改输入信号的任何系统的通用名称。在 SciPy 中，信号可以被看作是一个 NumPy 数组。有不同类型的滤波器用于不同类型的操作。有两种广泛的滤波操作：线性和非线性。线性滤波器总是可以简化为通过适当矩阵的乘积来实现另一个扁平化的 NumPy 数组。当然，这通常不是计算滤波器的最佳方式，因为涉及的矩阵和向量可能非常庞大。例如，使用这种方法对 (512 \times 512) 图像进行滤波将需要将 (512² \times 512²) 矩阵与 (512²) 向量相乘。仅尝试使用标准 NumPy 数组存储 (512² \times 512²) 矩阵将需要 (68,719,476,736) 个元素。每个元素占 4 字节，这将需要 (256\textrm{GB}) 内存。在大多数应用中，该矩阵的大多数元素为零，并且采用了不同的方法来计算滤波器的输出。

卷积/相关

许多线性滤波器还具有平移不变性的特性。这意味着滤波操作在信号的不同位置上是相同的，并且这意味着可以仅从矩阵的一行（或列）的知识构造滤波矩阵。在这种情况下，可以使用傅里叶变换来实现矩阵乘法。

让 (x\left[n\right]) 定义为以整数 (n) 为索引的一维信号。两个一维信号的完全卷积可以表示为

[y\left[n\right]=\sum_{k=-\infty}^{\infty}x\left[k\right]h\left[n-k\right].]

只有在我们将序列限制为可以存储在计算机中的有限支持序列时，才能直接实现此方程。选择 (n=0) 作为两个序列的起始点，让 (K+1) 是使得对所有 (n\geq K+1) 有 (x\left[n\right]=0) 的值，(M+1) 是使得对所有 (n\geq M+1) 有 (h\left[n\right]=0) 的值，那么离散卷积表达式为

[y\left[n\right]=\sum_{k=\max\left(n-M,0\right)}^{\min\left(n,K\right)}x\left[k\right]h\left[n-k\right].]

为了方便起见，假设 (K\geq M.) 那么，更明确地，此操作的输出是

\begin{eqnarray} y\left[0\right] & = & x\left[0\right]h\left[0\right]\ y\left[1\right] & = & x\left[0\right]h\left[1\right]+x\left[1\right]h\left[0\right]\ y\left[2\right] & = & x\left[0\right]h\left[2\right]+x\left[1\right]h\left[1\right]+x\left[2\right]h\left[0\right]\ \vdots & \vdots & \vdots\ y\left[M\right] & = & x\left[0\right]h\left[M\right]+x\left[1\right]h\left[M-1\right]+\cdots+x\left[M\right]h\left[0\right]\ y\left[M+1\right] & = & x\left[1\right]h\left[M\right]+x\left[2\right]h\left[M-1\right]+\cdots+x\left[M+1\right]h\left[0\right]\ \vdots & \vdots & \vdots\ y\left[K\right] & = & x\left[K-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[0\right]\ y\left[K+1\right] & = & x\left[K+1-M\right]h\left[M\right]+\cdots+x\left[K\right]h\left[1\right]\ \vdots & \vdots & \vdots\ y\left[K+M-1\right] & = & x\left[K-1\right]h\left[M\right]+x\left[K\right]h\left[M-1\right]\ y\left[K+M\right] & = & x\left[K\right]h\left[M\right].\end{eqnarray}

因此，两个长度分别为(K+1)和(M+1)的有限序列的完全离散卷积结果是一个长度为(K+M+1=\left(K+1\right)+\left(M+1\right)-1)的有限序列。

SciPy 中实现了 1-D 卷积，使用函数convolve。此函数接受信号(x)、(h)和两个可选标志‘mode’和‘method’作为输入，并返回信号(y)。

第一个可选标志‘mode’允许指定要返回的输出信号的哪个部分。默认值‘full’返回整个信号。如果标志的值为‘same’，则只返回长度与第一个输入相同的中间(K)个值，从(y\left[\left\lfloor \frac{M-1}{2}\right\rfloor \right])开始。如果标志的值为‘valid’，则只返回中间(K-M+1=\left(K+1\right)-\left(M+1\right)+1)个输出值，其中(z)取决于从(h\left[0\right])到(h\left[M\right])的最小输入的所有值。换句话说，只返回(y\left[M\right])到(y\left[K\right])（含）的值。

第二个可选标志‘method’确定如何计算卷积，可以通过傅里叶变换方法使用fftconvolve或直接方法。默认情况下，选择预期较快的方法。傅里叶变换方法的时间复杂度为(O(N\log N))，而直接方法的时间复杂度为(O(N²))。根据大 O 常数和(N)的值，其中一种方法可能更快。默认值‘auto’执行粗略计算并选择预期较快的方法，而值‘direct’和‘fft’则强制使用另外两种方法计算。

下面的代码展示了两个序列的卷积的简单示例：

>>> x = np.array([1.0, 2.0, 3.0])
>>> h = np.array([0.0, 1.0, 0.0, 0.0, 0.0])
>>> signal.convolve(x, h)
array([ 0.,  1.,  2.,  3.,  0.,  0.,  0.])
>>> signal.convolve(x, h, 'same')
array([ 2.,  3.,  0.])

这个相同的函数convolve实际上可以接受 N-D 数组作为输入，并返回两个数组的 N-D 卷积，如下面的代码示例所示。对于这种情况，相同的输入标志也是可用的。

>>> x = np.array([[1., 1., 0., 0.], [1., 1., 0., 0.], [0., 0., 0., 0.], [0., 0., 0., 0.]])
>>> h = np.array([[1., 0., 0., 0.], [0., 0., 0., 0.], [0., 0., 1., 0.], [0., 0., 0., 0.]])
>>> signal.convolve(x, h)
array([[ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
 [ 1.,  1.,  0.,  0.,  0.,  0.,  0.],
 [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
 [ 0.,  0.,  1.,  1.,  0.,  0.,  0.],
 [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
 [ 0.,  0.,  0.,  0.,  0.,  0.,  0.],
 [ 0.,  0.,  0.,  0.,  0.,  0.,  0.]])

相关性非常类似于卷积，只是减号变为加号。因此，

[w\left[n\right]=\sum_{k=-\infty}^{\infty}y\left[k\right]x\left[n+k\right],]

是信号(y)和(x)的(交)相关性。对于长度有限的信号，其中(y\left[n\right]=0)超出范围(\left[0,K\right])，(x\left[n\right]=0)超出范围(\left[0,M\right])，求和可以简化为

[w\left[n\right]=\sum_{k=\max\left(0,-n\right)}^{\min\left(K,M-n\right)}y\left[k\right]x\left[n+k\right].]

再次假设(K\geq M)，这是

\begin{eqnarray} w\left[-K\right] & = & y\left[K\right]x\left[0\right]\ w\left[-K+1\right] & = & y\left[K-1\right]x\left[0\right]+y\left[K\right]x\left[1\right]\ \vdots & \vdots & \vdots\ w\left[M-K\right] & = & y\left[K-M\right]x\left[0\right]+y\left[K-M+1\right]x\left[1\right]+\cdots+y\left[K\right]x\left[M\right]\ w\left[M-K+1\right] & = & y\left[K-M-1\right]x\left[0\right]+\cdots+y\left[K-1\right]x\left[M\right]\ \vdots & \vdots & \vdots\ w\left[-1\right] & = & y\left[1\right]x\left[0\right]+y\left[2\right]x\left[1\right]+\cdots+y\left[M+1\right]x\left[M\right]\ w\left[0\right] & = & y\left[0\right]x\left[0\right]+y\left[1\right]x\left[1\right]+\cdots+y\left[M\right]x\left[M\right]\ w\left[1\right] & = & y\left[0\right]x\left[1\right]+y\left[1\right]x\left[2\right]+\cdots+y\left[M-1\right]x\left[M\right]\ w\left[2\right] & = & y\left[0\right]x\left[2\right]+y\left[1\right]x\left[3\right]+\cdots+y\left[M-2\right]x\left[M\right]\ \vdots & \vdots & \vdots\ w\left[M-1\right] & = & y\left[0\right]x\left[M-1\right]+y\left[1\right]x\left[M\right]\ w\left[M\right] & = & y\left[0\right]x\left[M\right].\end{eqnarray}

SciPy 函数correlate实现了这个操作。同样的标志也适用于此操作，以返回完整的(K+M+1)长度序列（‘full’），或一个从(w\left[-K+\left\lfloor \frac{M-1}{2}\right\rfloor \right])开始大小相同的序列（‘same’），或一个值取决于最小序列所有值的序列（‘valid’）。最后一种选项返回从(w\left[M-K\right])到(w\left[0\right])（包括）的(K-M+1)个值。

函数correlate还可以接受任意的 N-D 数组作为输入，并返回两个数组的 N-D 卷积。

当 (N=2,) 可以使用correlate 和/或 convolve 构建任意图像滤波器，执行模糊化、增强和边缘检测等操作。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt 

>>> image = datasets.face(gray=True)
>>> w = np.zeros((50, 50))
>>> w[0][0] = 1.0
>>> w[49][25] = 1.0
>>> image_new = signal.fftconvolve(image, w) 

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

"此代码显示两个图。第一张图是熟悉的浣熊攀爬在棕榈树上的照片。第二张图应用了 FIR 滤波器，并且由于在滤波核定义中手动设置的双峰，照片被叠加显示了两份。"

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

"此代码显示两个图。第一张图是熟悉的浣熊攀爬在棕榈树上的照片。第二张图应用了 FIR 滤波器，并且由于在滤波核定义中手动设置的双峰，照片被叠加显示了两份。"

在时间域中计算卷积，通常用于当信号之一远小于另一个时进行滤波（ (K\gg M) ），否则提供由函数fftconvolve 提供的频域进行线性滤波计算更有效。默认情况下，convolve 使用choose_conv_method 估算最快的方法。

如果滤波函数 (w[n,m]) 可以按如下分解：

[h[n, m] = h_1[n] h_2[m],]

卷积可以通过函数sepfir2d 计算。例如，我们考虑一个高斯滤波器gaussian。

[h[n, m] \propto e^{-x²-y²} = e^{-x²} e^{-y²},]

通常用于模糊处理。

>>> import numpy as np
>>> from scipy import signal, datasets
>>> import matplotlib.pyplot as plt 

>>> image = np.asarray(datasets.ascent(), np.float64)
>>> w = signal.windows.gaussian(51, 10.0)
>>> image_new = signal.sepfir2d(image, w, w) 

>>> plt.figure()
>>> plt.imshow(image)
>>> plt.gray()
>>> plt.title('Original image')
>>> plt.show()

"此代码显示两个图。第一张图是从下方拍摄的两人攀爬木楼梯的灰度照片。第二张图应用了 2-D 高斯 FIR 窗口，看起来非常模糊。你仍然可以看出这是一张照片，但主题不太明确。"

>>> plt.figure()
>>> plt.imshow(image_new)
>>> plt.gray()
>>> plt.title('Filtered image')
>>> plt.show()

"此代码显示两个图。第一张图是从下方拍摄的两人攀爬木楼梯的灰度照片。第二张图应用了 2-D 高斯 FIR 窗口，看起来非常模糊。你仍然可以看出这是一张照片，但主题不太明确。"

差分方程滤波

一般的线性 1-D 滤波器类（包括卷积滤波器）由差分方程描述。

[\sum_{k=0}{N}a_{k}y\left[n-k\right]=\sum_{k=0}b_{k}x\left[n-k\right],]

其中 (x\left[n\right]) 是输入序列，(y\left[n\right]) 是输出序列。如果假定初始休息使得对于 (n<0)，(y\left[n\right]=0)，那么这种类型的滤波器可以使用卷积来实现。然而，如果 (a_{k}\neq0) 对于 (k\geq1)，则卷积滤波器序列 (h\left[n\right]) 可能是无限的。此外，这种一般的线性滤波器允许对 (y\left[n\right]) 的初始条件进行放置，以使得对于 (n<0) 的滤波器不能使用卷积来表示。

差分方程滤波器可以被视为在前一个值的基础上递归地找到 (y\left[n\right])。

[a_{0}y\left[n\right]=-a_{1}y\left[n-1\right]-\cdots-a_{N}y\left[n-N\right]+\cdots+b_{0}x\left[n\right]+\cdots+b_{M}x\left[n-M\right].]

通常情况下，选择 (a_{0}=1) 进行归一化。SciPy 中这种一般差分方程滤波器的实现比前述方程复杂一些。实现方式是只需要延迟一个信号。实际的实现方程为（假设 (a_{0}=1) ）：

\begin{eqnarray} y\left[n\right] & = & b_{0}x\left[n\right]+z_{0}\left[n-1\right]\ z_{0}\left[n\right] & = & b_{1}x\left[n\right]+z_{1}\left[n-1\right]-a_{1}y\left[n\right]\ z_{1}\left[n\right] & = & b_{2}x\left[n\right]+z_{2}\left[n-1\right]-a_{2}y\left[n\right]\ \vdots & \vdots & \vdots\ z_{K-2}\left[n\right] & = & b_{K-1}x\left[n\right]+z_{K-1}\left[n-1\right]-a_{K-1}y\left[n\right]\ z_{K-1}\left[n\right] & = & b_{K}x\left[n\right]-a_{K}y\left[n\right],\end{eqnarray}

这里 (K=\max\left(N,M\right))。注意，如果 (K>M) 则 (b_{K}=0)，如果 (K>N) 则 (a_{K}=0)。因此，输出在时间 (n) 只依赖于时间 (n) 的输入以及前一个时间的 (z_{0}) 的值。只要计算和存储每个时间步骤中的 (K) 值 (z_{0}\left[n-1\right]\ldots z_{K-1}\left[n-1\right])，就可以始终计算这些值。

在 SciPy 中，可以使用命令lfilter调用这种差分方程滤波器。该命令的输入为向量 (b)、向量 (a)、信号 (x)，并返回与上述方程计算出的长度相同的向量 (y)。如果 (x) 是 N-D，则沿提供的轴计算滤波器。如果需要，可以提供初始条件，以提供 (z_{0}\left[-1\right]) 到 (z_{K-1}\left[-1\right]) 的值，否则将假定它们全部为零。如果提供了初始条件，则中间变量的最终条件也会被返回。例如，这些条件可以用来在相同状态下重新启动计算。

有时，用信号(x\left[n\right])和(y\left[n\right])来表达初始条件更加方便。换句话说，也许你有(x\left[-M\right])到(x\left[-1\right])的值和(y\left[-N\right])到(y\left[-1\right])的值，希望确定应该作为差分方程滤波器初始条件的(z_{m}\left[-1\right])的值。不难证明，对于(0\leq m<K)，

[z_{m}\left[n\right]=\sum_{p=0}^{K-m-1}\left(b_{m+p+1}x\left[n-p\right]-a_{m+p+1}y\left[n-p\right]\right).]

使用这个公式，我们可以找到给定(y)（和(x))的初始条件向量(z_{0}\left[-1\right])到(z_{K-1}\left[-1\right])。命令lfiltic执行此功能。

例如，考虑以下系统：

[y[n] = \frac{1}{2} x[n] + \frac{1}{4} x[n-1] + \frac{1}{3} y[n-1]]

该代码计算给定信号(x[n])的信号(y[n])；首先对于初始条件(y[-1] = 0)（默认情况），然后通过lfiltic将其设置为(y[-1] = 2)。

>>> import numpy as np
>>> from scipy import signal 

>>> x = np.array([1., 0., 0., 0.])
>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.lfilter(b, a, x)
array([0.5, 0.41666667, 0.13888889, 0.0462963])
>>> zi = signal.lfiltic(b, a, y=[2.])
>>> signal.lfilter(b, a, x, zi=zi)
(array([ 1.16666667,  0.63888889,  0.21296296,  0.07098765]), array([0.02366]))

注意，输出信号(y[n])的长度与输入信号(x[n])的长度相同。

线性系统分析

描述线性差分方程的线性系统可以完全由上述系数向量(a)和(b)描述；另一种表示是通过其传递函数(H(z))提供因子(k)、(N_z)个零点(z_k)和(N_p)个极点(p_k)。

[H(z) = k \frac{ (z-z_1)(z-z_2)...(z-z_{N_z})}{ (z-p_1)(z-p_2)...(z-p_{N_p})}.]

这种替代表示可以通过 SciPy 函数tf2zpk获得；逆操作由zpk2tf提供。

对于上述例子，我们有

>>> b = np.array([1.0/2, 1.0/4])
>>> a = np.array([1.0, -1.0/3])
>>> signal.tf2zpk(b, a)
(array([-0.5]), array([ 0.33333333]), 0.5)

即，该系统在(z=-1/2)处有一个零点，在(z=1/3)处有一个极点。

SciPy 函数freqz允许计算由系数(a_k)和(b_k)描述的系统的频率响应。查看freqz函数的帮助以获取详细示例。

滤波器设计

时间离散滤波器可以分为有限响应（FIR）滤波器和无限响应（IIR）滤波器。FIR 滤波器可以提供线性相位响应，而 IIR 滤波器则不能。SciPy 提供了设计这两种类型滤波器的函数。

FIR 滤波器

函数 firwin 根据窗口方法设计滤波器。根据提供的参数，函数返回不同类型的滤波器（如低通、带通等）。

下面的示例分别设计了低通和带阻滤波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> b1 = signal.firwin(40, 0.5)
>>> b2 = signal.firwin(41, [0.3, 0.8])
>>> w1, h1 = signal.freqz(b1)
>>> w2, h2 = signal.freqz(b2) 

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w1, 20*np.log10(np.abs(h1)), 'b')
>>> plt.plot(w2, 20*np.log10(np.abs(h2)), 'r')
>>> plt.ylabel('Amplitude Response (dB)')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

"此代码显示了一个振幅响应在 Y 轴上 vs 频率在 X 轴上的 X-Y 图。第一个（低通）蓝色迹线从 0 dB 的通带开始，并在中途左右开始下降，阻带约下降 80 dB。第二个（带阻）红色迹线从 0 dB 开始并结束，但中间部分从峰值下降约 60 dB，滤波器会抑制信号时带有波动。"

注意 firwin 默认使用归一化频率，其中值 (1) 对应于 Nyquist 频率，而函数 freqz 定义为值 (\pi) 对应于 Nyquist 频率。

函数 firwin2 允许通过指定一组角频率和对应增益来设计几乎任意频率响应。

下面的示例设计了一个具有任意振幅响应的滤波器。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> b = signal.firwin2(150, [0.0, 0.3, 0.6, 1.0], [1.0, 2.0, 0.5, 0.0])
>>> w, h = signal.freqz(b) 

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, np.abs(h))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

"此代码显示了一个振幅响应在 Y 轴上 vs 频率在 X 轴上的 X-Y 图。单个迹线形状类似于心跳信号。"

请注意 firwin2 和 freqz 的 y 轴的线性缩放以及 Nyquist 频率的不同定义（如上所述）。

IIR 滤波器

SciPy 提供了两个函数来直接设计 IIR 滤波器 iirdesign 和 iirfilter，其中滤波器类型（如椭圆）作为参数传递，并提供了几个特定类型滤波器设计函数，例如 ellip。

下面的示例设计了一个椭圆低通滤波器，并分别定义了通带和阻带波动。请注意，为了达到相同的阻带衰减（约 60 dB），滤波器阶数（4 阶）比上面示例中的 FIR 滤波器低得多。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> b, a = signal.iirfilter(4, Wn=0.2, rp=5, rs=60, btype='lowpass', ftype='ellip')
>>> w, h = signal.freqz(b, a) 

>>> plt.title('Digital filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.title('Digital filter frequency response')
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency (rad/sample)')
>>> plt.grid()
>>> plt.show()

"此代码生成一个振幅响应在 X 轴上频率的 X-Y 图。单个迹线显示一个平滑的低通滤波器，左侧的第三个通带接近 0 dB。右侧的三分之二约为 60 dB，有两个锐利的狭窄谷，向下跌至-100 dB。"

滤波器系数

滤波器系数可以用几种不同的格式存储：

  • ‘ba’ 或 ‘tf’ = 传递函数系数

  • ‘zpk’ = 零点、极点和总增益

  • ‘ss’ = 状态空间系统表示

  • ‘sos’ = 二阶段传递函数系数

函数，如 tf2zpk 和 zpk2ss，可以在它们之间进行转换。

传递函数表示

ba 或 tf 格式是一个 2 元组 (b, a)，表示传递函数，其中 b 是长度为 M+1 的系数数组，表示 M 阶分子多项式，a 是长度为 N+1 的系数数组，表示 N 阶分母多项式，系数按传递函数变量的正降幂排列。因此，元组 (b = [b_0, b_1, ..., b_M]) 和 (a =[a_0, a_1, ..., a_N]) 可以表示形如下的模拟滤波器：

[H(s) = \frac {b_0 s^M + b_1 s^{(M-1)} + \cdots + b_M} {a_0 s^N + a_1 s^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i s^{(M-i)}} {\sum_{i=0}^N a_i s^{(N-i)}}]

或者离散时间滤波器的形式：

[H(z) = \frac {b_0 z^M + b_1 z^{(M-1)} + \cdots + b_M} {a_0 z^N + a_1 z^{(N-1)} + \cdots + a_N} = \frac {\sum_{i=0}^M b_i z^{(M-i)}} {\sum_{i=0}^N a_i z^{(N-i)}}.]

此“正幂”形式在控制工程中更常见。如果 M 和 N 相等（这对所有通过双线性变换生成的滤波器都是真的），则这等同于 DSP 中首选的“负幂”离散时间形式：

[H(z) = \frac {b_0 + b_1 z^{-1} + \cdots + b_M z^{-M}} {a_0 + a_1 z^{-1} + \cdots + a_N z^{-N}} = \frac {\sum_{i=0}^M b_i z^{-i}} {\sum_{i=0}^N a_i z^{-i}}.]

尽管这对于常见的滤波器是真的，但请记住，这在一般情况下并非如此。如果 M 和 N 不相等，则在找到极点和零点之前，必须先将离散时间传递函数系数转换为“正幂”形式。

此表示在高阶存在数值误差，因此在可能的情况下，更喜欢使用其他格式。

零点和极点表示

zpk 格式是一个 3 元组 (z, p, k)，其中 z 是长度为 M 的复零点数组 (z = [z_0, z_1, ..., z_{M-1}])，p 是长度为 N 的复极点数组 (p = [p_0, p_1, ..., p_{N-1}])，k 是标量增益。这些表示数字传递函数：

[H(z) = k \cdot \frac {(z - z_0) (z - z_1) \cdots (z - z_{(M-1)})} {(z - p_0) (z - p_1) \cdots (z - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (z - z_i)} {\prod_{i=0}^{N-1} (z - p_i)}]

或者模拟传递函数：

[H(s) = k \cdot \frac {(s - z_0) (s - z_1) \cdots (s - z_{(M-1)})} {(s - p_0) (s - p_1) \cdots (s - p_{(N-1)})} = k \frac {\prod_{i=0}^{M-1} (s - z_i)} {\prod_{i=0}^{N-1} (s - p_i)}.]

尽管根集存储为有序的 NumPy 数组，但它们的顺序并不重要： ([-1, -2], [-3, -4], 1) 和 ([-2, -1], [-4, -3], 1) 是相同的滤波器。

状态空间系统表示

ss 格式是一个四元组 (A, B, C, D)，表示一个 N 阶数字/离散时间系统的状态空间形式：

[\begin{split}\mathbf{x}[k+1] = A \mathbf{x}[k] + B \mathbf{u}[k]\ \mathbf{y}[k] = C \mathbf{x}[k] + D \mathbf{u}[k]\end{split}]

或者连续/模拟系统的形式：

[\begin{split}\dot{\mathbf{x}}(t) = A \mathbf{x}(t) + B \mathbf{u}(t)\ \mathbf{y}(t) = C \mathbf{x}(t) + D \mathbf{u}(t),\end{split}]

具有 P 输入，Q 输出和 N 状态变量的系统，其中：

  • x 是状态向量

  • y 是长度为 Q 的输出向量

  • u 是长度为 P 的输入向量

  • A 是状态矩阵，形状为 (N, N)

  • B 是形状为(N, P)的输入矩阵

  • C 是形状为(Q, N)的输出矩阵

  • D 是形状为(Q, P)的前馈或反馈矩阵。 （在系统没有直接前馈的情况下，D 中的所有值均为零。）

状态空间是最一般的表示形式，也是唯一允许多输入多输出（MIMO）系统的表示形式。 对于给定的传递函数，有多种状态空间表示。 具体来说，“可控规范形式”和“可观测规范形式”的系数与 tf 表示相同，因此会遭受相同的数值误差。

二阶段表示

sos 格式是形状为(n_sections, 6)的单个二维数组，表示串联的二阶传递函数序列，以最小的数值误差实现更高阶的滤波器。 每行对应于一个二阶 tf 表示，前三列提供分子系数，后三列提供分母系数：

[[b_0, b_1, b_2, a_0, a_1, a_2]]

系数通常是标准化的，使得 (a_0) 始终为 1。 在浮点计算中通常不重要，滤波器输出将相同，无论顺序如何。

滤波器转换

IIR 滤波器设计函数首先生成一个标准化截止频率为 1 rad/sec 的模拟低通滤波器原型。 然后使用以下替换将其转换为其他频率和带类型：

类型 转换
lp2lp (s \rightarrow \frac{s}{\omega_0})
lp2hp (s \rightarrow \frac{\omega_0}{s})
lp2bp (s \rightarrow \frac{s² + {\omega_0}²}{s \cdot \mathrm{BW}})
lp2bs (s \rightarrow \frac{s \cdot \mathrm{BW}}{s² + {\omega_0}²})

在这里，(\omega_0) 是新的截止频率或中心频率，而 (\mathrm{BW}) 是带宽。这些在对数频率轴上保持对称。

要将转换后的模拟滤波器转换为数字滤波器，使用 bilinear 变换，进行以下替换：

[s \rightarrow \frac{2}{T} \frac{z - 1}{z + 1},]

当中 T 是采样时间（采样频率的倒数）。

其他滤波器

信号处理包还提供了许多其他滤波器。

中值滤波器

当噪声明显非高斯或者希望保留边缘时，通常会应用中值滤波器。中值滤波器通过对感兴趣点周围的像素值进行排序来工作。这个邻域像素值列表的样本中位数被用作输出数组的值。如果邻域中有偶数个元素，则使用中间两个值的平均值作为中位数。一个适用于 N-D 数组的通用中值滤波器是 medfilt。一个专门针对 2-D 数组的版本是 medfilt2d。

顺序滤波器

中值滤波器是更一般的称为顺序滤波器的滤波器类的一个特定示例。要计算特定像素的输出，所有顺序滤波器使用周围区域内的数组值。这些数组值被排序，然后从中选择一个作为输出值。对于中值滤波器，数组值列表的样本中位数用作输出。一般顺序滤波器允许用户选择在排序值中哪一个作为输出。例如，可以选择在列表中选择最大值或最小值。顺序滤波器除了输入数组和区域掩码之外，还接受一个指定哪些邻居数组值在排序列表中应该作为输出的元素。执行顺序滤波器的命令是order_filter.

Wiener 滤波器

Wiener 滤波器是一种简单的去模糊滤波器，用于图像去噪。这不是通常描述的图像重建问题中的 Wiener 滤波器，而是一种简单的局部均值滤波器。设(x)为输入信号，则输出为

[\begin{split}y=\left{ \begin{array}{cc} \frac{\sigma{2}}{\sigma_{x}{2}}m_{x}+\left(1-\frac{\sigma{2}}{\sigma_{x}{2}}\right)x & \sigma_{x}{2}\geq\sigma,\ m_{x} & \sigma_{x}{2}<\sigma,\end{array}\right.\end{split}]

其中(m_{x})是均值的局部估计，(\sigma_{x}{2})是方差的局部估计。这些估计的窗口是一个可选的输入参数（默认为(3\times3)）。参数(\sigma)是一个阈值噪声参数。如果未提供(\sigma)，则将其估计为局部方差的平均值。

Hilbert 滤波器

Hilbert 变换从实信号构造复值解析信号。例如，如果(x=\cos\omega n)，那么(y=\textrm{hilbert}\left(x\right))会返回（除了接近边缘的地方）(y=\exp\left(j\omega n\right))。在频率域中，Hilbert 变换执行

[Y=X\cdot H,]

其中(H)对于正频率为(2)，对于负频率为(0)，对于零频率为(1)。

模拟滤波器设计

函数iirdesign、iirfilter以及特定滤波器类型的滤波器设计函数（例如ellip）都具有一个analog标志，允许设计模拟滤波器。

下面的示例设计了一个模拟（IIR）滤波器，通过tf2zpk获得极点和零点，并在复平面中绘制它们。在幅度响应中可以清楚地看到大约在(\omega \approx 150)和(\omega \approx 300)处的零点。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> b, a = signal.iirdesign(wp=100, ws=200, gpass=2.0, gstop=40., analog=True)
>>> w, h = signal.freqs(b, a) 

>>> plt.title('Analog filter frequency response')
>>> plt.plot(w, 20*np.log10(np.abs(h)))
>>> plt.ylabel('Amplitude Response [dB]')
>>> plt.xlabel('Frequency')
>>> plt.grid()
>>> plt.show()

这段代码显示了两个图。第一个图是 IIR 滤波器响应的 X-Y 图，其中 Y 轴是幅度响应，X 轴是频率。所示的低通滤波器在 0 到 100 Hz 的通带内具有 0 dB 响应，并且在约 175 Hz 到 1 KHz 之间的阻带下降约 40 dB。滤波器在约 175 Hz 和 300 Hz 附近有两个明显的不连续点。第二个图是在复平面上显示传递函数的 X-Y 图。Y 轴是实值，X 轴是复值。滤波器在[300+0j, 175+0j, -175+0j, -300+0j]附近有四个零点，显示为蓝色 X 标记。滤波器还在[50-30j, -50-30j, 100-8j, -100-8j]附近有四个极点，显示为红色点。

>>> z, p, k = signal.tf2zpk(b, a) 

>>> plt.plot(np.real(z), np.imag(z), 'ob', markerfacecolor='none')
>>> plt.plot(np.real(p), np.imag(p), 'xr')
>>> plt.legend(['Zeros', 'Poles'], loc=2) 

>>> plt.title('Pole / Zero Plot')
>>> plt.xlabel('Real')
>>> plt.ylabel('Imaginary')
>>> plt.grid()
>>> plt.show()

这段代码显示了两个图。第一个图是 IIR 滤波器响应的 X-Y 图，其中 Y 轴是幅度响应，X 轴是频率。所示的低通滤波器在 0 到 100 Hz 的通带内具有 0 dB 响应，并且在约 175 Hz 到 1 KHz 之间的阻带下降约 40 dB。滤波器在约 175 Hz 和 300 Hz 附近有两个明显的不连续点。第二个图是在复平面上显示传递函数的 X-Y 图。Y 轴是实值，X 轴是复值。滤波器在[300+0j, 175+0j, -175+0j, -300+0j]附近有四个零点，显示为蓝色 X 标记。滤波器还在[50-30j, -50-30j, 100-8j, -100-8j]附近有四个极点，显示为红色点。

频谱分析

周期图测量

scipy 函数periodogram提供了使用周期图方法估算频谱密度的方法。

下面的示例计算了白高斯噪声中正弦信号的周期图。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*np.sqrt(2)
>>> freq = 1270.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> x = amp*np.sin(2*np.pi*freq*time)
>>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

>>> f, Pper_spec = signal.periodogram(x, fs, 'flattop', scaling='spectrum') 

>>> plt.semilogy(f, Pper_spec)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD')
>>> plt.grid()
>>> plt.show()

这段代码显示了单个 X-Y 对数线性图，其中 Y 轴上的功率谱密度与 X 轴上的频率相对应。单条蓝色曲线显示了一个功率级别为 1e-3 的噪声底线，1270 Hz 处有一个功率为 1 的单峰。噪声底线的测量显示噪声波动下降到 1e-7。

使用 Welch 方法进行频谱分析

特别是在抗噪性方面改进的方法是 Welch 方法，由 scipy 函数welch实现。

下面的示例使用了Welch 方法来估算频谱，并使用与上述示例相同的参数。请注意频谱图的噪声底线要平滑得多。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> fs = 10e3
>>> N = 1e5
>>> amp = 2*np.sqrt(2)
>>> freq = 1270.0
>>> noise_power = 0.001 * fs / 2
>>> time = np.arange(N) / fs
>>> x = amp*np.sin(2*np.pi*freq*time)
>>> x += np.random.normal(scale=np.sqrt(noise_power), size=time.shape) 

>>> f, Pwelch_spec = signal.welch(x, fs, scaling='spectrum') 

>>> plt.semilogy(f, Pwelch_spec)
>>> plt.xlabel('frequency [Hz]')
>>> plt.ylabel('PSD')
>>> plt.grid()
>>> plt.show()

"此代码显示一个 X-Y 对数线性图，Y 轴上是功率谱密度，X 轴上是频率。单个蓝色追踪显示了一个平滑的噪声底部，功率水平为 6e-2，一个峰值功率高达 2，频率为 1270 Hz。"

朗伯-斯卡戈尔周期图 (lombscargle)

最小二乘谱分析（LSSA）[1] [2] 是一种基于对数据样本进行正弦拟合的最小二乘拟合来估计频谱的方法，类似于傅里叶分析。傅里叶分析是科学中使用最广泛的谱方法，通常会增加长周期噪声在长时间间隔记录中；LSSA 减少了这类问题。

朗伯-斯卡戈尔方法对不均匀采样数据执行频谱分析，已知是发现和检验弱周期信号的有效方法。

对于包含(N_{t})个测量(X_{j}\equiv X(t_{j}))的时间序列，采样时间为(t_{j})，假设已经进行了缩放和移位，使其平均值为零，方差为单位，频率为(f)的归一化朗伯-斯卡戈尔周期图为

[P_{n}(f) = \frac{1}{2}\left{\frac{\left[\sum_{j}{N_{t}}X_{j}\cos\omega(t_{j}-\tau)\right]{2}}{\sum_{j}{N_{t}}\cos\omega(t_{j}-\tau)}+\frac{\left[\sum_{j}{N_{t}}X_{j}\sin\omega(t_{j}-\tau)\right]{2}}{\sum_{j}{N_{t}}\sin\omega(t_{j}-\tau)}\right}.]

这里，(\omega \equiv 2\pi f) 是角频率。频率依赖的时间偏移 (\tau) 由以下公式给出

[\tan 2\omega\tau = \frac{\sum_{j}^{N_{t}}\sin 2\omega t_{j}}{\sum_{j}^{N_{t}}\cos 2\omega t_{j}}.]

lombscargle函数使用了由汤森德[3]提出的稍微修改的算法，允许在每个频率下通过输入数组的单次遍历计算周期图。

方程被重构为：

[P_{n}(f) = \frac{1}{2}\left[\frac{(c_{\tau}XC + s_{\tau}XS){2}}{c_{\tau}CC + 2c_{\tau}s_{\tau}CS + s_{\tau}^{2}SS} + \frac{(c_{\tau}XS - s_{\tau}XC){2}}{c_{\tau}SS - 2c_{\tau}s_{\tau}CS + s_{\tau}^{2}CC}\right]]

和

[\tan 2\omega\tau = \frac{2CS}{CC-SS}.]

这里，

[c_{\tau} = \cos\omega\tau,\qquad s_{\tau} = \sin\omega\tau,]

当求和时

[\begin{split}XC &= \sum_{j}^{N_{t}} X_{j}\cos\omega t_{j}\ XS &= \sum_{j}^{N_{t}} X_{j}\sin\omega t_{j}\ CC &= \sum_{j}^{N_{t}} \cos^{2}\omega t_{j}\ SS &= \sum_{j}^{N_{t}} \sin^{2}\omega t_{j}\ CS &= \sum_{j}^{N_{t}} \cos\omega t_{j}\sin\omega t_{j}.\end{split}]

需要进行(N_{f}(2N_{t}+3))次三角函数求值，相较于简单实现，速度提升约(\sim 2)倍。

[% LaTeX 宏使 LaTeX 公式更易读： \newcommand{\IC}{{\mathbb{C}}} % 复数集合 \newcommand{\IN}{{\mathbb{N}}} % 自然数集合 \newcommand{\IR}{{\mathbb{R}}} % 实数集合 \newcommand{\IZ}{{\mathbb{Z}}} % 整数集合 \newcommand{\jj}{{\mathbb{j}}} % 虚数单位 \newcommand{\e}{\operatorname{e}} % 自然对数的底 \newcommand{\dd}{\operatorname{d}} % 无穷小算子 \newcommand{\conj}[1]{\overline{#1}} % 复共轭 \newcommand{\conjT}[1]{\overline{#1^T}} % 转置复共轭 \newcommand{\inv}[1]{\left(#1\right)^{!-1}} % 逆 % 因为未加载物理包，我们自己定义宏： \newcommand{\vb}[1]{\mathbf{#1}} % 向量和矩阵加粗] ## 短时傅里叶变换

这一节提供了有关使用ShortTimeFFT类的一些背景信息：短时傅里叶变换（STFT）可用于分析信号随时间的频谱特性。它通过滑动窗口将信号分成重叠的块，并计算每个块的傅里叶变换。对于连续时间的复值信号 (x(t))，STFT 定义为[4]

[S(f, t) := \int_\IR x(\xi), \conj{w(\xi-t)},\e^{-\jj2\pi f \xi}\dd\xi\ ,]

其中 (w(t)) 是复值窗口函数，其复共轭为 (\conj{w(t)})。可以解释为确定 (x) 与窗口 (w) 的标量积，该窗口由时间 (t) 平移，然后被频率 (f) 调制（即频移）。对于采样信号 (x[k] := x(kT))，(k\in\IZ)，采样间隔 (T)（即采样频率 fs 的倒数），需要使用离散版本，即仅在离散网格点 (S[q, p] := S(q \Delta f, p\Delta t))，(q,p\in\IZ) 上进行评估。可以表达为

(1)[S[q,p] = \sum_{k=0}^{N-1} x[k],\conj{w[k-p h]}, \e^{-\jj2\pi q k / N}\ , \quad q,p\in\IZ\ ,]

其中 p 表示 (S) 的时间索引，时间间隔为 (\Delta t := h T)，(h\in\IN)（参见 delta_t），它可以表达为 (h) 个样本的 hop 大小。(q) 表示 (S) 的频率索引，步长为 (\Delta f := 1 / (N T))（参见 delta_f），使其与 FFT 兼容。(w[m] := w(mT))，(m\in\IZ) 是采样窗口函数。

要更符合ShortTimeFFT的实现，重构方程 (1) 为两个步骤是有意义的：

  1. 通过窗口 (w[m]) 提取由 (M) 个样本组成的第 (p) 个切片（参见m_num），其中心位于 (t[p] := p \Delta t = h T)（参见delta_t），即，

    (2)[x_p[m] = x!\big[m - \lfloor M/2\rfloor + h p\big], \conj{w[m]}\ , \quad m = 0, \ldots M-1\ ,]

    整数 (\lfloor M/2\rfloor) 表示 M//2，即窗口的中点（m_num_mid）。为了符号方便，假设对于 (k\not\in{0, 1, \ldots, N-1})，有 (x[k]:=0)。在子节滑动窗口中详细讨论了切片的索引。

  2. 然后执行离散傅立叶变换（即，一个 FFT）对 (x_p[m]) 进行操作。

    (3)[S[q, p] = \sum_{m=0}^{M-1} x_p[m] \exp!\big{% -2\jj\pi (q + \phi_m), m / M\big}\ .]

    注意可以指定线性相位 (\phi_m)（参见phase_shift），它对应于将输入向右移动 (\phi_m) 个样本。默认值为 (\phi_m = \lfloor M/2\rfloor)（与 phase_shift = 0 定义对应），这样可以抑制未移动信号的线性相位成分。此外，FFT 可能通过用零填充 (w[m]) 来进行过采样。这可以通过指定mfft大于窗口长度m_num来实现——这将把 (M) 设置为mfft（这意味着对于 (m\not\in{0, 1, \ldots, M-1})，也有 (w[m]:=0)）。

反短时傅立叶变换（istft）通过颠倒这两个步骤来实现：

  1. 执行逆离散傅立叶变换，即，

    [x_p[m] = \frac{1}{M}\sum_{q=0}^M S[q, p], \exp!\big{ 2\jj\pi (q + \phi_m), m / M\big}\ .]

  2. 加权求和偏移的切片，以重构原始信号，即，

    [x[k] = \sum_p x_p!\big[\mu_p(k)\big], w_d!\big[\mu_p(k)\big]\ ,\quad \mu_p(k) = k + \lfloor M/2\rfloor - h p]

    对于 (k \in [0, \ldots, n-1])。(w_d[m]) 是称为规范对偶窗口的 (w[m])，也由 (M) 个样本组成。

注意，并非所有窗口和跳跃大小都必然存在逆 STFT。对于给定的窗口 (w[m])，跳跃大小 (h) 必须足够小，以确保 (x[k]) 的每个样本都受到至少一个窗口片段的非零值影响。这有时被称为“非零重叠条件”（参见check_NOLA）。有关更多细节，请参阅小节 Inverse STFT and Dual Windows。

滑动窗口

本小节讨论了如何通过一个例子来索引ShortTimeFFT 中的滑动窗口：考虑一个长度为 6 的窗口，hop 间隔为 2，采样间隔为 T 为 1，例如 ShortTimeFFT (np.ones(6), 2, fs=1)。下图以示意图的方式展示了前四个窗口位置，也称为时间片段：

../_images/tutorial_stft_sliding_win_start.svg

x 轴表示时间 (t)，对应于底部蓝色框的样本索引 k。y 轴表示时间片段索引 (p)。信号 (x[k]) 从索引 (k=0) 开始，并以浅蓝色背景标出。按定义，零片段 ((p=0)) 居中于 (t=0)。每个片段的中心点 (m_num_mid)，这里是样本 6//2=3，用文本“mid”标记。默认情况下，stft 计算所有与信号有重叠的片段。因此，第一个片段在 p_min = -1，最低样本索引为 k_min = -5。第一个不被片段左边伸出影响的样本索引为 (p_{lb} = 2)，第一个不受边界影响的样本索引为 (k_{lb} = 5)。属性lower_border_end 返回元组 ((k_{lb}, p_{lb}))。

信号末端的行为如下所示，信号具有 (n=50) 个样本，如下所示，以蓝色背景标示：

../_images/tutorial_stft_sliding_win_stop.svg

这里最后一个片段索引为 (p=26)。因此，遵循 Python 约定的结束索引在范围之外，p_max = 27 表示第一个不接触信号的片段。相应的样本索引为 k_max = 55。第一个向右伸出的片段是 (p_{ub} = 24)，其第一个样本为 (k_{ub}=45)。函数upper_border_begin 返回元组 ((k_{ub}, p_{ub}))。 ### 逆 STFT 和双窗口

双窗口术语源于框架理论 [5]，在那里，框架是一种能够表示给定希尔伯特空间中任何函数的级数展开。在那里，展开式 ({g_k}) 和 ({h_k}) 如果对于给定希尔伯特空间 (\mathcal{H}) 中的所有函数 (f)，它们是双框架。

[f = \sum_{k\in\IN} \langle f, g_k\rangle h_k = \sum_{k\in\IN} \langle f, h_k\rangle g_k\ , \quad f \in \mathcal{H}\ ,

holds，其中 (\langle ., .\rangle) 表示 (\mathcal{H}) 的标量积。所有框架都有对偶框架 [5]。

在文献中，仅在离散网格点 (S(q \Delta f, p\Delta t)) 处评估的 STFT 被称为 “Gabor 框架” [4] [5]。由于窗口 (w[m]) 的支持仅限于有限区间，ShortTimeFFT 属于所谓的 “无痛非正交展开” 类别 [4]。在这种情况下，对偶窗口始终具有相同的支持，并且可以通过反转对角矩阵来计算。以下将粗略推导需要理解矩阵操作的过程：

由于 Eq. (1) 中给出的 STFT 是 (x[k]) 的线性映射，可以用向量-矩阵表示法表达。这使得通过线性最小二乘方法的形式解可以表达其逆（如 lstsq），从而导出了一个漂亮而简单的结果。

我们从重新表述 Eq. (2) 的窗函数开始。

(4)[\begin{split} \vb{x}p = \vb{W},\vb{x} = \begin{bmatrix} \cdots & 0 & w[0] & 0 & \cdots&&&\ & \cdots & 0 & w[1] & 0 & \cdots&&\ & & & & \ddots&&&\ &&\cdots & 0 & 0 & w[M-1] & 0 & \cdots \end{bmatrix}\begin{bmatrix} x[0]\ x[1]\ \vdots\ x[N-1] \end{bmatrix}\ ,\end{split}]

其中 (M \times N) 矩阵 (\vb{W}_{!p}) 在 ((ph)) 主对角线上具有非零条目，即，

(5)[\begin{split} W_p[m,k] = w[m], \delta_{m+ph,k}\ ,\quad \delta_{k,l} &= \begin{cases} 1 & \text{ for } k=l\ ,\ 0 & \text{ elsewhere ,} \end{cases}\end{split}]

其中 (\delta_{k,l}) 是 Kronecker Delta。Eq. (3) 可以表示为

[\vb{s}_p = \vb{F},\vb{x}_p \quad\text{with}\quad F[q,m] =\exp!\big{-2\jj\pi (q + \phi_m), m / M\big}\ ,]

允许将第 (p) 切片的 STFT 写为

(6)[\vb{s}p = \vb{F}\vb{W},\vb{x} =: \vb{G}_p,\vb{x} \quad\text{with}\quad s_p[q] = S[p,q]\ .

注意 (\vb{F}) 是酉的，即其逆等于其共轭转置，意味着 (\conjT{\vb{F}}\vb{F} = \vb{I})。

为了获得 STFT 的单个向量-矩阵方程，这些切片被堆叠成一个向量。

[\begin{split}\vb{s} := \begin{bmatrix} \vb{s}_0\ \vb{s}1\ \vdots\ \vb{s} \end{bmatrix} = \begin{bmatrix} \vb{G}_0\ \vb{G}1\ \vdots\ \vb{G} \end{bmatrix}, \vb{x} =: \vb{G}, \vb{x}\ ,

其中 (P) 是生成的 STFT 的列数。为了反演这个方程，可以利用 Moore-Penrose 逆 (\vb{G}^\dagger)。

(7)[\vb{x} = \inv{\conjT{\vb{G}}\vb{G}}, \conjT{\vb{G}} \vb{s} =: \vb{G}^\dagger \vb{s}\ ,]

如果存在

[\begin{split}\vb{D} := \conjT{\vb{G}}\vb{G} = \begin{bmatrix} \conjT{\vb{G}_0}& \conjT{\vb{G}1}& \cdots & \conjT{\vb{G}{P-1}} \end{bmatrix}^T \begin{bmatrix} \vb{G}0\ \vb{G}1\ \vdots\ \vb{G} \end{bmatrix} = \sum^{P-1} \conjT{\vb{G}_p}\vb{G}_p \ .\end{split}]

是可逆的。因此(\vb{x} = \vb{G}^\dagger\vb{G},\vb{x} = \inv{\conjT{\vb{G}}\vb{G}},\conjT{\vb{G}}\vb{G},\vb{x})显然成立。(\vb{D})始终是一个具有非负对角条目的对角矩阵。这在进一步简化(\vb{D})时变得清晰，

(8)[\vb{D} = \sum_{p=0}^{P-1} \conjT{\vb{G}p}\vb{G}p = \sum^{P-1} \conjT{(\vb{F},\vb{W})}, \vb{F},\vb{W}{!p} = \sum^{P-1} \conjT{\vb{W}{!p}}\vb{W} =: \sum_{p=0}^{P-1} \vb{D}_p]

由于(\vb{F})是酉矩阵。此外，

(9)[\begin{split}D_p[r,s] &= \sum_{m=0}^{M-1} \conj{W_p^T[r,m]},W_p[m,s] = \sum_{m=0}^{M-1} \left(\conj{w[m]}, \delta_{m+ph,r}\right) \left(w[m], \delta_{m+ph,s}\right)\ &= \sum_{m=0}^{M-1} \big|w[m]\big|², \delta_{r,s}, \delta_{r,m+ph}\ .\end{split}]

表明(\vb{D}_p)是一个具有非负实数条目的对角矩阵。因此，对(\vb{D}_p)求和保持这一性质。这允许进一步简化方程(7)，即

(10)[\begin{split}\vb{x} &= \vb{D}^{-1} \conjT{\vb{G}}\vb{s} = \sum_{p=0}^{P-1} \vb{D}^{-1}\conjT{\vb{W}{!p}}, \conjT{\vb{F}}\vb{s}p = \sum^{P-1} (\conj{\vb{W}\vb{D}{-1}})T, \conjT{\vb{F}}\vb{s}p\ &=: \sum^{P-1}\conjT{\vb{U}_p},\conjT{\vb{F}}\vb{s}_p\ .\end{split}]

利用方程(5)，(8)，(9)，(\vb{U}p=\vb{W}\vb{D}^{-1})可以表达为

(11)[\begin{split}U_p[m, k] &= W[m,k], D^{-1}[k,k] = \left(w[m] \delta_{m+ph,k}\right) \inv{\sum_{\eta=0}^{P-1} \vb{D}\eta[k,k]} \delta\ &= w[m] \inv{\sum_{\eta=0}{P-1}\sum_{\mu=0} \big|w[\mu]\big|²,\delta_{m+ph, \mu+\eta h}} \delta_{m+ph,k}\ &= w[m] \inv{\sum_{\eta=0}^{P-1} \big|w[m+(p-\eta)h]\big|²} \delta_{m+ph,k} \ .\end{split}]

这表明(\vb{U}_p)在方程(5)中与(\vb{W}_p)具有相同的结构，即只在(ph)次副对角线上有非零条目。逆中的求和项可以解释为将(w[\mu]²)滑动到(w[m])（带有反演），因此只有与(w[m])重叠的分量才会起作用。因此，所有远离边界的(U_p[m, k])都是相同的窗口。为了避免边界效应，(x[k])用零填充，扩大(\vb{U})，使所有接触到(x[k])的切片都包含相同的双窗口。

[w_d[m] = w[m] \inv{\sum_{\eta\in\IZ} \big|w[m + \eta, h]\big|²}\ .]

由于 (w[m] = 0) 对 (m \not\in{0, \ldots, M-1}) 成立，因此只需要对满足 (|\eta| < M/h) 的索引 (\eta) 进行求和。插入方程 (6) 到方程 (10)，可以证明名为双窗口的名称，

[\vb{x} = \sum_{p=0}^{P-1} \conjT{\vb{U}p},\conjT{\vb{F}}, \vb{F},\vb{W},\vb{x} = \left(\sum_{p=0}^{P-1} \conjT{\vb{U}p},\vb{W}\right)\vb{x}\ ,]

表明 (\vb{U}p) 和 (\vb{W}) 是可互换的。因此，(w_d[m]) 也是具有双窗口 (w[m]) 的有效窗口。注意 (w_d[m]) 不是唯一的双窗口，因为通常 (\vb{s}) 的条目比 (\vb{x}) 多。可以证明，(w_d[m]) 具有最小能量（或 (L_2) 范数）[4]，这就是被称为“规范双窗口”的原因。### 与旧版本实现的比较

函数 stft，istft 和 spectrogram 的实现早于 ShortTimeFFT。本节讨论了旧的“遗留”版本和新的 ShortTimeFFT 实现之间的关键区别。重写的主要动机是认识到集成 双窗口 不能在不破坏兼容性的情况下以合理的方式进行。这为重新思考代码结构和参数化提供了机会，从而使一些隐含的行为更加显式。

以下示例比较了具有负斜率的复值啁啾信号的两个 STFTs:

>>> import matplotlib.pyplot as plt
>>> import numpy as np
>>> from scipy.fft import fftshift
>>> from scipy.signal import stft, istft, spectrogram, ShortTimeFFT
...
>>> fs, N = 200, 1001  # 200 Hz sampling rate for 5 s signal
>>> t_z = np.arange(N) / fs  # time indexes for signal
>>> z = np.exp(2j*np.pi*70 * (t_z - 0.2*t_z**2))  # complex-valued chirp
...
>>> nperseg, noverlap = 50, 40
>>> win = ('gaussian', 1e-2 * fs)  # Gaussian with 0.01 s standard dev.
...
>>> # Legacy STFT:
>>> f0_u, t0, Sz0_u = stft(z, fs, win, nperseg, noverlap,
...                        return_onesided=False, scaling='spectrum')
>>> f0, Sz0 = fftshift(f0_u), fftshift(Sz0_u, axes=0)
...
>>> # New STFT:
>>> SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap,
...                                fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz1 = SFT.stft(z)
...
>>> # Plot results:
>>> fig1, axx = plt.subplots(2, 1, sharex='all', sharey='all',
...                          figsize=(6., 5.))  # enlarge figure a bit
>>> t_lo, t_hi, f_lo, f_hi = SFT.extent(N, center_bins=True)
>>> t_str0 = r"Legacy stft() produces $%d\times%d$ points" % Sz0.T.shape
>>> t_str1 = r"ShortTimeFFT produces $%d\times%d$ points" % Sz1.T.shape
>>> _ = axx[0].set(title=t_str0, xlim=(t_lo, t_hi), ylim=(f_lo, f_hi))
>>> _ = axx[1].set(title=t_str1, xlabel="Time $t$ in seconds " +
...                rf"($\Delta t= %g\,$s)" % SFT.delta_t)
...
>>> # Calculate extent of plot with centered bins since imshow
... # does not interpolate by default:
... dt2 = (nperseg-noverlap) / fs  / 2 # equals SFT.delta_t / 2
>>> df2 = fs / nperseg / 2 # equals SFT.delta_f / 2
>>> extent0 = (-dt2, t0[-1] + dt2, f0[0] - df2, f0[-1] - df2)
>>> extent1 = SFT.extent(N, center_bins=True)
...
>>> kw = dict(origin='lower', aspect='auto', cmap='viridis')
>>> im1a = axx[0].imshow(abs(Sz0), extent=extent0, **kw)
>>> im1b = axx[1].imshow(abs(Sz1), extent=extent1, **kw)
>>> fig1.colorbar(im1b, ax=axx, label="Magnitude $|S_z(t, f)|$")
>>> _ = fig1.supylabel(rf"Frequency $f$ in Hertz ($\Delta f = %g\,$Hz)" %
...                    SFT.delta_f)
>>> plt.show()

../_images/signal-10.png

ShortTimeFFT 生成比旧版本多 3 个时间切片是主要的不同点。正如在 滑动窗口 部分中所述，所有接触信号的切片都包含在新版本中。这具有将 STFT 切片和重新组装的优势，如 ShortTimeFFT 的代码示例所示。此外，使用所有接触切片使得 ISTFT 在某些位置为零的窗口情况下更加稳健。

注意具有相同时间戳的切片产生相等的结果（在数值精度上），即：

>>> np.allclose(Sz0, Sz1[:, 2:-1])
True

通常，这些额外的切片包含非零值。由于我们示例中的大重叠，它们非常小。例如：

>>> abs(Sz1[:, 1]).min(), abs(Sz1[:, 1]).max()
(6.925060911593139e-07, 8.00271269218721e-07)

ISTFT 可用于重建原始信号：

>>> t0_r, z0_r = istft(Sz0_u, fs, win, nperseg, noverlap,
...                    input_onesided=False, scaling='spectrum')
>>> z1_r = SFT.istft(Sz1, k1=N)
...
>>> len(z0_r), len(z)
(1010, 1001)
>>> np.allclose(z0_r[:N], z)
True
>>> np.allclose(z1_r, z)
True

请注意，旧实现返回的信号比原始信号更长。另一方面，新版 istft 允许明确指定重建信号的起始索引 k0 和结束索引 k1。旧实现中的长度差异是因为信号长度不是切片的倍数。

在这个示例中，新版和旧版之间的进一步差异包括：

  • 参数 fft_mode='centered' 确保在图中的双边 FFT 中，零频率垂直居中。在旧实现中，需要使用 fftshift。fft_mode='twosided' 产生与旧版本相同的行为。

  • 参数 phase_shift=None 确保两个版本的相位相同。ShortTimeFFT 的默认值 0 产生带有额外线性相位项的 STFT 切片。

谱图被定义为 STFT 的绝对平方 [4]。spectrogram 提供的是与该定义一致的内容，即：

>>> np.allclose(SFT.spectrogram(z), abs(Sz1)**2)
True

另一方面，旧版 spectrogram 提供了另一种带有不同信号边界处理的 STFT 实现。以下示例展示了如何使用 ShortTimeFFT 来获取与旧版 spectrogram 生成的相同 SFT：

>>> # Legacy spectrogram (detrending for complex signals not useful):
>>> f2_u, t2, Sz2_u = spectrogram(z, fs, win, nperseg, noverlap,
...                               detrend=None, return_onesided=False,
...                               scaling='spectrum', mode='complex')
>>> f2, Sz2 = fftshift(f2_u), fftshift(Sz2_u, axes=0)
...
>>> # New STFT:
... SFT = ShortTimeFFT.from_window(win, fs, nperseg, noverlap,
...                                fft_mode='centered',
...                                scale_to='magnitude', phase_shift=None)
>>> Sz3 = SFT.stft(z, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
>>> t3 = SFT.t(N, p0=0, p1=(N-noverlap)//SFT.hop, k_offset=nperseg//2)
...
>>> np.allclose(t2, t3)
True
>>> np.allclose(f2, SFT.f)
True
>>> np.allclose(Sz2, Sz3)
True

与其他 STFT 不同的是，时间切片不是从 0 开始，而是从 nperseg//2 开始，即：

>>> t2
array([0.125, 0.175, 0.225, 0.275, 0.325, 0.375, 0.425, 0.475, 0.525,
 ...
 4.625, 4.675, 4.725, 4.775, 4.825, 4.875])

此外，只返回不突出到右侧的切片，将最后一个切片居中于 4.875 秒，使其比默认的 stft 参数化更短。

使用 mode 参数，传统 频谱图 还可以返回 ‘角度’、‘相位’、‘PSD’ 或 ‘幅度’。传统 频谱图 的 缩放 行为并不简单，因为它取决于参数 mode、scaling 和 return_onesided。在 短时傅里叶变换 中，并没有所有组合的直接对应关系，因为它仅提供窗口的 ‘幅度’、‘PSD’ 或根本没有 缩放。以下表格展示了这些对应关系：

传统 频谱图 短时傅里叶变换
模式 缩放 return_onesided
--- --- ---
PSD 密度 真实
PSD 密度 错误
幅度 谱 真实
幅度 谱 错误
复数 谱 真实
复数 谱 错误
PSD 谱 真实
PSD 谱 错误
复数 密度 真实
复数 密度 错误
幅度 密度 真实
幅度 密度 错误
— — —

当在复值输入信号上使用onesided输出时，旧的spectrogram切换到two-sided模式。ShortTimeFFT引发了一个TypeError，因为所使用的rfft函数只接受实值输入。这个评论讨论了旧的spectrogram参数在单个余弦输入时的变化。

去趋势

SciPy 提供了函数detrend，用于去除数据系列中的常数或线性趋势，以查看更高阶效果。

下面的示例去除了二阶多项式时间序列的常数和线性趋势，并绘制了剩余的信号分量。

>>> import numpy as np
>>> import scipy.signal as signal
>>> import matplotlib.pyplot as plt 

>>> t = np.linspace(-10, 10, 20)
>>> y = 1 + t + 0.01*t**2
>>> yconst = signal.detrend(y, type='constant')
>>> ylin = signal.detrend(y, type='linear') 

>>> plt.plot(t, y, '-rx')
>>> plt.plot(t, yconst, '-bo')
>>> plt.plot(t, ylin, '-k+')
>>> plt.grid()
>>> plt.legend(['signal', 'const. detrend', 'linear detrend'])
>>> plt.show()

“此代码生成一个没有单位的 X-Y 图。红色轨迹对应于原始信号曲线，从左下到右上。蓝色轨迹具有常数去趋势，并位于红色轨迹下方，Y 轴偏移为零。最后的黑色轨迹具有线性去趋势，并且从左到右几乎平坦，突出显示原始信号的曲线。这最后的轨迹具有平均斜率为零，并且看起来非常不同。”

参考资料

一些进一步阅读和相关软件：

线性代数（scipy.linalg）

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

当使用优化过的 ATLAS LAPACK 和 BLAS 库构建 SciPy 时，它具有非常快的线性代数功能。如果你深入挖掘，所有原始的 LAPACK 和 BLAS 库都可以用于更快的计算。在本节中，描述了这些例程的一些更易于使用的接口。

所有这些线性代数例程都期望一个可以转换为二维数组的对象。这些例程的输出也是一个二维数组。

scipy.linalg vs numpy.linalg

scipy.linalg 包含了所有 numpy.linalg 中的函数，还包括一些在 numpy.linalg 中不包含的更高级的函数。

使用 scipy.linalg 而不是 numpy.linalg 的另一个优势是，它始终编译了 BLAS/LAPACK 支持，而对于 NumPy 这是可选的。因此，根据 NumPy 的安装方式，SciPy 版本可能更快。

因此，除非你不希望将 scipy 作为 numpy 程序的依赖项，否则请使用 scipy.linalg 而不是 numpy.linalg。

numpy.matrix vs 2-D numpy.ndarray

表示矩阵的类和基本操作，如矩阵乘法和转置，是 numpy 的一部分。为了方便起见，我们在这里总结了 numpy.matrix 和 numpy.ndarray 之间的区别。

numpy.matrix 是一个矩阵类，其接口比 numpy.ndarray 更方便用于矩阵操作。例如，这个类支持类似 MATLAB 的通过分号创建语法，* 运算符默认进行矩阵乘法，并包含 I 和 T 成员，分别用作求逆和转置的快捷方式：

>>> import numpy as np
>>> A = np.asmatrix('[1 2;3 4]')
>>> A
matrix([[1, 2],
 [3, 4]])
>>> A.I
matrix([[-2\. ,  1\. ],
 [ 1.5, -0.5]])
>>> b = np.asmatrix('[5 6]')
>>> b
matrix([[5, 6]])
>>> b.T
matrix([[5],
 [6]])
>>> A*b.T
matrix([[17],
 [39]])

尽管 numpy.matrix 类很方便，但不建议使用，因为它并没有添加任何不能通过二维 numpy.ndarray 对象完成的功能，而且可能导致使用的类混淆。例如，上述代码可以重写为：

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
 [3, 4]])
>>> linalg.inv(A)
array([[-2\. ,  1\. ],
 [ 1.5, -0.5]])
>>> b = np.array([[5,6]]) #2D array
>>> b
array([[5, 6]])
>>> b.T
array([[5],
 [6]])
>>> A*b #not matrix multiplication!
array([[ 5, 12],
 [15, 24]])
>>> A.dot(b.T) #matrix multiplication
array([[17],
 [39]])
>>> b = np.array([5,6]) #1D array
>>> b
array([5, 6])
>>> b.T  #not matrix transpose!
array([5, 6])
>>> A.dot(b)  #does not matter for multiplication
array([17, 39])

scipy.linalg 的操作可以同样应用于 numpy.matrix 或二维 numpy.ndarray 对象。

基本例程

求逆

矩阵 (\mathbf{A}) 的逆矩阵是矩阵 (\mathbf{B})，使得 (\mathbf{AB}=\mathbf{I})，其中 (\mathbf{I}) 是由主对角线上的元素为 1 组成的单位矩阵。通常，(\mathbf{B}) 被表示为 (\mathbf{B}=\mathbf{A}^{-1})。在 SciPy 中，可以使用 linalg.inv (A) 获得 NumPy 数组 A 的逆矩阵，或者如果 A 是一个矩阵，则使用 A.I。例如，令

[\begin{split}\mathbf{A} = \left[\begin{array}{ccc} 1 & 3 & 5\ 2 & 5 & 1\ 2 & 3 & 8\end{array}\right],\end{split}]

然后

[\begin{split}\mathbf{A^{-1}} = \frac{1}{25} \left[\begin{array}{ccc} -37 & 9 & 22 \ 14 & 2 & -9 \ 4 & -3 & 1 \end{array}\right] = % \left[\begin{array}{ccc} -1.48 & 0.36 & 0.88 \ 0.56 & 0.08 & -0.36 \ 0.16 & -0.12 & 0.04 \end{array}\right].\end{split}]

以下示例展示了 SciPy 中的计算过程

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1,3,5],[2,5,1],[2,3,8]])
>>> A
array([[1, 3, 5],
 [2, 5, 1],
 [2, 3, 8]])
>>> linalg.inv(A)
array([[-1.48,  0.36,  0.88],
 [ 0.56,  0.08, -0.36],
 [ 0.16, -0.12,  0.04]])
>>> A.dot(linalg.inv(A)) #double check
array([[  1.00000000e+00,  -1.11022302e-16,  -5.55111512e-17],
 [  3.05311332e-16,   1.00000000e+00,   1.87350135e-16],
 [  2.22044605e-16,  -1.11022302e-16,   1.00000000e+00]])

解线性系统

使用 SciPy 命令linalg.solve可以直接解线性方程组。该命令需要输入一个矩阵和一个右手边向量，然后计算解向量。提供了一个输入对称矩阵的选项，适用时可以加快处理速度。例如，假设要解以下同时方程组：

\begin{eqnarray} x + 3y + 5z & = & 10 \ 2x + 5y + z & = & 8 \ 2x + 3y + 8z & = & 3 \end{eqnarray}

我们可以通过矩阵求逆找到解向量：

[\begin{split}\left[\begin{array}{c} x\ y\ z\end{array}\right]=\left[\begin{array}{ccc} 1 & 3 & 5\ 2 & 5 & 1\ 2 & 3 & 8\end{array}\right]^{-1}\left[\begin{array}{c} 10\ 8\ 3\end{array}\right]=\frac{1}{25}\left[\begin{array}{c} -232\ 129\ 19\end{array}\right]=\left[\begin{array}{c} -9.28\ 5.16\ 0.76\end{array}\right].\end{split}]

然而，最好使用 linalg.solve 命令，这样可以更快速和更稳定地数值计算。在这种情况下，如下示例所示，它将给出相同的答案：

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1, 2], [3, 4]])
>>> A
array([[1, 2],
 [3, 4]])
>>> b = np.array([[5], [6]])
>>> b
array([[5],
 [6]])
>>> linalg.inv(A).dot(b)  # slow
array([[-4\. ],
 [ 4.5]])
>>> A.dot(linalg.inv(A).dot(b)) - b  # check
array([[  8.88178420e-16],
 [  2.66453526e-15]])
>>> np.linalg.solve(A, b)  # fast
array([[-4\. ],
 [ 4.5]])
>>> A.dot(np.linalg.solve(A, b)) - b  # check
array([[ 0.],
 [ 0.]])

求行列式

方阵 (\mathbf{A}) 的行列式通常表示为 (\left|\mathbf{A}\right|)，在线性代数中经常使用。假设 (a_{ij}) 是矩阵 (\mathbf{A}) 的元素，并且令 (M_{ij}=\left|\mathbf{A}_{ij}\right|) 是通过从 (\mathbf{A}) 中去除第 (i^{\textrm{th}}) 行和第 (j^{\textrm{th}}) 列后留下的矩阵的行列式。那么，对于任意行 (i),

[\left|\mathbf{A}\right|=\sum_{j}\left(-1\right)^{i+j}a_{ij}M_{ij}.]

这是定义行列式的递归方式，其中基本情况是接受 (1\times1) 矩阵的行列式仅为矩阵元素。在 SciPy 中，可以使用linalg.det计算行列式。例如，矩阵的行列式为

[\begin{split}\mathbf{A=}\left[\begin{array}{ccc} 1 & 3 & 5\ 2 & 5 & 1\ 2 & 3 & 8\end{array}\right]\end{split}]

是

\begin{eqnarray} \left|\mathbf{A}\right| & = & 1\left|\begin{array}{cc} 5 & 1\ 3 & 8\end{array}\right|-3\left|\begin{array}{cc} 2 & 1\ 2 & 8\end{array}\right|+5\left|\begin{array}{cc} 2 & 5\ 2 & 3\end{array}\right|\ & = & 1\left(5\cdot8-3\cdot1\right)-3\left(2\cdot8-2\cdot1\right)+5\left(2\cdot3-2\cdot5\right)=-25.\end{eqnarray}.

在 SciPy 中，计算如下所示：

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
 [3, 4]])
>>> linalg.det(A)
-2.0

计算范数

也可以使用 SciPy 计算矩阵和向量的范数。使用linalg.norm函数可以根据 order 参数（默认为 2）计算所需阶数的向量或矩阵范数。这个函数接受一个秩为 1（向量）或秩为 2（矩阵）的数组以及一个可选的 order 参数。

对于向量x，参数可以是任意实数，包括inf或-inf。计算得到的范数为

[\begin{split}\left\Vert \mathbf{x}\right\Vert =\left{ \begin{array}{cc} \max\left|x_{i}\right| & \textrm{ord}=\textrm{inf}\ \min\left|x_{i}\right| & \textrm{ord}=-\textrm{inf}\ \left(\sum_{i}\left|x_{i}\right|{\textrm{ord}}\right){1/\textrm{ord}} & \left|\textrm{ord}\right|<\infty.\end{array}\right.\end{split}]

对于矩阵(\mathbf{A})，范数的有效值只有(\pm2,\pm1,) (\pm) inf 和‘fro’（或‘f’）。因此，

[\begin{split}\left\Vert \mathbf{A}\right\Vert =\left{ \begin{array}{cc} \max_{i}\sum_{j}\left|a_{ij}\right| & \textrm{ord}=\textrm{inf}\ \min_{i}\sum_{j}\left|a_{ij}\right| & \textrm{ord}=-\textrm{inf}\ \max_{j}\sum_{i}\left|a_{ij}\right| & \textrm{ord}=1\ \min_{j}\sum_{i}\left|a_{ij}\right| & \textrm{ord}=-1\ \max\sigma_{i} & \textrm{ord}=2\ \min\sigma_{i} & \textrm{ord}=-2\ \sqrt{\textrm{trace}\left(\mathbf{A}^{H}\mathbf{A}\right)} & \textrm{ord}=\textrm{'fro'}\end{array}\right.\end{split}]

其中(\sigma_{i})是矩阵(\mathbf{A})的奇异值。

例子：

>>> import numpy as np
>>> from scipy import linalg
>>> A=np.array([[1,2],[3,4]])
>>> A
array([[1, 2],
 [3, 4]])
>>> linalg.norm(A)
5.4772255750516612
>>> linalg.norm(A,'fro') # frobenius norm is the default
5.4772255750516612
>>> linalg.norm(A,1) # L1 norm (max column sum)
6
>>> linalg.norm(A,-1)
4
>>> linalg.norm(A,np.inf) # L inf norm (max row sum)
7

解决线性最小二乘问题和伪逆

线性最小二乘问题在应用数学的许多分支中出现。在这个问题中，寻找一组线性缩放系数，使得模型能够拟合数据。特别地，假设数据(y_{i})通过系数(c_{j})和模型函数(f_{j}\left(\mathbf{x}{i}\right))与数据(\mathbf{x})相关联，模型表达式为

[y_{i}=\sum_{j}c_{j}f_{j}\left(\mathbf{x}{i}\right)+\epsilon,]

其中(\epsilon_{i})表示数据中的不确定性。最小二乘策略是选择系数(c_{j})使得

[J\left(\mathbf{c}\right)=\sum_{i}\left|y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right|^{2}.]

理论上，全局最小值发生在

[\frac{\partial J}{\partial c_{n}{*}}=0=\sum_{i}\left(y_{i}-\sum_{j}c_{j}f_{j}\left(x_{i}\right)\right)\left(-f_{n}\left(x_{i}\right)\right)]

或

\begin{eqnarray} \sum_{j}c_{j}\sum_{i}f_{j}\left(x_{i}\right)f_{n}^{}\left(x_{i}\right) & = & \sum_{i}y_{i}f_{n}^{}\left(x_{i}\right)\ \mathbf{A}^{H}\mathbf{Ac} & = & \mathbf{A}^{H}\mathbf{y}\end{eqnarray},

其中

[\left{ \mathbf{A}\right} {ij}=f\left(x_{i}\right).]

当(\mathbf{A^{H}A})可逆时，

[\mathbf{c}=\left(\mathbf{A}{H}\mathbf{A}\right)\mathbf{A}{H}\mathbf{y}=\mathbf{A}\mathbf{y},]

其中(\mathbf{A}^{\dagger})被称为(\mathbf{A})的伪逆。请注意，使用此定义的(\mathbf{A})模型可以写成

[\mathbf{y}=\mathbf{Ac}+\boldsymbol{\epsilon}.]

命令linalg.lstsq将解决线性最小二乘问题，给定(\mathbf{A})和(\mathbf{y})。此外，linalg.pinv将给出(\mathbf{A}^{\dagger})，给定(\mathbf{A})。

下面的示例和图表演示了使用linalg.lstsq和linalg.pinv解决数据拟合问题。下面显示的数据是使用以下模型生成的:

[y_{i}=c_{1}e^{-x_{i}}+c_{2}x_{i},]

当(i=1\ldots10)时，(x_{i}=0.1i)，(c_{1}=5)，(c_{2}=4)。在(y_{i})和系数(c_{1})、(c_{2})上添加噪声，并使用线性最小二乘法估计。

>>> import numpy as np
>>> from scipy import linalg
>>> import matplotlib.pyplot as plt
>>> rng = np.random.default_rng() 

>>> c1, c2 = 5.0, 2.0
>>> i = np.r_[1:11]
>>> xi = 0.1*i
>>> yi = c1*np.exp(-xi) + c2*xi
>>> zi = yi + 0.05 * np.max(yi) * rng.standard_normal(len(yi)) 

>>> A = np.c_[np.exp(-xi)[:, np.newaxis], xi[:, np.newaxis]]
>>> c, resid, rank, sigma = linalg.lstsq(A, zi) 

>>> xi2 = np.r_[0.1:1.0:100j]
>>> yi2 = c[0]*np.exp(-xi2) + c[1]*xi2 

>>> plt.plot(xi,zi,'x',xi2,yi2)
>>> plt.axis([0,1.1,3.0,5.5])
>>> plt.xlabel('$x_i$')
>>> plt.title('Data fitting with linalg.lstsq')
>>> plt.show()

" "

广义逆

使用命令linalg.pinv计算广义逆。设(\mathbf{A})为一个(M\times N)矩阵，若(M>N)，则广义逆为

[\mathbf{A}{\dagger}=\left(\mathbf{A}\mathbf{A}\right){-1}\mathbf{A},]

如果(M<N)，则矩阵的广义逆是

[\mathbf{A}{#}=\mathbf{A}\left(\mathbf{A}\mathbf{A}{H}\right).]

如果(M=N)，则

[\mathbf{A}{\dagger}=\mathbf{A}{#}=\mathbf{A}^{-1},]

只要(\mathbf{A})是可逆的。

分解

在许多应用中，使用其他表示法对矩阵进行分解是有用的。SciPy 支持几种分解方法。

特征值和特征向量

特征值-特征向量问题是最常用的线性代数操作之一。在一个流行的形式中，特征值-特征向量问题是为某个方阵(\mathbf{A})找到标量(\lambda)和相应向量(\mathbf{v})，使得

[\mathbf{Av}=\lambda\mathbf{v}.]

对于一个(N\times N)矩阵，有(N)（不一定是不同的）特征值 — 是多项式的根

[\left|\mathbf{A}-\lambda\mathbf{I}\right|=0.]

特征向量(\mathbf{v})有时也被称为右特征向量，以区别于满足以下条件的另一组左特征向量

[\mathbf{v}_{L}{H}\mathbf{A}=\lambda\mathbf{v}_{L}]

或者

[\mathbf{A}{H}\mathbf{v}_{L}=\lambda\mathbf{v}_{L}.]

默认情况下，命令 linalg.eig 返回 (\lambda) 和 (\mathbf{v}.) 然而，它也可以只返回 (\mathbf{v}_{L}) 和 (\lambda) 本身（ linalg.eigvals 也只返回 (\lambda)）。

此外，linalg.eig 也可以解决更一般的特征值问题。

\begin{eqnarray} \mathbf{Av} & = & \lambda\mathbf{Bv}\ \mathbf{A}^{H}\mathbf{v}{L} & = & \lambda{*}\mathbf{B}\mathbf{v}\end{eqnarray}

对于方阵 (\mathbf{A}) 和 (\mathbf{B}.) 标准特征值问题是 (\mathbf{B}=\mathbf{I}.) 的一个示例。当可以解决广义特征值问题时，它提供了 (\mathbf{A}) 的分解如下：

[\mathbf{A}=\mathbf{BV}\boldsymbol{\Lambda}\mathbf{V}^{-1},]

其中 (\mathbf{V}) 是特征向量的集合列，(\boldsymbol{\Lambda}) 是特征值的对角线矩阵。

按定义，特征向量仅在一个常数比例因子内定义。在 SciPy 中，特征向量的缩放因子被选择为 (\left\Vert \mathbf{v}\right\Vert {2}=\sum_{i}v_{i}=1.)

例如，考虑找到矩阵的特征值和特征向量

[\begin{split}\mathbf{A}=\left[\begin{array}{ccc} 1 & 5 & 2\ 2 & 4 & 1\ 3 & 6 & 2\end{array}\right].\end{split}]

特征多项式是

\begin{eqnarray} \left|\mathbf{A}-\lambda\mathbf{I}\right| & = & \left(1-\lambda\right)\left[\left(4-\lambda\right)\left(2-\lambda\right)-6\right]-\ & & 5\left[2\left(2-\lambda\right)-3\right]+2\left[12-3\left(4-\lambda\right)\right]\ & = & -\lambda{3}+7\lambda+8\lambda-3.\end{eqnarray}

这个多项式的根是 (\mathbf{A}) 的特征值：

\begin{eqnarray} \lambda_{1} & = & 7.9579\ \lambda_{2} & = & -1.2577\ \lambda_{3} & = & 0.2997.\end{eqnarray}

与每个特征值对应的特征向量可以使用原方程找到。然后可以找到与这些特征值相关的特征向量。

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1, 2], [3, 4]])
>>> la, v = linalg.eig(A)
>>> l1, l2 = la
>>> print(l1, l2)   # eigenvalues
(-0.3722813232690143+0j) (5.372281323269014+0j)
>>> print(v[:, 0])   # first eigenvector
[-0.82456484  0.56576746]
>>> print(v[:, 1])   # second eigenvector
[-0.41597356 -0.90937671]
>>> print(np.sum(abs(v**2), axis=0))  # eigenvectors are unitary
[1\. 1.]
>>> v1 = np.array(v[:, 0]).T
>>> print(linalg.norm(A.dot(v1) - l1*v1))  # check the computation
3.23682852457e-16

奇异值分解

奇异值分解（SVD）可以看作是对不是方阵的矩阵的特征值问题的扩展。设(\mathbf{A})是一个(M\times N)矩阵，其中(M)和(N)是任意的。矩阵(\mathbf{A}{H}\mathbf{A})和(\mathbf{A}\mathbf{A})是大小分别为(N\times N)和(M\times M)的方正定矩阵。已知方正定矩阵的特征值是实数且非负。此外，(\mathbf{A}{H}\mathbf{A})和(\mathbf{A}\mathbf{A})最多有(\min\left(M,N\right))个相同的非零特征值。将这些正特征值定义为(\sigma_{i}{2})。这些的平方根称为(\mathbf{A})的奇异值。(\mathbf{A}\mathbf{A})的特征向量按列收集在一个(N\times N)的酉矩阵(\mathbf{V})中，而(\mathbf{A}\mathbf{A}^{H})的特征向量按列收集在酉矩阵(\mathbf{U})中，奇异值按主对角线项设置为奇异值的(M\times N)零矩阵(\mathbf{\boldsymbol{\Sigma}})。然后

[\mathbf{A=U}\boldsymbol{\Sigma}\mathbf{V}^{H}]

是矩阵(\mathbf{A})的奇异值分解。每个矩阵都有一个奇异值分解。有时，奇异值被称为(\mathbf{A})的频谱。使用命令linalg.svd将返回(\mathbf{U})，(\mathbf{V}^{H})和(\sigma_{i})作为奇异值数组。要获取矩阵(\boldsymbol{\Sigma})，请使用linalg.diagsvd。以下示例说明了linalg.svd的使用：

>>> import numpy as np
>>> from scipy import linalg
>>> A = np.array([[1,2,3],[4,5,6]])
>>> A
array([[1, 2, 3],
 [4, 5, 6]])
>>> M,N = A.shape
>>> U,s,Vh = linalg.svd(A)
>>> Sig = linalg.diagsvd(s,M,N)
>>> U, Vh = U, Vh
>>> U
array([[-0.3863177 , -0.92236578],
 [-0.92236578,  0.3863177 ]])
>>> Sig
array([[ 9.508032  ,  0\.        ,  0\.        ],
 [ 0\.        ,  0.77286964,  0\.        ]])
>>> Vh
array([[-0.42866713, -0.56630692, -0.7039467 ],
 [ 0.80596391,  0.11238241, -0.58119908],
 [ 0.40824829, -0.81649658,  0.40824829]])
>>> U.dot(Sig.dot(Vh)) #check computation
array([[ 1.,  2.,  3.],
 [ 4.,  5.,  6.]])

LU 分解

LU 分解将(M\times N)矩阵(\mathbf{A})表示为

[\mathbf{A}=\mathbf{P},\mathbf{L},\mathbf{U},]

其中(\mathbf{P})是一个(M\times M)的置换矩阵（单位矩阵的行置换），(\mathbf{L})是一个(M\times K)的下三角或梯形矩阵（(K=\min\left(M,N\right))），具有单位对角线，(\mathbf{U})是一个上三角或梯形矩阵。SciPy 中进行这种分解的命令是linalg.lu。

这样的分解通常用于解决许多同时方程组，其中左边不变，但右边会发生变化。例如，假设我们要解决

[\mathbf{A}\mathbf{x}{i}=\mathbf{b}]

对于许多不同的(\mathbf{b}_{i})，LU 分解允许将其写为

[\mathbf{PLUx}{i}=\mathbf{b}.]

因为 (\mathbf{L}) 是下三角形矩阵，这个方程可以通过前向和后向替换非常快速地解出 (\mathbf{U}\mathbf{x}_{i})，最终，通过 linalg.lu_factor 命令求解每个新右手边的系统。

Cholesky 分解

Cholesky 分解是 LU 分解适用于 Hermitian 正定矩阵的特例。当 (\mathbf{A}=\mathbf{A}^{H}) 并且对于所有的 (\mathbf{x})，(\mathbf{x}^{H}\mathbf{Ax}\geq0)，则可以找到 (\mathbf{A}) 的分解，使得

\begin{eqnarray} \mathbf{A} & = & \mathbf{U}^{H}\mathbf{U}\ \mathbf{A} & = & \mathbf{L}\mathbf{L}^{H}\end{eqnarray}，

其中 (\mathbf{L}) 是下三角形矩阵，(\mathbf{U}) 是上三角形矩阵。注意，(\mathbf{L}=\mathbf{U}^{H})。Cholesky 分解的命令 linalg.cholesky 计算 Cholesky 因子分解。对于使用 Cholesky 因子分解来解方程组，还有 linalg.cho_factor 和 linalg.cho_solve 类似于它们的 LU 分解对应的例程。

QR 分解

QR 分解（有时称为极分解）适用于任何 (M\times N) 数组，并找到一个 (M\times M) 的酉矩阵 (\mathbf{Q}) 和一个 (M\times N) 的上梯形矩阵 (\mathbf{R})，使得

[\mathbf{A=QR}]。

注意，如果已知 (\mathbf{A}) 的 SVD 分解，则可以找到 QR 分解。

[\mathbf{A}=\mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^{H}=\mathbf{QR}]。

意味着 (\mathbf{Q}=\mathbf{U}) 和 (\mathbf{R}=\boldsymbol{\Sigma}\mathbf{V}^{H})。然而，在 SciPy 中使用独立的算法来找到 QR 和 SVD 分解。QR 分解的命令是 linalg.qr。

Schur 分解

对于一个方阵 (N\times N)，(\mathbf{A})，Schur 分解找到（不一定唯一的）矩阵 (\mathbf{T}) 和 (\mathbf{Z})，使得

[\mathbf{A}=\mathbf{ZT}\mathbf{Z}^{H}]，

其中 (\mathbf{Z}) 是一个酉矩阵，(\mathbf{T}) 是上三角或准上三角，具体取决于是否请求实舒尔形式或复舒尔形式。当 (\mathbf{A}) 是实值时，实舒尔形式中 (\mathbf{T}) 和 (\mathbf{Z}) 都是实数值。当 (\mathbf{A}) 是实值矩阵时，实舒尔形式只是准上三角，因为 (2\times2) 块从主对角线伸出，对应于任何复值特征值。命令 linalg.schur 找到舒尔分解，而命令 linalg.rsf2csf 将 (\mathbf{T}) 和 (\mathbf{Z}) 从实舒尔形式转换为复舒尔形式。舒尔形式在计算矩阵函数时特别有用。

以下示例说明了舒尔分解：

>>> from scipy import linalg
>>> A = np.asmatrix('[1 3 2; 1 4 5; 2 3 6]')
>>> T, Z = linalg.schur(A)
>>> T1, Z1 = linalg.schur(A, 'complex')
>>> T2, Z2 = linalg.rsf2csf(T, Z)
>>> T
array([[ 9.90012467,  1.78947961, -0.65498528],
 [ 0\.        ,  0.54993766, -1.57754789],
 [ 0\.        ,  0.51260928,  0.54993766]])
>>> T2
array([[ 9.90012467+0.00000000e+00j, -0.32436598+1.55463542e+00j,
 -0.88619748+5.69027615e-01j],
 [ 0\.        +0.00000000e+00j,  0.54993766+8.99258408e-01j,
 1.06493862+3.05311332e-16j],
 [ 0\.        +0.00000000e+00j,  0\.        +0.00000000e+00j,
 0.54993766-8.99258408e-01j]])
>>> abs(T1 - T2) # different
array([[  1.06604538e-14,   2.06969555e+00,   1.69375747e+00],  # may vary
 [  0.00000000e+00,   1.33688556e-15,   4.74146496e-01],
 [  0.00000000e+00,   0.00000000e+00,   1.13220977e-15]])
>>> abs(Z1 - Z2) # different
array([[ 0.06833781,  0.88091091,  0.79568503],    # may vary
 [ 0.11857169,  0.44491892,  0.99594171],
 [ 0.12624999,  0.60264117,  0.77257633]])
>>> T, Z, T1, Z1, T2, Z2 = map(np.asmatrix,(T,Z,T1,Z1,T2,Z2))
>>> abs(A - Z*T*Z.H)  # same
matrix([[  5.55111512e-16,   1.77635684e-15,   2.22044605e-15],
 [  0.00000000e+00,   3.99680289e-15,   8.88178420e-16],
 [  1.11022302e-15,   4.44089210e-16,   3.55271368e-15]])
>>> abs(A - Z1*T1*Z1.H)  # same
matrix([[  4.26993904e-15,   6.21793362e-15,   8.00007092e-15],
 [  5.77945386e-15,   6.21798014e-15,   1.06653681e-14],
 [  7.16681444e-15,   8.90271058e-15,   1.77635764e-14]])
>>> abs(A - Z2*T2*Z2.H)  # same
matrix([[  6.02594127e-16,   1.77648931e-15,   2.22506907e-15],
 [  2.46275555e-16,   3.99684548e-15,   8.91642616e-16],
 [  8.88225111e-16,   8.88312432e-16,   4.44104848e-15]])

插值分解

scipy.linalg.interpolative 包含用于计算矩阵的插值分解（ID）的例程。对于秩为 (k \leq \min { m, n }) 的矩阵 (A \in \mathbb{C}^{m \times n})，这是一种分解

[A \Pi = \begin{bmatrix} A \Pi_{1} & A \Pi_{2} \end{bmatrix} = A \Pi_{1} \begin{bmatrix} I & T \end{bmatrix},]

其中 (\Pi = [\Pi_{1}, \Pi_{2}]) 是一个带有 (\Pi_{1} \in { 0, 1 }^{n \times k}) 的置换矩阵，即 (A \Pi_{2} = A \Pi_{1} T)。这等效地可以写成 (A = BP)，其中 (B = A \Pi_{1}) 和 (P = [I, T] \Pi^{\mathsf{T}}) 是 骨架 和 插值 矩阵。

另请参阅

scipy.linalg.interpolative — 获取更多信息。

矩阵函数

考虑具有泰勒级数展开的函数 (f\left(x\right))

[f\left(x\right)=\sum_{k=0}{\infty}\frac{f\left(0\right)}{k!}x^{k}.]

可以使用这个泰勒级数为方阵 (\mathbf{A}) 定义矩阵函数

[f\left(\mathbf{A}\right)=\sum_{k=0}{\infty}\frac{f\left(0\right)}{k!}\mathbf{A}^{k}.]

注意

尽管这提供了矩阵函数的一个有用表示，但这很少是计算矩阵函数的最佳方法。特别是，如果矩阵不可对角化，结果可能不准确。

指数和对数函数

矩阵指数是更常见的矩阵函数之一。实现矩阵指数的首选方法是使用缩放和 Padé 近似来计算 (e^{x})。此算法实现为 linalg.expm。

矩阵指数的逆是作为矩阵对数定义的矩阵对数：

[\mathbf{A}\equiv\exp\left(\log\left(\mathbf{A}\right)\right).]

可以使用 linalg.logm 获得矩阵对数。

三角函数

矩阵中的三角函数，(\sin), (\cos), 和 (\tan), 在 linalg.sinm, linalg.cosm, 和 linalg.tanm 中实现。矩阵正弦和余弦可以使用欧拉恒等式定义为

\begin{eqnarray} \sin\left(\mathbf{A}\right) & = & \frac{e{j\mathbf{A}}-e{-j\mathbf{A}}}{2j}\ \cos\left(\mathbf{A}\right) & = & \frac{e{j\mathbf{A}}+e{-j\mathbf{A}}}{2}.\end{eqnarray}

切线是

[\tan\left(x\right)=\frac{\sin\left(x\right)}{\cos\left(x\right)}=\left[\cos\left(x\right)\right]^{-1}\sin\left(x\right)]

因此，矩阵切线被定义为

[\left[\cos\left(\mathbf{A}\right)\right]^{-1}\sin\left(\mathbf{A}\right).]

超 bolic 三角函数

超 bolic 三角函数，(\sinh), (\cosh), 和 (\tanh), 也可以用熟悉的定义为矩阵定义：

\begin{eqnarray} \sinh\left(\mathbf{A}\right) & = & \frac{e{\mathbf{A}}-e{-\mathbf{A}}}{2}\ \cosh\left(\mathbf{A}\right) & = & \frac{e{\mathbf{A}}+e{-\mathbf{A}}}{2}\ \tanh\left(\mathbf{A}\right) & = & \left[\cosh\left(\mathbf{A}\right)\right]^{-1}\sinh\left(\mathbf{A}\right).\end{eqnarray}

这些矩阵函数可以在 linalg.sinhm, linalg.coshm, 和 linalg.tanhm 中找到。

任意函数

最后，任意接受一个复数并返回一个复数的函数可以被称为矩阵函数，使用命令 linalg.funm。此命令接受矩阵和任意的 Python 函数。然后，它使用 Golub 和 Van Loan 的书“Matrix Computations”中的算法，使用 Schur 分解计算应用于矩阵的函数。请注意函数需要接受复数作为输入才能使用此算法。例如，以下代码计算应用于矩阵的零阶贝塞尔函数。

>>> from scipy import special, linalg
>>> rng = np.random.default_rng()
>>> A = rng.random((3, 3))
>>> B = linalg.funm(A, lambda x: special.jv(0, x))
>>> A
array([[0.06369197, 0.90647174, 0.98024544],
 [0.68752227, 0.5604377 , 0.49142032],
 [0.86754578, 0.9746787 , 0.37932682]])
>>> B
array([[ 0.6929219 , -0.29728805, -0.15930896],
 [-0.16226043,  0.71967826, -0.22709386],
 [-0.19945564, -0.33379957,  0.70259022]])
>>> linalg.eigvals(A)
array([ 1.94835336+0.j, -0.72219681+0.j, -0.22270006+0.j])
>>> special.jv(0, linalg.eigvals(A))
array([0.25375345+0.j, 0.87379738+0.j, 0.98763955+0.j])
>>> linalg.eigvals(B)
array([0.25375345+0.j, 0.87379738+0.j, 0.98763955+0.j])

注意，由于矩阵解析函数的定义方式，贝塞尔函数已经作用于矩阵的特征值。

特殊矩阵

SciPy 和 NumPy 提供了几个用于创建在工程和科学中经常使用的特殊矩阵的函数。

类型 函数 描述
块对角 scipy.linalg.block_diag 从提供的数组创建一个块对角矩阵。
循环 scipy.linalg.circulant 创建一个循环矩阵。
伴随 scipy.linalg.companion 创建一个伴随矩阵。
卷积 scipy.linalg.convolution_matrix 创建一个卷积矩阵。
离散傅立叶 scipy.linalg.dft 创建一个离散傅立叶变换矩阵。
菲德勒 scipy.linalg.fiedler 创建一个对称的菲德勒矩阵。
菲德勒伴随 scipy.linalg.fiedler_companion 创建一个菲德勒伴随矩阵。
哈达玛德 scipy.linalg.hadamard 创建一个哈达玛德矩阵。
汉克尔 scipy.linalg.hankel 创建一个汉克尔矩阵。
赫尔默特 scipy.linalg.helmert 创建一个赫尔默特矩阵。
希尔伯特 scipy.linalg.hilbert 创建一个希尔伯特矩阵。
逆希尔伯特 scipy.linalg.invhilbert 创建希尔伯特矩阵的逆矩阵。
莱斯利 scipy.linalg.leslie 创建一个莱斯利矩阵。
帕斯卡 scipy.linalg.pascal 创建一个帕斯卡矩阵。
逆帕斯卡 scipy.linalg.invpascal 创建帕斯卡矩阵的逆矩阵。
托普利茨 scipy.linalg.toeplitz 创建一个托普利茨矩阵。
范德蒙德 numpy.vander 创建一个范德蒙德矩阵。

用这些函数的示例，请参阅它们各自的文档字符串。

稀疏数组（scipy.sparse）

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

简介

scipy.sparse 及其子模块提供了用于处理稀疏数组的工具。稀疏数组是只有数组中少数位置包含任何数据的数组，大多数位置被视为“空”。稀疏数组很有用，因为它们允许用于线性代数（scipy.sparse.linalg）或基于图的计算（scipy.sparse.csgraph）的算法更简单、更快速或内存消耗较少，但是它们通常在像切片、重塑或赋值等操作上不太灵活。本指南将介绍scipy.sparse中稀疏数组的基础知识，解释稀疏数据结构的独特之处，并引导用户查看用户指南中解释稀疏线性代数和图方法的其他部分。

入门稀疏数组

稀疏数组是一种特殊类型的数组，其中数组中只有少数位置包含数据。这允许使用压缩表示数据，仅记录存在数据的位置。有许多不同的稀疏数组格式，每种格式在压缩和功能之间进行不同的权衡。首先，让我们构建一个非常简单的稀疏数组，即坐标（COO）数组 (coo_array) 并将其与密集数组进行比较：

>>> import scipy as sp
>>> import numpy
>>> dense = numpy.array([[1, 0, 0, 2], [0, 4, 1, 0], [0, 0, 5, 0]])
>>> sparse = sp.sparse.coo_array(dense)
>>> dense
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])
>>> sparse
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 5 stored elements in COOrdinate format>

注意，在我们的密集数组中，有五个非零值。例如，2 在位置 0,3，4 在位置 1,1。所有其他值都为零。稀疏数组显式记录这五个值（参见 COOrdinate format 中的 5 个存储元素），然后将所有其余的零表示为隐式值。

大多数稀疏数组方法的工作方式与密集数组方法类似：

>>> sparse.max()
5
>>> dense.max()
5
>>> sparse.argmax()
10
>>> dense.argmax()
10
>>> sparse.mean()
1.0833333333333333
>>> dense.mean()
1.0833333333333333

稀疏数组还具有一些“额外”的属性，例如 .nnz，它返回存储值的数量：

>>> sparse.nnz
5

大多数减少操作，例如 .mean()、.sum() 或 .max()，在应用到稀疏数组的轴上时将返回一个 numpy 数组：

>>> sparse.mean(axis=1)
array([0.75, 1.25, 1.25])

这是因为稀疏数组上的减少操作通常是密集的。

理解稀疏数组格式

不同类型的稀疏数组具有不同的功能。例如，COO 数组不能被索引或切片：

>>> dense[2, 2]
5
>>> sparse[2, 2]
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
TypeError: 'coo_array' object is not subscriptable

但是，其他格式，例如压缩稀疏行（CSR）csr_array 支持切片和元素索引：

>>> sparse.tocsr()[2, 2]
5

有时，scipy.sparse会返回与输入稀疏矩阵格式不同的稀疏矩阵格式。例如，COO 格式的两个稀疏数组的点积将是 CSR 格式数组：

>>> sparse @ sparse.T
<3x3 sparse array of type '<class 'numpy.int64'>'
 with 5 stored elements in Compressed Sparse Row format>

这种改变是因为scipy.sparse会改变输入稀疏数组的格式，以使用最有效的计算方法。

scipy.sparse模块包含以下格式，每种格式都有自己独特的优势和劣势：

  • 块状稀疏行（BSR）数组scipy.sparse.bsr_array，在数组的数据部分以连续的块出现时最合适。

  • 坐标（COO）数组scipy.sparse.coo_array提供了一种简单的构建稀疏数组和原地修改它们的方法。COO 也可以快速转换为其他格式，如 CSR、CSC 或 BSR。

  • 压缩稀疏行（CSR）数组scipy.sparse.csr_array，最适用于快速算术运算、向量乘积和按行切片。

  • 压缩稀疏列（CSC）数组scipy.sparse.csc_array最适用于快速算术运算、向量乘积和按列切片。

  • 对角线（DIA）数组scipy.sparse.dia_array，对于有效存储和快速算术运算很有用，只要数据主要沿着数组的对角线出现。

  • 键值字典（DOK）数组scipy.sparse.dok_array，对于快速构建和单元素访问很有用。

  • 列表列表（LIL）数组scipy.sparse.lil_array，对于快速构建和修改稀疏数组很有用。

每种稀疏数组格式的优势和劣势的更多信息可以在它们的文档中找到。

所有scipy.sparse数组格式都可以直接从numpy.ndarray构建。然而，某些稀疏格式也可以以不同的方式构建。每个稀疏数组格式都有不同的优势，并且这些优势在每个类中都有文档记录。例如，构建稀疏数组最常见的方法之一是从单独的row、column和data值构建稀疏数组。对于我们之前的数组：

>>> dense
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])

row、column和data数组描述了稀疏数组中条目的行、列和值：

>>> row = [0,0,1,1,2]
>>> col = [0,3,1,2,2]
>>> data = [1,2,4,1,5]

使用这些，我们现在可以定义一个稀疏数组而不需要首先构建一个密集数组：

>>> csr = sp.sparse.csr_array((data, (row, col)))
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 5 stored elements in Compressed Sparse Row format>

不同的类有不同的构造函数，但是scipy.sparse.csr_array、scipy.sparse.csc_array和scipy.sparse.coo_array允许使用这种构造方式。

稀疏数组，隐式零和重复

稀疏数组很有用，因为它们表示大部分值是隐式的，而不是存储一个实际的占位值。在scipy.sparse中，用于表示“无数据”的值是隐式零。当需要显式零时，这可能会令人困惑。例如，在图方法中，我们经常需要能够区分（A）连接节点 i 和 j 的权重为零的链接和（B）i 和 j 之间没有链接。稀疏矩阵可以做到这一点，只要我们记住显式和隐式零即可。

例如，在我们之前的csr数组中，我们可以通过将它包含在data列表中来包含一个显式零。让我们把底行和最后一列数组中的最后一个条目视为显式零：

>>> row = [0,0,1,1,2,2]
>>> col = [0,3,1,2,2,3]
>>> data = [1,2,4,1,5,0]

那么，我们的稀疏数组将有六个存储元素，而不是五个：

>>> csr = sp.sparse.csr_array((data, (row, col)))
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 6 stored elements in Compressed Sparse Row format>

“额外”的元素就是我们的显式零。当转换回密集数组时，两者仍然是相同的，因为密集数组将所有东西都显式地表示：

>>> csr.todense()
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])
>>> dense
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])

但是，对于稀疏算术、线性代数和图方法，位置2,3处的值将被视为显式零。要去除此显式零，我们可以使用csr.eliminate_zeros()方法。这个方法在稀疏数组中原地操作，并移除任何零值存储元素：

>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 6 stored elements in Compressed Sparse Row format>
>>> csr.eliminate_zeros()
>>> csr
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 5 stored elements in Compressed Sparse Row format>

在csr.eliminate_zeros()之前，有六个存储元素。之后，只有五个存储元素。

另一个复杂性点源于在构建稀疏数组时处理 重复项 的方式。当我们在构建稀疏数组时在 row,col 处有两个或更多条目时，就会发生 重复项。例如，我们可能用重复值在 1,1 处表示先前的数组：

>>> row = [0,0,1,1,1,2]
>>> col = [0,3,1,1,2,2]
>>> data = [1,2,1,3,1,5]

在这种情况下，我们可以看到有 两个 data 值对应于我们最终数组中的 1,1 位置。scipy.sparse 将单独存储这些值：

>>> dupes = sp.sparse.coo_array((data, (row, col)))
>>> dupes
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 6 stored elements in COOrdinate format>

请注意，这个稀疏数组中有六个存储的元素，尽管只有五个唯一的数据位置。当这些数组转换回密集数组时，重复值将被求和。因此，在位置 1,1 处，密集数组将包含重复存储条目的总和，即 1 + 3：

>>> dupes.todense()
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])

要删除稀疏数组本身中的重复值，从而减少存储元素的数量，可以使用 .sum_duplicates() 方法：

>>> dupes.sum_duplicates()
>>> dupes
<3x4 sparse array of type '<class 'numpy.int64'>'
 with 5 stored elements in COOrdinate format>

现在我们的稀疏数组中只有五个存储的元素，且与本指南中一直使用的数组相同：

>>> dupes.todense()
array([[1, 0, 0, 2],
 [0, 4, 1, 0],
 [0, 0, 5, 0]])

规范格式

几种稀疏数组格式具有“规范格式”，以实现更高效的操作。通常这些格式包括像增加限制这样的额外条件：

  • 任何值都没有重复条目

  • 已排序的索引

具有规范形式的类包括：coo_array，csr_array，csc_array 和 bsr_array。详细信息请参阅这些类的文档字符串，了解每种规范表示的细节。

要检查这些类的实例是否处于规范形式，请使用 .has_canonical_format 属性：

>>> coo = sp.sparse.coo_array(([1, 1, 1], ([0, 2, 1], [0, 1, 2])))
>>> coo.has_canonical_format
False

要将实例转换为规范形式，请使用 .sum_duplicates() 方法：

>>> coo.sum_duplicates()
>>> coo.has_canonical_format
True

稀疏数组的下一步操作

当处理大型、几乎为空的数组时，稀疏数组类型最为有用。特别是在这些情况下，稀疏线性代数和稀疏图方法的效率显著提高。

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

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

介绍

ARPACK [1] 是一个 Fortran 包，提供快速查找大稀疏矩阵少数特征值/特征向量的例程。为了找到这些解，它只需左乘问题矩阵。此操作通过反向通信接口执行。这种结构的结果是 ARPACK 能够找到任何线性函数映射向量到向量的特征值和特征向量。

ARPACK 提供的所有功能都包含在两个高级接口中 scipy.sparse.linalg.eigs 和 scipy.sparse.linalg.eigsh。 eigs 提供实现接口，用于查找实数或复数非对称方阵的特征值/特征向量，而 eigsh 提供接口，用于实对称或复共轭矩阵。

基本功能

ARPACK 可以解决形如

[A \mathbf{x} = \lambda \mathbf{x}]

或一般的特征值问题形式

[A \mathbf{x} = \lambda M \mathbf{x}.]

ARPACK 的优势在于它可以计算特定子集的特征值/特征向量对。这是通过关键字which实现的。以下which值可用：

  • which = 'LM' : 最大幅度的特征值 (eigs, eigsh), 即复数的欧几里德范数中的最大特征值.

  • which = 'SM' : 最小幅度的特征值 (eigs, eigsh), 即复数的欧几里德范数中的最小特征值.

  • which = 'LR' : 最大实部的特征值 (eigs).

  • which = 'SR' : 最小实部的特征值 (eigs).

  • which = 'LI' : 最大虚部的特征值 (eigs).

  • which = 'SI' : 最小虚部的特征值 (eigs).

  • which = 'LA' : 最大代数值的特征值 (eigsh), 即包含任何负号的最大特征值.

  • which = 'SA' : 最小代数值的特征值 (eigsh), 即包含任何负号的最小特征值.

  • which = 'BE' : 光谱两端的特征值 (eigsh).

注意，ARPACK 通常更擅长找到极端特征值，即具有较大幅度的特征值。特别是，使用which = 'SM'可能导致执行时间缓慢和/或异常结果。更好的方法是使用转移反演模式。

转移反演模式

Shift-invert mode 依赖于以下观察。对于广义特征值问题

[A \mathbf{x} = \lambda M \mathbf{x},]

可以证明

[(A - \sigma M)^{-1} M \mathbf{x} = \nu \mathbf{x},]

其中

[\nu = \frac{1}{\lambda - \sigma}.]

举例

假设您想要查找大矩阵的最小和最大特征值以及相应的特征向量。ARPACK 可以处理多种输入形式：如密集矩阵，例如numpy.ndarray 实例，稀疏矩阵，例如scipy.sparse.csr_matrix，或者从scipy.sparse.linalg.LinearOperator 派生的一般线性操作员。为了简单起见，在本例中，我们将构建一个对称的正定矩阵。

>>> import numpy as np
>>> from scipy.linalg import eig, eigh
>>> from scipy.sparse.linalg import eigs, eigsh
>>> np.set_printoptions(suppress=True)
>>> rng = np.random.default_rng()
>>>
>>> X = rng.random((100, 100)) - 0.5
>>> X = np.dot(X, X.T)  # create a symmetric matrix

现在我们有一个对称矩阵 X，用来测试这些程序。首先，使用 eigh 计算标准特征值分解：

>>> evals_all, evecs_all = eigh(X)

随着 X 的维度增长，这个例程变得非常慢。特别是，如果只需要少量特征向量和特征值，ARPACK 可能是一个更好的选择。首先，计算 X 的最大特征值 (which = 'LM') 并将其与已知结果进行比较：

>>> evals_large, evecs_large = eigsh(X, 3, which='LM')
>>> print(evals_all[-3:])
[29.22435321 30.05590784 30.58591252]
>>> print(evals_large)
[29.22435321 30.05590784 30.58591252]
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1\.  0\.  0.],       # may vary (signs)
 [ 0\.  1\.  0.],
 [-0\.  0\. -1.]])

结果如预期。ARPACK 恢复了所需的特征值，并且它们与先前已知的结果相匹配。此外，特征向量是正交的，这是我们预期的。现在，让我们尝试解最小幅度特征值的问题：

>>> evals_small, evecs_small = eigsh(X, 3, which='SM')
Traceback (most recent call last):       # may vary (convergence)
...
scipy.sparse.linalg._eigen.arpack.arpack.ArpackNoConvergence:
ARPACK error -1: No convergence (1001 iterations, 0/3 eigenvectors converged)

糟糕。我们看到，如上所述，ARPACK 并不太擅长找到小特征值。可以通过几种方法来解决这个问题。我们可以增加容差 (tol) 以加快收敛速度：

>>> evals_small, evecs_small = eigsh(X, 3, which='SM', tol=1E-2)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 0.99999999  0.00000024 -0.00000049],    # may vary (signs)
 [-0.00000023  0.99999999  0.00000056],
 [ 0.00000031 -0.00000037  0.99999852]])

这种方法有效，但结果的精度会降低。另一种选择是将最大迭代次数 (maxiter) 从 1000 增加到 5000：

>>> evals_small, evecs_small = eigsh(X, 3, which='SM', maxiter=5000)
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1\.  0\.  0.],           # may vary (signs)
 [-0\.  1\.  0.],
 [ 0\.  0\. -1.]])

我们得到了预期的结果，但计算时间要长得多。幸运的是，ARPACK 包含一个模式，可以快速确定非外部特征值：shift-invert mode。如上所述，这种模式涉及将特征值问题转换为具有不同特征值的等价问题。在这种情况下，我们希望找到接近零的特征值，因此我们选择 sigma = 0。然后转换后的特征值将满足 (\nu = 1/(\lambda - \sigma) = 1/\lambda)，因此我们的小特征值 (\lambda) 变为大特征值 (\nu)。

>>> evals_small, evecs_small = eigsh(X, 3, sigma=0, which='LM')
>>> evals_all[:3]
array([0.00053181, 0.00298319, 0.01387821])
>>> evals_small
array([0.00053181, 0.00298319, 0.01387821])
>>> np.dot(evecs_small.T, evecs_all[:,:3])
array([[ 1\.  0\.  0.],    # may vary (signs)
 [ 0\. -1\. -0.],
 [-0\. -0\.  1.]])

我们得到了我们希望的结果，计算时间大大减少。请注意，从 (\nu \to \lambda) 的转换完全在后台进行。用户不必担心细节。

移位-反转模式提供的不仅仅是获取少量小特征值的快速方法。比如，您希望找到内部特征值和特征向量，例如那些接近(\lambda = 1)的。只需设置sigma = 1，ARPACK 将处理其余部分：

>>> evals_mid, evecs_mid = eigsh(X, 3, sigma=1, which='LM')
>>> i_sort = np.argsort(abs(1. / (1 - evals_all)))[-3:]
>>> evals_all[i_sort]
array([0.94164107, 1.05464515, 0.99090277])
>>> evals_mid
array([0.94164107, 0.99090277, 1.05464515])
>>> print(np.dot(evecs_mid.T, evecs_all[:,i_sort]))
array([[-0\.  1\.  0.],     # may vary (signs)
 [-0\. -0\.  1.],
 [ 1\.  0\.  0.]]

特征值的顺序不同，但它们都在那里。请注意，移位-反转模式需要内部解决矩阵的逆问题。这由eigsh和eigs自动处理，但用户也可以指定该操作。有关详细信息，请参阅scipy.sparse.linalg.eigsh和scipy.sparse.linalg.eigs的文档字符串。

使用 LinearOperator

现在我们考虑一个情况，您想避免创建密集矩阵，而是使用scipy.sparse.linalg.LinearOperator。我们的第一个线性算子应用于输入向量和用户提供给算子本身的向量(\mathbf{d})之间的逐元素乘法。这个算子模拟了一个对角矩阵，其主对角线上的元素是(\mathbf{d})，它的主要优点在于前向和伴随操作都是简单的逐元素乘法，而不是矩阵-向量乘法。对于对角矩阵，我们期望的特征值等于沿主对角线的元素，即(\mathbf{d})。使用eigsh得到的特征值和特征向量与应用于密集矩阵时使用eigh得到的进行比较：

>>> from scipy.sparse.linalg import LinearOperator
>>> class Diagonal(LinearOperator):
...     def __init__(self, diag, dtype='float32'):
...         self.diag = diag
...         self.shape = (len(self.diag), len(self.diag))
...         self.dtype = np.dtype(dtype)
...     def _matvec(self, x):
...         return self.diag*x
...     def _rmatvec(self, x):
...         return self.diag*x 

>>> N = 100
>>> rng = np.random.default_rng()
>>> d = rng.normal(0, 1, N).astype(np.float64)
>>> D = np.diag(d)
>>> Dop = Diagonal(d, dtype=np.float64) 

>>> evals_all, evecs_all = eigh(D)
>>> evals_large, evecs_large = eigsh(Dop, 3, which='LA', maxiter=1e3)
>>> evals_all[-3:]
array([1.53092498, 1.77243671, 2.00582508])
>>> evals_large
array([1.53092498, 1.77243671, 2.00582508])
>>> print(np.dot(evecs_large.T, evecs_all[:,-3:]))
array([[-1\.  0\.  0.],     # may vary (signs)
 [-0\. -1\.  0.],
 [ 0\.  0\. -1.]]

在这种情况下，我们创建了一个快速且简单的Diagonal算子。外部库PyLops提供了与Diagonal算子类似的功能，以及其他几个算子。

最后，我们考虑一个线性算子，模仿一阶导数模板的应用。在这种情况下，该算子等效于一个实非对称矩阵。再次，我们将估计的特征值和特征向量与将相同的一阶导数应用于输入信号的密集矩阵进行比较：

>>> class FirstDerivative(LinearOperator):
...     def __init__(self, N, dtype='float32'):
...         self.N = N
...         self.shape = (self.N, self.N)
...         self.dtype = np.dtype(dtype)
...     def _matvec(self, x):
...         y = np.zeros(self.N, self.dtype)
...         y[1:-1] = (0.5*x[2:]-0.5*x[0:-2])
...         return y
...     def _rmatvec(self, x):
...         y = np.zeros(self.N, self.dtype)
...         y[0:-2] = y[0:-2] - (0.5*x[1:-1])
...         y[2:] = y[2:] + (0.5*x[1:-1])
...         return y 

>>> N = 21
>>> D = np.diag(0.5*np.ones(N-1), k=1) - np.diag(0.5*np.ones(N-1), k=-1)
>>> D[0] = D[-1] = 0 # take away edge effects
>>> Dop = FirstDerivative(N, dtype=np.float64) 

>>> evals_all, evecs_all = eig(D)
>>> evals_large, evecs_large = eigs(Dop, 4, which='LI')
>>> evals_all_imag = evals_all.imag
>>> isort_imag = np.argsort(np.abs(evals_all_imag))
>>> evals_all_imag = evals_all_imag[isort_imag]
>>> evals_large_imag = evals_large.imag
>>> isort_imag = np.argsort(np.abs(evals_large_imag))
>>> evals_large_imag = evals_large_imag[isort_imag]
>>> evals_all_imag[-4:]
array([-0.95105652, 0.95105652, -0.98768834, 0.98768834])
>>> evals_large_imag
array([0.95105652, -0.95105652, 0.98768834, -0.98768834]) # may vary

注意，这个算子的特征值都是虚数。此外，scipy.sparse.linalg.eigs函数的关键字which='LI'会产生具有最大绝对虚部（正负都有）的特征值。同样，在PyLops库中有一个更高级的一阶导数算子的实现，名为FirstDerivative算子。

参考资料

posted @ 2024-06-27 17:12  绝不原创的飞龙  阅读(124)  评论(0)    收藏  举报