线性代数 (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 与 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 中，NumPy 数组 A 的矩阵逆通过使用 linalg.inv (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 参数，可以使用多种范数定义。此函数接受一个秩 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.lstsq 和 linalg.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 矩阵，则广义逆矩阵为

\[\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 分解是 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}.\]

注意

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

指数和对数函数

矩阵指数是更常见的矩阵函数之一。实现矩阵指数的首选方法是使用缩放和 \(e^{x}\) 的 Padé 近似。此算法实现为 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 的 书 “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	创建离散傅里叶变换矩阵。
Fiedler	scipy.linalg.fiedler	创建对称 Fiedler 矩阵。
Fiedler 伴随	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 矩阵。
逆 Hilbert	scipy.linalg.invhilbert	创建 Hilbert 矩阵的逆。
Leslie	scipy.linalg.leslie	创建 Leslie 矩阵。
Pascal	scipy.linalg.pascal	创建 Pascal 矩阵。
逆 Pascal	scipy.linalg.invpascal	创建 Pascal 矩阵的逆。
Toeplitz	scipy.linalg.toeplitz	创建 Toeplitz 矩阵。
Van der Monde	numpy.vander	创建 Van der Monde 矩阵。

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

高级功能

批量支持

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