线性代数 (scipy.linalg)#

当 SciPy 使用优化的 ATLAS LAPACK 和 BLAS 库构建时,它具有非常快速的线性代数功能。如果深入挖掘,所有原始的 LAPACK 和 BLAS 库都可以供您使用,以获得更快的速度。在本节中,将介绍一些更易于使用这些例程的接口。

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

scipy.linalg 与 numpy.linalg#

scipy.linalg 包含了 numpy.linalg 中的所有函数,外加一些 numpy.linalg 中没有的更高级的函数。

使用 scipy.linalg 而非 numpy.linalg 的另一个优势是它总是编译有 BLAS/LAPACK 支持,而对于 NumPy 来说这是可选的。因此,SciPy 版本可能会更快,具体取决于 NumPy 的安装方式。

因此,除非您不想将 scipy 添加为您 numpy 程序的依赖项,否则请使用 scipy.linalg 代替 numpy.linalg

numpy.matrix 与二维 numpy.ndarray#

表示矩阵的类以及基本的矩阵乘法和转置等操作是 numpy 的一部分。为方便起见,我们在此总结 numpy.matrixnumpy.ndarray 之间的区别。

numpy.matrix 是一个矩阵类,对于矩阵操作,它比 numpy.ndarray 具有更方便的接口。例如,该类支持通过分号实现类似 MATLAB 的创建语法,默认将 * 运算符作为矩阵乘法,并包含作为逆和转置快捷方式的 IT 成员。

>>> 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 中,NumPy 数组 A 的矩阵逆通过 linalg.inv (A) 获取,如果 A 是 Matrix,则使用 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\) 行和第 \(j\) 列后留下的矩阵的行列式。那么,对于任何行 \(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 参数使用不同的参数,可以使用多种范数定义。此函数接受秩为 1(向量)或秩为 2(矩阵)的数组以及一个可选的 order 参数(默认为 2)。基于这些输入,计算出所需阶数的向量或矩阵范数。

对于向量 x,order 参数可以是任何实数,包括 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.0
>>> linalg.norm(A, -1)
4.0
>>> linalg.norm(A, np.inf) # L inf norm (max row sum)
7.0

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

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

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

其中 \(\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_{j}\left(x_{i}\right).\]

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

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

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

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

命令 linalg.lstsq 将在给定 \(\mathbf{A}\)\(\mathbf{y}\) 的情况下求解 \(\mathbf{c}\) 的线性最小二乘问题。此外,linalg.pinv 将在给定 \(\mathbf{A}\) 的情况下求出 \(\mathbf{A}^{\dagger}.\)

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

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

其中 \(x_{i}=0.1i\)\(i=1, \ldots, 10\)\(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}^{H}\mathbf{A}\right)^{-1}\mathbf{A}^{H},\]

而如果 \(M<N\) 矩阵,广义逆为

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

\(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}^{H}\]

或者:

\[\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}^{H}\mathbf{v}_{L}\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}^{2}=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^{2}+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\) 矩阵。矩阵 \(\mathbf{A}^{H}\mathbf{A}\)\(\mathbf{A}\mathbf{A}^{H}\) 分别是大小为 \(N\times N\)\(M\times M\) 的正方埃尔米特矩阵 [1]。众所周知,埃尔米特矩阵的特征值是实数且非负的。此外,\(\mathbf{A}^{H}\mathbf{A}\)\(\mathbf{A}\mathbf{A}^{H}\) 最多有 \(\min\left(M,N\right)\) 个相同的非零特征值。将这些正特征值定义为 \(\sigma_{i}^{2}.\) 它们的平方根称为 \(\mathbf{A}\) 的奇异值。\(\mathbf{A}^{H}\mathbf{A}\) 的特征向量按列收集到 \(N\times N\) 酉矩阵 [2] \(\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}_{i}\]

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

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

由于 \(\mathbf{L}\) 是下三角矩阵,该方程可以针对 \(\mathbf{U}\mathbf{x}_{i}\) 求解,并最终使用前向和后向代入非常快速地求解 \(\mathbf{x}_{i}\)。花费最初的时间对 \(\mathbf{A}\) 进行因子分解,可以在未来非常快速地求解类似的方程组。如果执行 LU 分解的目的是为了求解线性系统,则应使用命令 linalg.lu_factor,随后重复应用命令 linalg.lu_solve 为每个新的右侧求解系统。

Cholesky 分解#

Cholesky 分解是 LU 分解的一种特殊情况,适用于埃尔米特正定矩阵。当 \(\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}.\) 命令 linalg.cholesky 计算 Cholesky 分解。对于使用 Cholesky 分解求解方程组,还有 linalg.cho_factorlinalg.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}\) 是上三角或准上三角矩阵,取决于是否请求实 Schur 形式或复 Schur 形式。对于实 Schur 形式,当 \(\mathbf{A}\) 是实值时,\(\mathbf{T}\)\(\mathbf{Z}\) 都是实值的。当 \(\mathbf{A}\) 是实值矩阵时,实 Schur 形式仅为准上三角,因为对应于任何复数值特征值的 \(2\times2\) 块会从主对角线挤出。命令 linalg.schur 查找 Schur 分解,而命令 linalg.rsf2csf\(\mathbf{T}\)\(\mathbf{Z}\) 从实 Schur 形式转换为复 Schur 形式。Schur 形式在计算矩阵函数时特别有用。

以下示例说明了 Schur 分解

>>> 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(k\right)}\left(0\right)}{k!}x^{k}.\]

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

