线性代数 (scipy.linalg)

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

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

scipy.linalg vs 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 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 类，因为它没有添加任何 2D 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 或 2D numpy.ndarray 对象。

基本例程

求逆矩阵

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

\[\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 参数的不同参数，可以使用各种范数定义。此函数采用秩 1（向量）或秩 2（矩阵）数组和一个可选的 order 参数（默认为 2）。基于这些输入，计算请求的阶数的向量或矩阵范数。

对于向量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.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.lstsq 和 linalg.pinv 解决数据拟合问题。下面显示的数据是使用模型生成的：

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

其中 \(x_{i}=0.1i\)，对于 \(i=1\ldots10\)，\(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\) 矩阵，其中 \(M\) 和 \(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.]])

[1]

一个厄米矩阵 \(\mathbf{D}\) 满足 \(\mathbf{D}^{H}=\mathbf{D}.\)

[2]

一个酉矩阵 \(\mathbf{D}\) 满足 \(\mathbf{D}^{H}\mathbf{D}=\mathbf{I}=\mathbf{D}\mathbf{D}^{H}\)，因此 \(\mathbf{D}^{-1}=\mathbf{D}^{H}.\)

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)\)），其对角线元素为单位 1，而 \(\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 分解是适用于 Hermitian 正定矩阵的 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_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}\) 是上三角矩阵或准上三角矩阵，具体取决于是否请求实 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.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).\]

双曲三角函数

双曲三角函数 \(\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 的著作《矩阵计算》中的算法，使用 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	创建范德蒙矩阵。

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