线性代数 (`scipy.linalg`)#

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

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

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 类很方便，但我们不鼓励使用它，因为它没有提供 2-D 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}\) 是主对角线上全为一的单位矩阵。通常，\(\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\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.]])

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

线性代数 (scipy.linalg)#