\[f\left(\mathbf{A}\right)=\sum_{k=0}^{\infty}\frac{f^{\left(k\right)}\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.sinmlinalg.cosmlinalg.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).\]

双曲三角函数#

双曲三角函数 \(\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.sinhmlinalg.coshmlinalg.tanhm 找到。

任意函数#

最后,任何接受一个复数并返回一个复数的任意函数都可以使用命令 linalg.funm 作为矩阵函数调用。此命令接受矩阵和一个任意的 Python 函数。然后它实现 Golub 和 Van Loan 的《矩阵计算》中的一种算法,以使用 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

创建离散傅里叶变换矩阵。

Fiedler

scipy.linalg.fiedler

创建对称 Fiedler 矩阵。

Fiedler Companion

scipy.linalg.fiedler_companion

创建 Fiedler 伴随矩阵。

Hadamard

scipy.linalg.hadamard

创建 Hadamard 矩阵。

Hankel

scipy.linalg.hankel

创建 Hankel 矩阵。

Helmert

scipy.linalg.helmert

创建 Helmert 矩阵。

Hilbert

scipy.linalg.hilbert

创建 Hilbert 矩阵。

Inverse Hilbert

scipy.linalg.invhilbert

创建 Hilbert 矩阵的逆矩阵。

Leslie

scipy.linalg.leslie

创建 Leslie 矩阵。

Pascal

scipy.linalg.pascal

创建 Pascal 矩阵。

Inverse Pascal

scipy.linalg.invpascal

创建 Pascal 矩阵的逆矩阵。

Toeplitz

scipy.linalg.toeplitz

创建 Toeplitz 矩阵。

Van der Monde

numpy.vander

创建 Van der Monde 矩阵。

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

高级特性#

批处理支持#

SciPy 的一些线性代数函数可以在给定 N 维数组输入时处理标量、一维或二维数组的批次。例如,通常接受(二维)矩阵的线性代数函数可能接受形状为 (4, 3, 2) 的数组,它将其解释为四个 3x2 矩阵的批次。在这种情况下,我们称输入的*核心形状*为 (3, 2),而*批次形状*为 (4,)。同样,通常接受(一维)向量的线性代数函数会将 (4, 3, 2) 数组视为 (4, 3) 的向量批次,在这种情况下,输入的*核心形状*为 (2,),而*批次形状*为 (4, 3)。核心形状的长度也称为*核心维度*。在这些情况下,输出的最终形状是输入的批次形状与输出的核心形状连接(即当输入的批次形状为 () 时输出的形状)。更多信息,请参阅 批处理线性运算

overwrite_* 参数#

许多线性代数函数具有 overwrite_* 参数,用于表明操作允许就地进行并覆盖输入数组(命名约定是数组参数 a 有一个对应的 overwrite_ab 参数有一个对应的 overwrite_b,依此类推)。

默认情况下,scipy.linalg 函数会保留它们的输入并在内部进行复制。在某些情况下,就地工作可能会提高性能或避免耗尽内存——但建议进行测量以确保您确实在特定工作负载中看到了改进。

但请注意,设置 overwrite_a=True 仅表示它*可能*就地工作;它是否真的这样做取决于是否满足附加要求。对于大多数函数,就地操作要求以下所有条件为真

  • 输入数组的数据类型是 LAPACK 兼容的:只有单精度和双精度浮点数(float{32,64})和复数(complex{64,128})数组是兼容的;

  • 数组是 Fortran 顺序的;

  • 目前只有二维数组(ndim == 2)是兼容的;批处理数组使用核心形状数组大小的内部缓冲区。但这在未来的 SciPy 版本中可能会改变。

如果违反了任何这些条件,则会进行内部复制,无论 overwrite_aTrue 还是 False

某些函数可能会放宽其中一些要求或施加额外的就地工作约束。

请注意,这一切都是完全不透明的:如果您设置 overwrite_a=True 但函数必须在副本上操作,则不会生成警告或诊断输出。

作为演示,考虑一个简单的求逆操作

>>> import numpy as np
>>> from scipy.linalg import inv
>>> a = np.asarray([[2, 1], [0, 1]])
>>> a_inv = inv(a)
>>> a_inv
array([[ 0.5, -0.5],
       [ 0. ,  1. ]])
>>> a
array([[2, 1],
       [0, 1]])

由于 a 是整数数据类型,scipy.linalg.inv 函数在内部制作了一个浮点副本,并将其提供给 LAPACK 的求逆例程。因此将 overwrite_a 设置为 True 没有区别

>>> a_inv = inv(a, overwrite_a=True)
>>> a_inv
array([[ 0.5, -0.5],
       [ 0. ,  1. ]])
>>> a
array([[2, 1],
       [0, 1]])

要强制它就地工作,我们需要手动将 a 数组转换为 Fortran 顺序且数据类型兼容

>>> af = np.asarray(a, order='F', dtype=np.float64)
>>> a_inv = inv(af, overwrite_a=True)
>>> a_inv
array([[ 0.5, -0.5],
       [ 0. ,  1. ]])
>>> af
array([[ 0.5, -0.5],
       [ 0. ,  1. ]])

注意输入数组 af 现在已被覆盖。我们可以检查输入和输出数组是否确实使用了相同的内存缓冲区

>>> np.shares_memory(af, a_inv)
True

通常,a_inv is af 可能成立也可能不成立:唯一的保证是输入和输出共享内存缓冲区的*某个*部分。

最后,我们建议您仅在性能/内存测量显示有显著改进时才使用这些 overwrite_a 参数。在绝大多数使用场景中,添加它们属于过早优化。