线性代数 (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.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 是 Matrix,则使用 A.I。例如,令
那么
下面的示例演示了在 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 可以直接求解线性方程组。此命令期望输入一个矩阵和一个右侧向量,然后计算出解向量。此外还提供了一个输入对称矩阵的选项,在适用的情况下可以加快处理速度。作为一个示例,假设需要求解以下联立方程
我们可以使用矩阵逆找到解向量
然而,最好使用 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,\)
这是一种定义行列式的递归方法,其中基本情况定义为接受 \(1\times1\) 矩阵的行列式即为唯一的矩阵元素。在 SciPy 中,行列式可以用 linalg.det 计算。例如,以下矩阵的行列式
为
在 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。计算出的范数为
对于矩阵 \(\mathbf{A}\),范数的唯一有效值为 \(\pm2,\pm1,\) \(\pm\) inf 和 ‘fro’(或 ‘f’)。因此,
其中 \(\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}\) 相关联,模型如下
其中 \(\epsilon_{i}\) 表示数据中的不确定性。最小二乘法的策略是选择系数 \(c_{j}\) 以最小化
理论上,全局最小值出现在
或者:
其中
当 \(\mathbf{A^{H}A}\) 可逆时,
其中 \(\mathbf{A}^{\dagger}\) 称为 \(\mathbf{A}\) 的伪逆。注意,使用这个 \(\mathbf{A}\) 的定义,模型可以写成
命令 linalg.lstsq 将在给定 \(\mathbf{A}\) 和 \(\mathbf{y}\) 的情况下求解 \(\mathbf{c}\) 的线性最小二乘问题。此外,linalg.pinv 将在给定 \(\mathbf{A}\) 的情况下求出 \(\mathbf{A}^{\dagger}.\)
以下示例和图形演示了使用 linalg.lstsq 和 linalg.pinv 求解数据拟合问题。下面显示的数据是使用以下模型生成的
其中 \(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\),广义逆为
而如果 \(M<N\) 矩阵,广义逆为
在 \(M=N\) 的情况下,
只要 \(\mathbf{A}\) 是可逆的。
分解#
在许多应用中,使用其他表示法分解矩阵很有用。SciPy 支持几种分解。
特征值和特征向量#
特征值-特征向量问题是常用的线性代数运算之一。在一种流行形式中,特征值-特征向量问题是针对某个方阵 \(\mathbf{A}\) 寻找标量 \(\lambda\) 和相应的向量 \(\mathbf{v}\),使得
对于 \(N\times N\) 矩阵,有 \(N\) 个(不一定不同)特征值——(特征)多项式的根
特征向量 \(\mathbf{v}\) 有时也称为右特征向量,以区别于满足以下条件的另一组左特征向量
或者:
使用其默认可选参数,命令 linalg.eig 返回 \(\lambda\) 和 \(\mathbf{v}.\) 但是,它也可以只返回 \(\mathbf{v}_{L}\) 和 \(\lambda\)(linalg.eigvals 也只返回 \(\lambda\))。
此外,linalg.eig 还可以求解更一般的特征值问题
用于方阵 \(\mathbf{A}\) 和 \(\mathbf{B}.\) 标准特征值问题是 \(\mathbf{B}=\mathbf{I}\) 的一般特征值问题的一个例子。当可以求解广义特征值问题时,它提供了 \(\mathbf{A}\) 的分解如下
其中 \(\mathbf{V}\) 是以特征向量为列的集合,\(\boldsymbol{\Lambda}\) 是特征值的对角矩阵。
根据定义,特征向量仅定义到常数比例因子。在 SciPy 中,特征向量的缩放因子被选择为 \(\left\Vert \mathbf{v}\right\Vert ^{2}=\sum_{i}v_{i}^{2}=1.\)
作为一个示例,考虑求以下矩阵的特征值和特征向量
特征多项式为
这个多项式的根是 \(\mathbf{A}\) 的特征值
对应于每个特征值的特征向量可以使用原始方程找到。然后可以找到与这些特征值相关联的特征向量。
>>> 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}\) 的奇异值分解。每个矩阵都有奇异值分解。有时,奇异值被称为 \(\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{P}\) 是 \(M\times M\) 置换矩阵(单位矩阵行的置换),\(\mathbf{L}\) 是 \(M\times K\) 下三角或梯形矩阵(\(K=\min\left(M,N\right)\)),对角线为单位元,\(\mathbf{U}\) 是上三角或梯形矩阵。此分解的 SciPy 命令是 linalg.lu。
这种分解通常对于求解许多联立方程很有用,其中左侧不变但右侧在变。例如,假设我们要求解
针对许多不同的 \(\mathbf{b}_{i}\)。LU 分解允许将其写为
由于 \(\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}\) 的分解,使得
其中 \(\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}\) 的 SVD,则可以找到 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{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}\),这是一个因子分解
其中 \(\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)\)
可以使用此泰勒级数将矩阵函数定义为方阵 \(\mathbf{A}\)
注意
虽然这可作为矩阵函数的一种有用表示,但它很少是计算矩阵函数的最佳方式。特别地,如果矩阵不可对角化,结果可能不准确。
指数和对数函数#
矩阵指数是更常见的矩阵函数之一。实现矩阵指数的首选方法是使用缩放和 Padé 近似来计算 \(e^{x}\)。此算法实现为 linalg.expm。
矩阵指数的逆是矩阵对数,定义为矩阵指数的逆
矩阵对数可以通过 linalg.logm 获得。
三角函数#
三角函数 \(\sin\)、\(\cos\) 和 \(\tan\) 已针对 linalg.sinm、linalg.cosm 和 linalg.tanm 中的矩阵实现。矩阵正弦和余弦可以使用欧拉恒等式定义为
正切是
因此矩阵正切定义为
双曲三角函数#
双曲三角函数 \(\sinh\)、\(\cosh\) 和 \(\tanh\) 也可以使用熟悉的定义为矩阵定义
这些矩阵函数可以使用 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 提供了几种用于创建工程和科学中常用的特殊矩阵的函数。
类型 |
函数 |
描述 |
|---|---|---|
分块对角 |
从提供的数组创建分块对角矩阵。 |
|
循环 |
创建循环矩阵。 |
|
伴随 |
创建伴随矩阵。 |
|
卷积 |
创建卷积矩阵。 |
|
离散傅里叶 |
创建离散傅里叶变换矩阵。 |
|
Fiedler |
创建对称 Fiedler 矩阵。 |
|
Fiedler Companion |
创建 Fiedler 伴随矩阵。 |
|
Hadamard |
创建 Hadamard 矩阵。 |
|
Hankel |
创建 Hankel 矩阵。 |
|
Helmert |
创建 Helmert 矩阵。 |
|
Hilbert |
创建 Hilbert 矩阵。 |
|
Inverse Hilbert |
创建 Hilbert 矩阵的逆矩阵。 |
|
Leslie |
创建 Leslie 矩阵。 |
|
Pascal |
创建 Pascal 矩阵。 |
|
Inverse Pascal |
创建 Pascal 矩阵的逆矩阵。 |
|
Toeplitz |
创建 Toeplitz 矩阵。 |
|
Van der Monde |
创建 Van der Monde 矩阵。 |
有关这些函数的使用示例,请参阅它们各自的文档字符串。
高级特性#
批处理支持#
SciPy 的一些线性代数函数可以在给定 N 维数组输入时处理标量、一维或二维数组的批次。例如,通常接受(二维)矩阵的线性代数函数可能接受形状为 (4, 3, 2) 的数组,它将其解释为四个 3x2 矩阵的批次。在这种情况下,我们称输入的*核心形状*为 (3, 2),而*批次形状*为 (4,)。同样,通常接受(一维)向量的线性代数函数会将 (4, 3, 2) 数组视为 (4, 3) 的向量批次,在这种情况下,输入的*核心形状*为 (2,),而*批次形状*为 (4, 3)。核心形状的长度也称为*核心维度*。在这些情况下,输出的最终形状是输入的批次形状与输出的核心形状连接(即当输入的批次形状为 () 时输出的形状)。更多信息,请参阅 批处理线性运算。
overwrite_* 参数#
许多线性代数函数具有 overwrite_* 参数,用于表明操作允许就地进行并覆盖输入数组(命名约定是数组参数 a 有一个对应的 overwrite_a;b 参数有一个对应的 overwrite_b,依此类推)。
默认情况下,scipy.linalg 函数会保留它们的输入并在内部进行复制。在某些情况下,就地工作可能会提高性能或避免耗尽内存——但建议进行测量以确保您确实在特定工作负载中看到了改进。
但请注意,设置 overwrite_a=True 仅表示它*可能*就地工作;它是否真的这样做取决于是否满足附加要求。对于大多数函数,就地操作要求以下所有条件为真
输入数组的数据类型是 LAPACK 兼容的:只有单精度和双精度浮点数(
float{32,64})和复数(complex{64,128})数组是兼容的;数组是 Fortran 顺序的;
目前只有二维数组(
ndim == 2)是兼容的;批处理数组使用核心形状数组大小的内部缓冲区。但这在未来的 SciPy 版本中可能会改变。
如果违反了任何这些条件,则会进行内部复制,无论 overwrite_a 是 True 还是 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 参数。在绝大多数使用场景中,添加它们属于过早优化。