线性代数 (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}\)，norm 的唯一有效值为 \(\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}\) 与数据 \(\mathbf{x}_{i}\) 通过一组系数 \(c_{j}\) 和模型函数 \(f_{j}\left(\mathbf{x}_{i}\right)\) 通过模型相关联

\[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{c}\) 求解线性最小二乘问题，给定 \(\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},\]

其中 \(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)\)）且对角线为单位，而 \(\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 分解的一个特例，适用于 Hermitian 正定矩阵。当 \(\mathbf{A}=\mathbf{A}^{H}\) 且 \(\mathbf{x}^{H}\mathbf{Ax}\geq0\) 对所有 \(\mathbf{x}\) 成立时，可以找到 \(\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 的书“矩阵计算”中的一个算法，使用 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	创建范德蒙德矩阵。

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