SciPy-1-12-中文文档-七-
SciPy 1.12 中文文档(七)
scipy.linalg.eigvals
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.eigvals.html#scipy.linalg.eigvals
scipy.linalg.eigvals(a, b=None, overwrite_a=False, check_finite=True, homogeneous_eigvals=False)
从普通或广义特征值问题计算特征值。
查找一般矩阵的特征值:
a vr[:,i] = w[i] b vr[:,i]
参数:
a(M, M) array_like
将计算其特征值和特征向量的复数或实数矩阵。
b(M, M) array_like,可选
广义特征值问题中的右手边矩阵。如果省略,则假定为单位矩阵。
overwrite_a布尔型,可选
是否覆盖 a 中的数据(可能会提高性能)
check_finite布尔型,可选
是否检查输入矩阵只包含有限数。禁用可能会提高性能,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
homogeneous_eigvals布尔型,可选
如果为 True,则以齐次坐标返回特征值。在这种情况下,w
是一个(2, M)数组,以便:
w[1,i] a vr[:,i] = w[0,i] b vr[:,i]
默认为 False。
返回:
w(M,) 或 (2, M) 双精度或复数 ndarray
每个特征值根据其重复次数重复,但不按任何特定顺序。形状为(M,),除非homogeneous_eigvals=True
。
引发:
LinAlgError
如果特征值计算不收敛
另请参见
一般数组的特征值和右特征向量。
对称或厄米矩阵的特征值
对称/厄米带状矩阵的特征值
对称/厄米三对角矩阵的特征值
示例
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[0., -1.], [1., 0.]])
>>> linalg.eigvals(a)
array([0.+1.j, 0.-1.j])
>>> b = np.array([[0., 1.], [1., 1.]])
>>> linalg.eigvals(a, b)
array([ 1.+0.j, -1.+0.j])
>>> a = np.array([[3., 0., 0.], [0., 8., 0.], [0., 0., 7.]])
>>> linalg.eigvals(a, homogeneous_eigvals=True)
array([[3.+0.j, 8.+0.j, 7.+0.j],
[1.+0.j, 1.+0.j, 1.+0.j]])
scipy.linalg.eigh
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.eigh.html#scipy.linalg.eigh
scipy.linalg.eigh(a, b=None, *, lower=True, eigvals_only=False, overwrite_a=False, overwrite_b=False, turbo=<object object>, eigvals=<object object>, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None)
求解复共轭厄米特或实对称矩阵的标准或广义特征值问题。
寻找数组a
的特征值数组w
,并可选地找到数组v
的特征向量,其中b
是正定的,以便对每个特征值λ(w 的第 i 个条目)及其特征向量vi
(v 的第 i 列)满足以下条件:
a @ vi = λ * b @ vi
vi.conj().T @ a @ vi = λ
vi.conj().T @ b @ vi = 1
在标准问题中,假定b
是单位矩阵。
参数:
a(M, M) array_like
将计算其特征值和特征向量的复共轭厄米特或实对称矩阵。
b(M, M) array_like, optional
一个复共轭厄米特或实对称明确正定矩阵。如果省略,则假定单位矩阵。
lowerbool, optional
是否从a
和(如果适用)b
的下三角或上三角中取相关数组数据。(默认:下三角)
eigvals_onlybool, optional
是否仅计算特征值而不计算特征向量。(默认:两者都计算)
subset_by_indexiterable, optional
如果提供,则这个两元素迭代器定义了所需特征值的起始和结束索引(升序且从 0 开始计数)。要返回第二小到第五小的特征值,使用[1, 4]
。[n-3, n-1]
返回最大的三个。仅在“evr”、“evx”和“gvx”驱动器中可用。通过int()
直接转换为整数。
subset_by_valueiterable, optional
如果提供,则这个两元素迭代器定义了半开区间(a, b]
,仅返回这些值之间的特征值。仅在“evr”、“evx”和“gvx”驱动器中可用。使用np.inf
表示无约束的端点。
driverstr, optional
定义应使用哪个 LAPACK 驱动程序。标准问题的有效选项为“ev”、“evd”、“evr”、“evx”,广义问题的有效选项为“gv”、“gvd”、“gvx”。请参阅备注部分。标准问题的默认值为“evr”。对于广义问题,使用“gvd”进行完整设置,“gvx”进行请求的子集案例。
typeint, optional
对于广义问题,此关键字指定要为w
和v
解决的问题类型(只接受 1、2、3 作为可能的输入):
1 => a @ v = w @ b @ v
2 => a @ b @ v = w @ v
3 => b @ a @ v = w @ v
对于标准问题,此关键字被忽略。
overwrite_abool, optional
是否覆盖a
中的数据(可能提高性能)。默认为 False。
overwrite_bbool, optional
是否覆盖b
中的数据(可能提高性能)。默认为 False。
check_finitebool, optional
是否检查输入矩阵仅包含有限数。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,则可能导致问题(崩溃,非终止)。
turbobool, optional, deprecated
自版本 1.5.0 起不建议使用:eigh
关键字参数 turbo 已被 driver=gvd 关键字取代,并将在 SciPy 1.14.0 中删除。
eigvalstuple (lo, hi),可选,已废弃
自版本 1.5.0 起不建议使用:eigh
关键字参数 eigvals
已被 subset_by_index 关键字取代,并将在 SciPy 1.14.0 中删除。
Returns:
w(N,) ndarray
选择的 N (N<=M) 个特征值,按升序排列,根据其重复次数重复。
v(M, N) ndarray
对称/Hermitian 三对角矩阵的归一化特征向量对应于特征值 w[i]
的列 v[:,i]
。仅当 eigvals_only=False
时返回。
Raises:
LinAlgError
如果特征值计算不收敛,发生错误或 b 矩阵不是正定的。请注意,如果输入矩阵不对称或 Hermitian,则不会报告错误,但结果将是错误的。
See also
eigvalsh
对称或 Hermitian 数组的特征值
eig
非对称数组的特征值和右特征向量
eigh_tridiagonal
对称/Hermitian 三对角矩阵的特征值和右特征向量
Notes
为了允许表示仅具有其上/下三角部分的数组,此函数不会检查输入数组是否为 Hermitian/对称。还要注意,尽管不考虑,但有限性检查适用于整个数组,并且不受“lower”关键字的影响。
此函数在所有可能的关键字组合中使用 LAPACK 驱动程序进行计算,如果数组是实数,则以 sy
为前缀,如果是复数,则以 he
为前缀。例如,使用 “evr” 驱动程序求解浮点数组的问题将通过 “syevr” 解决,使用 “gvx” 驱动程序求解复数数组的问题将通过 “hegvx” 解决等等。
简而言之,最慢且最稳健的驱动程序是经典的 <sy/he>ev
,它使用对称 QR。对于最一般的情况,<sy/he>evr
被视为最佳选择。然而,有些情况下,<sy/he>evd
在更多内存使用的情况下计算更快。<sy/he>evx
,虽然比 <sy/he>ev
快,但在大数组中请求很少的特征值时性能通常比其他情况差,尽管仍然没有性能保证。
对于广义问题,根据给定类型参数进行归一化:
type 1 and 3 : v.conj().T @ a @ v = w
type 2 : inv(v).conj().T @ a @ inv(v) = w
type 1 or 2 : v.conj().T @ b @ v = I
type 3 : v.conj().T @ inv(b) @ v = I
Examples
>>> import numpy as np
>>> from scipy.linalg import eigh
>>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
>>> w, v = eigh(A)
>>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
True
仅请求特征值
>>> w = eigh(A, eigvals_only=True)
请求小于 10 的特征值。
>>> A = np.array([[34, -4, -10, -7, 2],
... [-4, 7, 2, 12, 0],
... [-10, 2, 44, 2, -19],
... [-7, 12, 2, 79, -34],
... [2, 0, -19, -34, 29]])
>>> eigh(A, eigvals_only=True, subset_by_value=[-np.inf, 10])
array([6.69199443e-07, 9.11938152e+00])
请求第二小的特征值及其特征向量
>>> w, v = eigh(A, subset_by_index=[1, 1])
>>> w
array([9.11938152])
>>> v.shape # only a single column is returned
(5, 1)
scipy.linalg.eigvalsh
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.eigvalsh.html#scipy.linalg.eigvalsh
scipy.linalg.eigvalsh(a, b=None, *, lower=True, overwrite_a=False, overwrite_b=False, turbo=<object object>, eigvals=<object object>, type=1, check_finite=True, subset_by_index=None, subset_by_value=None, driver=None)
解决复共轭厄米特或实对称矩阵的标准或广义特征值问题。
查找数组 a 的特征值数组 w,其中 b 是正定的,使得每个特征值λ(w 的第 i 个条目)及其特征向量 vi(v 的第 i 列)满足:
a @ vi = λ * b @ vi
vi.conj().T @ a @ vi = λ
vi.conj().T @ b @ vi = 1
在标准问题中,假定 b 为单位矩阵。
参数:
a(M, M)数组类型
将计算其特征值的复共轭厄米特或实对称矩阵。
b(M, M)数组类型,可选
复共轭厄米特或实对称正定矩阵 a。如果省略,则假定为单位矩阵。
lowerbool, optional
是否从 a 和(如果适用)b 的下三角形或上三角形获取相关数组数据。(默认值:lower)
overwrite_abool,可选
是否覆盖 a 中的数据(可能会提高性能)。默认为 False。
overwrite_bbool,可选
是否覆盖 b 中的数据(可能会提高性能)。默认为 False。
typeint, optional
对于广义问题,此关键字指定要为 w 和 v 解决的问题类型(只接受 1、2、3 作为可能的输入):
1 => a @ v = w @ b @ v
2 => a @ b @ v = w @ v
3 => b @ a @ v = w @ v
此关键字在标准问题中被忽略。
check_finitebool,可选
是否检查输入矩阵仅包含有限数字。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃、非终止)。
subset_by_indexiterable, optional
如果提供,此两元素可迭代对象定义所需特征值的起始和结束索引(升序和从 0 开始索引)。例如,返回第二小到第五小的特征值,使用[1, 4]
。返回最大的三个特征值使用[n-3, n-1]
。仅适用于“evr”、“evx”和“gvx”驱动程序。这些条目通过int()
直接转换为整数。
subset_by_value可迭代对象,可选
如果提供,此两元素可迭代对象定义半开区间‘(a, b]’,如果有的话,仅返回介于这些值之间的特征值。仅适用于“evr”、“evx”和“gvx”驱动程序。使用np.inf
表示无约束端点。
driver字符串,可选
定义应使用哪个 LAPACK 驱动程序。有效选项为“ev”、“evd”、“evr”、“evx”(标准问题)和“gv”、“gvd”、“gvx”(广义问题,其中 b 不为 None)。参见scipy.linalg.eigh
的注释部分。
turbobool,可选,已弃用
自 1.5.0 版本起已弃用:‘eigvalsh’关键字参数turbo已弃用,推荐使用driver=gvd
选项,并将在 SciPy 1.14.0 中移除。
eigvals元组(lo,hi),可选
自版本 1.5.0 起不推荐使用:‘eigvalsh’关键字参数eigvals
已被废弃,建议使用subset_by_index选项,将在 SciPy 1.14.0 中移除。
返回结果:
w(N,) ndarray
N(N<=M)个选定特征值,按升序排列,每个特征值根据其重复次数重复。
引发:
LinAlgError
如果特征值计算不收敛,发生错误或 b 矩阵不是正定的。请注意,如果输入矩阵不对称或不是厄米矩阵,不会报告错误,但结果将是错误的。
另请参阅
eigh
对称/厄米矩阵的特征值和右特征向量
eigvals
一般数组的特征值
eigvals_banded
对称/厄米带状矩阵的特征值
eigvalsh_tridiagonal
对称/厄米三对角矩阵的特征值
注意事项
为了允许表示仅具有其上/下三角部分的数组,此函数不会检查输入数组是否为厄米矩阵/对称矩阵。
这个函数作为scipy.linalg.eigh
的一个一行缩写,选项eigvals_only=True
用于获取特征值而不是特征向量。这里保留它作为一个传统便利功能。使用主函数可以有更多的控制,并且更符合 Python 风格。
示例
更多示例请参见scipy.linalg.eigh
。
>>> import numpy as np
>>> from scipy.linalg import eigvalsh
>>> A = np.array([[6, 3, 1, 5], [3, 0, 5, 1], [1, 5, 6, 2], [5, 1, 2, 2]])
>>> w = eigvalsh(A)
>>> w
array([-3.74637491, -0.76263923, 6.08502336, 12.42399079])
scipy.linalg.eig_banded
scipy.linalg.eig_banded(a_band, lower=False, eigvals_only=False, overwrite_a_band=False, select='a', select_range=None, max_ev=0, check_finite=True)
求解实对称或复共轭厄米特带矩阵的特征值问题。
找到矩阵 a 的特征值 w 和可选的右特征向量 v:
a v[:,i] = w[i] v[:,i]
v.H v = identity
矩阵 a 以下带或上带排序形式存储在 a_band 中:
如果为上三角形式,则 a_band[u + i - j, j] == a[i,j](如果 i <= j);如果为下三角形式,则 a_band[ i - j, j] == a[i,j](如果 i >= j)。
u 是对角线上方的波段数量。
例如 a_band(a 的形状为(6,6),u=2):
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
用*标记的单元格未使用。
参数:
a_band(u+1, M) 类似数组
MxM 矩阵 a 的波段。
lowerbool, optional
矩阵是否以下带形式存储(默认为上带形式)。
eigvals_onlybool, optional
仅计算特征值而不计算特征向量。(默认:也计算特征向量)
overwrite_a_bandbool, optional
丢弃 a_band 中的数据(可能提升性能)。
select,可选
计算哪些特征值
select | 计算 |
---|---|
‘a’ | 所有特征值 |
‘v’ | 区间(min, max]内的特征值 |
‘i’ | 索引 min <= i <= max 的特征值 |
select_range(min, max),可选
选择的特征值范围。
max_evint, optional
对于 select==’v’,预期最大特征值数。对于 select 的其他值,无意义。
如果有疑问,请不要改动此参数。
check_finitebool, optional
是否检查输入矩阵仅包含有限数。禁用可能会提高性能,但如果输入包含无限或 NaN,可能会导致问题(崩溃、非终止)。
返回:
w(M,) ndarray
特征值按升序排列,每个按其重复次数重复。
v(M, M) float or complex ndarray
对应于特征值 w[i]的归一化特征向量是列 v[:,i]。仅当eigvals_only=False
时才返回。
引发:
LinAlgError
如果特征值计算不收敛。
参见
对称/厄米特带矩阵的特征值。
一般数组的特征值和右特征向量。
对称/厄米特阵列的特征值和右特征向量。
对称/厄米特三对角矩阵的特征值和右特征向量。
示例
>>> import numpy as np
>>> from scipy.linalg import eig_banded
>>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
>>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
>>> w, v = eig_banded(Ab, lower=True)
>>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
True
>>> w = eig_banded(Ab, lower=True, eigvals_only=True)
>>> w
array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
仅请求介于[-3, 4]
之间的特征值。
>>> w, v = eig_banded(Ab, lower=True, select='v', select_range=[-3, 4])
>>> w
array([-2.22987175, 3.95222349])
scipy.linalg.eigvals_banded
scipy.linalg.eigvals_banded(a_band, lower=False, overwrite_a_band=False, select='a', select_range=None, check_finite=True)
解决实对称或复厄米特带矩阵特征值问题。
查找矩阵 a 的特征值 w:
a v[:,i] = w[i] v[:,i]
v.H v = identity
矩阵 a 存储在 a_band 中,可以是下三角或上三角顺序:
a_band[u + i - j, j] == a[i,j](如果为上三角形式;i <= j) a_band[ i - j, j] == a[i,j](如果为下三角形式;i >= j)
其中 u 是对角线上方带的数量。
a_band 的示例(a 的形状为 (6,6),u=2):
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
标有 * 的单元格未被使用。
参数:
a_band(u+1, M) array_like
M × M 矩阵 a 的带。
lowerbool, 可选
矩阵是下三角形式。(默认为上三角形式)
overwrite_a_bandbool, 可选
丢弃 a_band 中的数据(可能提高性能)
select, 可选
要计算的特征值
select | 计算 |
---|---|
‘a’ | 所有特征值 |
‘v’ | 特征值在区间 (min, max] 内 |
‘i’ | 特征值在 min <= i <= max 的索引处 |
select_range(min, max), 可选
选择特征值的范围
check_finitebool, 可选
是否检查输入矩阵仅包含有限数字。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
w(M,) ndarray
特征值,按升序排列,每个按其重数重复。
异常:
LinAlgError
如果特征值计算不收敛。
另请参阅
eig_banded
对称/厄米特带矩阵的特征值和右特征向量
eigvalsh_tridiagonal
对称/厄米特三对角矩阵的特征值
eigvals
一般数组的特征值
eigh
对称/厄米特数组的特征值和右特征向量
eig
非对称数组的特征值和右特征向量
示例
>>> import numpy as np
>>> from scipy.linalg import eigvals_banded
>>> A = np.array([[1, 5, 2, 0], [5, 2, 5, 2], [2, 5, 3, 5], [0, 2, 5, 4]])
>>> Ab = np.array([[1, 2, 3, 4], [5, 5, 5, 0], [2, 2, 0, 0]])
>>> w = eigvals_banded(Ab, lower=True)
>>> w
array([-4.26200532, -2.22987175, 3.95222349, 12.53965359])
scipy.linalg.eigh_tridiagonal
scipy.linalg.eigh_tridiagonal(d, e, eigvals_only=False, select='a', select_range=None, check_finite=True, tol=0.0, lapack_driver='auto')
求解实对称三对角矩阵的特征值问题。
查找矩阵a
的特征值w和可选的右特征向量v:
a v[:,i] = w[i] v[:,i]
v.H v = identity
对于具有对角元素d和非对角元素e的实对称矩阵a
。
参数:
dndarray,形状(ndim,)
数组的对角元素。
endarray,形状(ndim-1,)
数组的非对角元素。
eigvals_onlybool,可选
仅计算特征值,不计算特征向量。(默认:同时计算特征向量)
select,可选
要计算的特征值
select | 计算 |
---|---|
‘a’ | 所有特征值 |
‘v’ | 特征值在区间(min, max]内 |
‘i’ | 特征值满足 min <= i <= max 的条件 |
select_range(最小值, 最大值), 可选
选择的特征值范围
check_finitebool,可选
是否检查输入矩阵仅包含有限数值。禁用此选项可能提升性能,但如果输入包含无穷大或 NaN 可能会导致问题(崩溃、不终止)。
tolfloat
每个特征值所需的绝对容差(仅在‘stebz’为LAPACK 驱动器时使用)。如果 <= 0.(默认),则使用机器精度 eps 乘以矩阵a
的 1-范数,其中 eps 是机器精度,|a|
是矩阵a
的 1-范数。
lapack_driverstr
LAPACK 函数的使用,可以是‘auto’,‘stemr’,‘stebz’,‘sterf’或‘stev’。当‘auto’(默认)时,如果select='a'
,则使用‘stemr’,否则使用‘stebz’来找特征值。当使用‘stebz’来找特征值且eigvals_only=False
时,会调用第二次 LAPACK 函数(?STEIN
)来找对应的特征向量。只有当eigvals_only=True
且select='a'
时才能使用‘sterf’。只有当select='a'
时才能使用‘stev’。
返回:
w(M,) 数组
特征值按升序排列,每个根据其重复次数重复。
v(M, M) 数组
与特征值w[i]
对应的归一化特征向量是列v[:,i]
。仅当eigvals_only=False
时返回。
Raises:
LinAlgError
如果特征值计算不收敛。
另见
对称/Hermitian 三对角矩阵的特征值
非对称数组的特征值和右特征向量
对称/Hermitian 数组的特征值和右特征向量
eig_banded
对称/Hermitian 带状矩阵的特征值和右特征向量
注意事项
此函数利用了 LAPACK S/DSTEMR
例程。
示例
>>> import numpy as np
>>> from scipy.linalg import eigh_tridiagonal
>>> d = 3*np.ones(4)
>>> e = -1*np.ones(3)
>>> w, v = eigh_tridiagonal(d, e)
>>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
>>> np.allclose(A @ v - v @ np.diag(w), np.zeros((4, 4)))
True
scipy.linalg.eigvalsh_tridiagonal
scipy.linalg.eigvalsh_tridiagonal(d, e, select='a', select_range=None, check_finite=True, tol=0.0, lapack_driver='auto')
解实对称三对角矩阵的特征值问题。
计算 a
的特征值 w:
a v[:,i] = w[i] v[:,i]
v.H v = identity
对于实对称矩阵 a
,其对角元素为 d,非对角元素为 e。
参数:
dndarray,形状为 (ndim,)
数组的对角元素。
endarray,形状为 (ndim-1,)
数组的非对角元素。
select,可选
要计算的特征值
select | 计算的 |
---|---|
‘a’ | 所有特征值 |
‘v’ | 区间 (min, max] 中的特征值 |
‘i’ | 具有指数 min <= i <= max 的特征值 |
select_range(min, max),可选
选择的特征值范围
check_finitebool,可选
是否检查输入矩阵是否仅包含有限数。禁用可能会带来性能提升,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃、非终止)。
tolfloat
每个特征值所需的绝对容差(仅在 lapack_driver='stebz'
时使用)。如果一个特征值(或簇)位于此宽度的区间内,则认为其已收敛。如果 <= 0(默认),则使用 eps*|a|
的值,其中 eps 是机器精度,而 |a|
是矩阵 a
的 1-范数。
lapack_driverstr
要使用的 LAPACK 函数,可以是 ‘auto’、‘stemr’、‘stebz’、‘sterf’ 或 ‘stev’。当 select='a'
时,默认使用 ‘stemr’,否则使用 ‘stebz’。只有当 select='a'
时才能使用 ‘sterf’ 和 ‘stev’。
返回:
w(M,) ndarray
按升序排列的特征值,每个根据其重数重复。
引发:
LinAlgError
如果特征值计算不收敛。
另请参见
eigh_tridiagonal
对称/Hermitian 三对角矩阵的特征值和右特征向量
示例
>>> import numpy as np
>>> from scipy.linalg import eigvalsh_tridiagonal, eigvalsh
>>> d = 3*np.ones(4)
>>> e = -1*np.ones(3)
>>> w = eigvalsh_tridiagonal(d, e)
>>> A = np.diag(d) + np.diag(e, k=1) + np.diag(e, k=-1)
>>> w2 = eigvalsh(A) # Verify with other eigenvalue routines
>>> np.allclose(w - w2, np.zeros(4))
True
scipy.linalg.lu
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.lu.html#scipy.linalg.lu
scipy.linalg.lu(a, permute_l=False, overwrite_a=False, check_finite=True, p_indices=False)
计算带有部分枢轴的矩阵的 LU 分解。
分解满足:
A = P @ L @ U
其中P
是一个排列矩阵,L
是具有单位对角线元素的下三角矩阵,U
是上三角矩阵。如果将permute_l设置为True
,则L
已排列并且满足A = L @ U
。
参数:
a(M, N) array_like
要分解的数组
permute_lbool,可选
执行乘法 P*L(默认情况下不排列)
overwrite_abool,可选
是否覆盖数据中的数据(可能提高性能)
check_finitebool,可选
是否检查输入矩阵仅包含有限数值。禁用可能会提高性能,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
p_indicesbool,可选
如果为True
,则返回排列信息作为行索引。出于向后兼容性的原因,默认为False
。
返回:
(如果 permute_l
是 False
)
p(…, M, M) ndarray
取决于p_indices的排列数组或向量
l(…, M, K) ndarray
具有单位对角线的下三角或梯形数组。K = min(M, N)
u(…, K, N) ndarray
上三角或梯形数组
(如果 permute_l
是 True
)
pl(…, M, K) ndarray
排列后的 L 矩阵。K = min(M, N)
u(…, K, N) ndarray
上三角或梯形数组
注释
排列矩阵成本高昂,因为它们只是L
的行重新排序,因此强烈建议使用索引,如果需要排列。在 2D 情况下,关系简单地变成A = L[P, :] @ U
。在更高维度中,最好使用permute_l以避免复杂的索引技巧。
在 2D 情况下,如果出于某种原因需要索引,则仍然需要排列矩阵,可以通过np.eye(M)[P, :]
构造。
示例
>>> import numpy as np
>>> from scipy.linalg import lu
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> p, l, u = lu(A)
>>> np.allclose(A, p @ l @ u)
True
>>> p # Permutation matrix
array([[0., 1., 0., 0.], # Row index 1
[0., 0., 0., 1.], # Row index 3
[1., 0., 0., 0.], # Row index 0
[0., 0., 1., 0.]]) # Row index 2
>>> p, _, _ = lu(A, p_indices=True)
>>> p
array([1, 3, 0, 2]) # as given by row indices above
>>> np.allclose(A, l[p, :] @ u)
True
我们也可以使用 nd 数组,例如,一个 4D 数组的演示:
>>> rng = np.random.default_rng()
>>> A = rng.uniform(low=-4, high=4, size=[3, 2, 4, 8])
>>> p, l, u = lu(A)
>>> p.shape, l.shape, u.shape
((3, 2, 4, 4), (3, 2, 4, 4), (3, 2, 4, 8))
>>> np.allclose(A, p @ l @ u)
True
>>> PL, U = lu(A, permute_l=True)
>>> np.allclose(A, PL @ U)
True
scipy.linalg.lu_factor
scipy.linalg.lu_factor(a, overwrite_a=False, check_finite=True)
计算矩阵的置换 LU 分解。
分解是:
A = P L U
当 P 是一个置换矩阵时,L 是单位对角元的下三角矩阵,而 U 是上三角矩阵。
参数:
a(M, N) 数组样式
要分解的矩阵
overwrite_a布尔型,可选
是否覆盖 A 中的数据(可能提高性能)
check_finite布尔型,可选
是否检查输入矩阵仅包含有限数。禁用可能会提高性能,但如果输入包含无穷大或 NaN,则可能会导致问题(崩溃,非终止)。
返回:
lu(M, N) ndarray
包含 U 在其上三角形中,L 在其下三角形中的矩阵。L 的单位对角线元素未存储。
piv(K,) ndarray
表示置换矩阵 P 的枢轴索引:矩阵的第 i 行与行 piv[i]互换。形状为(K,)
,其中K = min(M, N)
。
另请参见
提供更用户友好的 LU 因子分解格式
使用矩阵的 LU 分解解方程系统
注意事项
这是来自 LAPACK 的*GETRF
例程的包装器。与lu
不同,它将 L 和 U 因子输出到单个数组中,并返回枢轴索引而不是置换矩阵。
虽然底层的*GETRF
例程返回基于 1 的枢轴索引,但lu_factor
返回的piv
数组包含基于 0 的索引。
示例
>>> import numpy as np
>>> from scipy.linalg import lu_factor
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> lu, piv = lu_factor(A)
>>> piv
array([2, 2, 3, 3], dtype=int32)
将 LAPACK 的piv
数组转换为 NumPy 索引并测试置换。
>>> def pivot_to_permutation(piv):
... perm = np.arange(len(piv))
... for i in range(len(piv)):
... perm[i], perm[piv[i]] = perm[piv[i]], perm[i]
... return perm
...
>>> p_inv = pivot_to_permutation(piv)
>>> p_inv
array([2, 0, 3, 1])
>>> L, U = np.tril(lu, k=-1) + np.eye(4), np.triu(lu)
>>> np.allclose(A[p_inv] - L @ U, np.zeros((4, 4)))
True
P L U 中的 P 矩阵由逆置换定义,并且可以使用 argsort 恢复:
>>> p = np.argsort(p_inv)
>>> p
array([1, 3, 0, 2])
>>> np.allclose(A - L[p] @ U, np.zeros((4, 4)))
True
或者:
>>> P = np.eye(4)[p]
>>> np.allclose(A - P @ L @ U, np.zeros((4, 4)))
True
scipy.linalg.lu_solve
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.lu_solve.html#scipy.linalg.lu_solve
scipy.linalg.lu_solve(lu_and_piv, b, trans=0, overwrite_b=False, check_finite=True)
解方程系统 a x = b,给定矩阵 a 的 LU 分解
参数:
(lu, piv)
系数矩阵 a 的因式分解,由 lu_factor 给出。特别是 piv 是从 0 开始的枢轴索引。
b数组
右手边
trans,可选
要解决的系统类型:
trans | system |
---|---|
0 | a x = b |
1 | a^T x = b |
2 | a^H x = b |
overwrite_b布尔型,可选
是否覆盖 b 中的数据(可能提高性能)
check_finite布尔型,可选
是否检查输入矩阵仅包含有限数。禁用可能会带来性能提升,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
x数组
系统的解决方案
另请参阅
LU 分解矩阵
示例
>>> import numpy as np
>>> from scipy.linalg import lu_factor, lu_solve
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> b = np.array([1, 1, 1, 1])
>>> lu, piv = lu_factor(A)
>>> x = lu_solve((lu, piv), b)
>>> np.allclose(A @ x - b, np.zeros((4,)))
True
scipy.linalg.svd
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.svd.html#scipy.linalg.svd
scipy.linalg.svd(a, full_matrices=True, compute_uv=True, overwrite_a=False, check_finite=True, lapack_driver='gesdd')
奇异值分解。
将矩阵a因子分解为两个单位矩阵U
和Vh
,以及奇异值(实数、非负)的一维数组s
,使得a == U @ S @ Vh
,其中S
是具有主对角线s
的适当形状的零矩阵。
参数:
a:(M, N)的 array_like
要分解的矩阵。
full_matrices:bool,可选
如果为 True(默认),U和Vh的形状为(M, M)
,(N, N)
。如果为 False,则形状为(M, K)
和(K, N)
,其中K = min(M, N)
。
compute_uv:bool,可选
是否计算U
和Vh
以及s
。默认为 True。
overwrite_a:bool,可选
是否覆盖a;可能提高性能。默认为 False。
check_finite:bool,可选
是否检查输入矩阵只包含有限数。禁用可能提高性能,但如果输入包含无穷大或 NaN,则可能导致问题(崩溃、不终止)。
lapack_driver:{‘gesdd’, ‘gesvd’},可选
是否使用更高效的分而治之方法('gesdd'
)或一般的矩形方法('gesvd'
)来计算 SVD。MATLAB 和 Octave 使用'gesvd'
方法。默认为'gesdd'
。
0.18 版中的新功能。
返回:
U:ndarray
单位矩阵,左奇异向量作为列。形状为(M, M)
或(M, K)
,取决于full_matrices。
s:ndarray
奇异值,按非增顺序排序。形状为(K,),其中K = min(M, N)
。
Vh:ndarray
单位矩阵,右奇异向量作为行。形状为(N, N)
或(K, N)
,取决于full_matrices。
对于compute_uv=False
,仅返回s
。
引发:
LinAlgError
如果奇异值分解计算不收敛。
另请参阅
计算矩阵的奇异值。
构造 Sigma 矩阵,给定向量 s。
示例
>>> import numpy as np
>>> from scipy import linalg
>>> rng = np.random.default_rng()
>>> m, n = 9, 6
>>> a = rng.standard_normal((m, n)) + 1.j*rng.standard_normal((m, n))
>>> U, s, Vh = linalg.svd(a)
>>> U.shape, s.shape, Vh.shape
((9, 9), (6,), (6, 6))
从分解重建原始矩阵:
>>> sigma = np.zeros((m, n))
>>> for i in range(min(m, n)):
... sigma[i, i] = s[i]
>>> a1 = np.dot(U, np.dot(sigma, Vh))
>>> np.allclose(a, a1)
True
或者,使用full_matrices=False
(注意此时U
的形状为(m, n)
而不是(m, m)
):
>>> U, s, Vh = linalg.svd(a, full_matrices=False)
>>> U.shape, s.shape, Vh.shape
((9, 6), (6,), (6, 6))
>>> S = np.diag(s)
>>> np.allclose(a, np.dot(U, np.dot(S, Vh)))
True
>>> s2 = linalg.svd(a, compute_uv=False)
>>> np.allclose(s, s2)
True
scipy.linalg.svdvals
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.svdvals.html#scipy.linalg.svdvals
scipy.linalg.svdvals(a, overwrite_a=False, check_finite=True)
计算矩阵的奇异值。
参数:
a(M, N) 数组样式
要分解的矩阵。
overwrite_abool,可选
是否覆盖 a;可能会提高性能。默认为 False。
check_finitebool,可选
是否检查输入矩阵仅包含有限数值。禁用此选项可能提高性能,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
s(min(M, N),) ndarray
按降序排序的奇异值。
引发:
LinAlgError
如果 SVD 计算不收敛。
另请参阅
svd
计算矩阵的完全奇异值分解。
diagsvd
根据向量 s 构造 Sigma 矩阵。
注意
svdvals(a)
与 svd(a, compute_uv=False)
的唯一区别在于对空矩阵 a
的边缘情况处理,它返回一个空序列:
>>> import numpy as np
>>> a = np.empty((0, 2))
>>> from scipy.linalg import svdvals
>>> svdvals(a)
array([], dtype=float64)
示例
>>> import numpy as np
>>> from scipy.linalg import svdvals
>>> m = np.array([[1.0, 0.0],
... [2.0, 3.0],
... [1.0, 1.0],
... [0.0, 2.0],
... [1.0, 0.0]])
>>> svdvals(m)
array([ 4.28091555, 1.63516424])
我们可以通过计算 m 点乘平面 (x,y) 中所有单位向量 u 的最大长度来验证 m 的最大奇异值。我们用一个大样本近似“所有”单位向量。由于线性性质,我们只需考虑角度在 [0, pi] 内的单位向量。
>>> t = np.linspace(0, np.pi, 2000)
>>> u = np.array([np.cos(t), np.sin(t)])
>>> np.linalg.norm(m.dot(u), axis=0).max()
4.2809152422538475
p 是一个秩为 1 的投影矩阵。在精确算术中,它的奇异值将为 [1, 0, 0, 0]。
>>> v = np.array([0.1, 0.3, 0.9, 0.3])
>>> p = np.outer(v, v)
>>> svdvals(p)
array([ 1.00000000e+00, 2.02021698e-17, 1.56692500e-17,
8.15115104e-34])
正交矩阵的奇异值都为 1。在这里,我们通过使用 scipy.stats.ortho_group
的 rvs() 方法创建一个随机正交矩阵。
>>> from scipy.stats import ortho_group
>>> orth = ortho_group.rvs(4)
>>> svdvals(orth)
array([ 1., 1., 1., 1.])
scipy.linalg.diagsvd
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.diagsvd.html#scipy.linalg.diagsvd
scipy.linalg.diagsvd(s, M, N)
从奇异值和大小 M、N 构造 SVD 中的 sigma 矩阵。
参数:
s(M,) 或 (N,) array_like
奇异值
M整数
矩阵其奇异值为s的大小。
N整数
矩阵其奇异值为s的大小。
返回:
S(M, N) ndarray
在奇异值分解中的 S 矩阵
另请参阅
矩阵的奇异值分解
计算矩阵的奇异值。
示例
>>> import numpy as np
>>> from scipy.linalg import diagsvd
>>> vals = np.array([1, 2, 3]) # The array representing the computed svd
>>> diagsvd(vals, 3, 4)
array([[1, 0, 0, 0],
[0, 2, 0, 0],
[0, 0, 3, 0]])
>>> diagsvd(vals, 4, 3)
array([[1, 0, 0],
[0, 2, 0],
[0, 0, 3],
[0, 0, 0]])
scipy.linalg.orth
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.orth.html#scipy.linalg.orth
scipy.linalg.orth(A, rcond=None)
使用 SVD 构造 A 的范围的正交基
参数:
A(M, N) 类似数组
输入数组
rcond 浮点数,可选
相对条件数。奇异值s
小于rcond * max(s)
被视为零。默认值:浮点数 eps * max(M,N)。
返回:
Q(M, K) ndarray
A 的范围的正交基。K = 由 rcond 确定的 A 的有效秩
另请参见
svd
矩阵的奇异值分解
null_space
矩阵的零空间
示例
>>> import numpy as np
>>> from scipy.linalg import orth
>>> A = np.array([[2, 0, 0], [0, 5, 0]]) # rank 2 array
>>> orth(A)
array([[0., 1.],
[1., 0.]])
>>> orth(A.T)
array([[0., 1.],
[1., 0.],
[0., 0.]])
scipy.linalg.null_space
scipy.linalg.null_space(A, rcond=None)
使用 SVD 构造 A 的零空间的标准正交基
参数:
A(M, N) array_like
输入数组
rcondfloat, optional
相对条件数。比rcond * max(s)
小的奇异值s
被认为是零。默认值:浮点数 eps * max(M, N)。
返回:
Z(N, K) ndarray
A 的零空间的标准正交基。K = 有效零空间的维度,由 rcond 确定。
另请参阅
svd
矩阵的奇异值分解
orth
矩阵范围
示例
1-D 零空间:
>>> import numpy as np
>>> from scipy.linalg import null_space
>>> A = np.array([[1, 1], [1, 1]])
>>> ns = null_space(A)
>>> ns * np.copysign(1, ns[0,0]) # Remove the sign ambiguity of the vector
array([[ 0.70710678],
[-0.70710678]])
2-D 零空间:
>>> from numpy.random import default_rng
>>> rng = default_rng()
>>> B = rng.random((3, 5))
>>> Z = null_space(B)
>>> Z.shape
(5, 2)
>>> np.allclose(B.dot(Z), 0)
True
基向量是标准正交的(舍入误差):
>>> Z.T.dot(Z)
array([[ 1.00000000e+00, 6.92087741e-17],
[ 6.92087741e-17, 1.00000000e+00]])
scipy.linalg.ldl
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.ldl.html#scipy.linalg.ldl
scipy.linalg.ldl(A, lower=True, hermitian=True, overwrite_a=False, check_finite=True)
计算对称/ Hermitian 矩阵的 LDLt 或 Bunch-Kaufman 分解。
此函数返回一个块对角矩阵 D,其中每个块的大小最多为 2x2,并且可能会返回一个可能排列的单位下三角矩阵 L
,使得分解 A = L D L^H
或 A = L D L^T
成立。如果 lower 为 False,则返回(再次可能排列的)上三角矩阵作为外因子。
排列数组可以通过行洗牌简单地将外因子三角化,即 lu[perm, :]
是一个上/下三角矩阵。这也等同于与置换矩阵 P
的乘积 P.dot(lu)
,其中 P
是列置换的单位矩阵 I[:, perm]
。
根据布尔值 lower 的值,仅引用输入数组的上三角或下三角部分。因此,输入一个三角矩阵会得到与提供完整矩阵相同的结果。
参数:
A:array_like
方阵输入数组
lower:bool, 可选
这会在因子分解的下三角或上三角外因子之间切换。下三角(lower=True
)是默认值。
hermitian:bool, 可选
对于复数数组,这会定义是否假设 A = A.conj().T
或 A = A.T
。对于实数数组,此切换无效。
overwrite_a:bool, 可选
允许重写 A 中的数据(可能会提升性能)。默认值为 False。
check_finite:bool, 可选
是否检查输入矩阵仅包含有限数。禁用可能会带来性能提升,但如果输入包含无穷大或 NaN,则可能导致问题(崩溃、非终止)。
返回:
lu:LU 分解后的数组
因子分解的(可能)排列的上/下三角外因子。
d:数组 d
因子分解的块对角乘积。
perm:数组 perm
将 lu 变为三角形形式的行置换索引数组。
异常:
值错误
如果输入数组不是方阵。
ComplexWarning
如果给定一个具有非零虚部对角线的复数数组,并且 hermitian 设置为 True。
另请参见
注意
该函数使用来自 LAPACK 的对称矩阵的 ?SYTRF
例程和 Hermitian 矩阵的 ?HETRF
例程。详见 [1] 获取算法细节。
根据 lower 关键字的值,只引用输入数组的下三角或上三角部分。此关键字还定义了因子分解的外因子的结构。
新版本 1.1.0 中引入。
参考文献
[1]
J.R. Bunch, L. Kaufman, 计算惯性和解决对称线性系统的一些稳定方法, Math. Comput. Vol.31, 1977. DOI:10.2307/2005787
Examples
给定一个代表带有其条目的完整对称数组的上三角数组 a
,获取 l
,d
和置换向量 perm:
>>> import numpy as np
>>> from scipy.linalg import ldl
>>> a = np.array([[2, -1, 3], [0, 2, 0], [0, 0, 1]])
>>> lu, d, perm = ldl(a, lower=0) # Use the upper part
>>> lu
array([[ 0\. , 0\. , 1\. ],
[ 0\. , 1\. , -0.5],
[ 1\. , 1\. , 1.5]])
>>> d
array([[-5\. , 0\. , 0\. ],
[ 0\. , 1.5, 0\. ],
[ 0\. , 0\. , 2\. ]])
>>> perm
array([2, 1, 0])
>>> lu[perm, :]
array([[ 1\. , 1\. , 1.5],
[ 0\. , 1\. , -0.5],
[ 0\. , 0\. , 1\. ]])
>>> lu.dot(d).dot(lu.T)
array([[ 2., -1., 3.],
[-1., 2., 0.],
[ 3., 0., 1.]])
scipy.linalg.cholesky
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.cholesky.html#scipy.linalg.cholesky
scipy.linalg.cholesky(a, lower=False, overwrite_a=False, check_finite=True)
计算矩阵的乔列斯基分解。
返回埃尔米特正定矩阵 A 的乔列斯基分解,(A = L L^) 或 (A = U^ U)。
参数:
a(M, M) array_like
要分解的矩阵
lowerbool, 可选
是否计算上三角或下三角的乔列斯基分解。默认为上三角。
overwrite_abool, 可选
是否覆盖a中的数据(可能提高性能)。
check_finitebool, 可选
是否检查输入矩阵是否仅包含有限数。禁用可能提高性能,但如果输入确实包含无穷大或 NaN,则可能导致问题(崩溃、非终止)。
返回:
c(M, M) ndarray
a的上三角或下三角乔列斯基因子。
引发:
LinAlgError如果分解失败。
示例
>>> import numpy as np
>>> from scipy.linalg import cholesky
>>> a = np.array([[1,-2j],[2j,5]])
>>> L = cholesky(a, lower=True)
>>> L
array([[ 1.+0.j, 0.+0.j],
[ 0.+2.j, 1.+0.j]])
>>> L @ L.T.conj()
array([[ 1.+0.j, 0.-2.j],
[ 0.+2.j, 5.+0.j]])
scipy.linalg.cholesky_banded
scipy.linalg.cholesky_banded(ab, overwrite_ab=False, lower=False, check_finite=True)
Cholesky 分解一个带状 Hermitian 正定矩阵
矩阵 a 以 lower-diagonal 或 upper-diagonal 形式存储在 ab 中:
ab[u + i - j, j] == a[i,j] (if upper form; i <= j)
ab[ i - j, j] == a[i,j] (if lower form; i >= j)
ab 的示例(a 的形状为(6,6),u=2):
upper form:
* * a02 a13 a24 a35
* a01 a12 a23 a34 a45
a00 a11 a22 a33 a44 a55
lower form:
a00 a11 a22 a33 a44 a55
a10 a21 a32 a43 a54 *
a20 a31 a42 a53 * *
参数:
ab(u + 1, M) array_like
带状矩阵
overwrite_abbool, 可选
在 ab 中丢弃数据(可能增强性能)
lowerbool, 可选
矩阵是否以 lower 形式表示(默认为 upper 形式)
check_finitebool, 可选
是否检查输入矩阵仅包含有限数。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
c(u + 1, M) ndarray
a 的 Cholesky 分解,与 ab 具有相同的带状格式
参见
cho_solve_banded
解决线性方程组,给定一个带状厄米特矩阵的 Cholesky 分解。
示例
>>> import numpy as np
>>> from scipy.linalg import cholesky_banded
>>> from numpy import allclose, zeros, diag
>>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]])
>>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1)
>>> A = A + A.conj().T + np.diag(Ab[2, :])
>>> c = cholesky_banded(Ab)
>>> C = np.diag(c[0, 2:], k=2) + np.diag(c[1, 1:], k=1) + np.diag(c[2, :])
>>> np.allclose(C.conj().T @ C - A, np.zeros((5, 5)))
True
scipy.linalg.cho_factor
scipy.linalg.cho_factor(a, lower=False, overwrite_a=False, check_finite=True)
计算矩阵的 Cholesky 分解,以在 cho_solve
中使用
返回包含 Hermitian 正定矩阵 a 的 Cholesky 分解 A = L L*
或 A = U* U
的矩阵。返回值可以直接用作 cho_solve
的第一个参数。
警告
返回的矩阵还在未使用 Cholesky 分解的条目中包含随机数据。如果需要将这些条目清零,请改用函数 cholesky
。
参数:
a(M, M) 类似数组
要分解的矩阵
lower布尔值,可选
是否计算上三角或下三角的 Cholesky 分解(默认为:上三角)
overwrite_a布尔值,可选
是否覆盖数据(可能会提高性能)
check_finite布尔值,可选
是否检查输入矩阵只包含有限数字。禁用此选项可能会提高性能,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃、非终止)。
返回:
c(M, M) 数组
矩阵的上三角形或下三角形包含矩阵 a 的 Cholesky 因子。矩阵的其他部分包含随机数据。
lower布尔值
指示因子位于下三角形还是上三角形的标志
引发:
线性代数错误
如果分解失败,则引发错误。
另请参阅
cho_solve
使用矩阵的 Cholesky 分解解线性方程组。
示例
>>> import numpy as np
>>> from scipy.linalg import cho_factor
>>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
>>> c, low = cho_factor(A)
>>> c
array([[3\. , 1\. , 0.33333333, 1.66666667],
[3\. , 2.44948974, 1.90515869, -0.27216553],
[1\. , 5\. , 2.29330749, 0.8559528 ],
[5\. , 1\. , 2\. , 1.55418563]])
>>> np.allclose(np.triu(c).T @ np. triu(c) - A, np.zeros((4, 4)))
True
scipy.linalg.cho_solve
scipy.linalg.cho_solve(c_and_lower, b, overwrite_b=False, check_finite=True)
解线性方程 A x = b,给定 A 的 Cholesky 分解。
参数:
(c, lower)元组,(数组,布尔值)
给定 cho_factor 给出的 a 的 Cholesky 分解
b数组
右侧
overwrite_b布尔值,可选
是否覆盖 b 中的数据(可能提高性能)
check_finite布尔值,可选
是否检查输入矩阵仅包含有限数值。禁用可能会提高性能,但如果输入包含无穷大或 NaN,则可能导致问题(崩溃、非终止)。
返回:
x数组
方程组 A x = b 的解
另见
cho_factor
矩阵 a 的 Cholesky 分解
示例
>>> import numpy as np
>>> from scipy.linalg import cho_factor, cho_solve
>>> A = np.array([[9, 3, 1, 5], [3, 7, 5, 1], [1, 5, 9, 2], [5, 1, 2, 6]])
>>> c, low = cho_factor(A)
>>> x = cho_solve((c, low), [1, 1, 1, 1])
>>> np.allclose(A @ x - [1, 1, 1, 1], np.zeros(4))
True
scipy.linalg.cho_solve_banded
scipy.linalg.cho_solve_banded(cb_and_lower, b, overwrite_b=False, check_finite=True)
解线性方程组A x = b
,给定带状 Hermitian 矩阵A
的 Cholesky 分解。
参数:
(cb, lower)元组,(ndarray, bool)
cb是由 cholesky_banded 给出的 A 的 Cholesky 分解。lower必须与传递给 cholesky_banded 的值相同。
b类数组
右侧向量
overwrite_b布尔值,可选
如果为 True,函数将覆盖b中的值。
check_finite布尔值,可选
是否检查输入矩阵仅包含有限数。禁用可能会提高性能,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
x数组
系统 A x = b 的解
另请参见
带状矩阵的 Cholesky 分解
注意事项
自 0.8.0 版本开始。
示例
>>> import numpy as np
>>> from scipy.linalg import cholesky_banded, cho_solve_banded
>>> Ab = np.array([[0, 0, 1j, 2, 3j], [0, -1, -2, 3, 4], [9, 8, 7, 6, 9]])
>>> A = np.diag(Ab[0,2:], k=2) + np.diag(Ab[1,1:], k=1)
>>> A = A + A.conj().T + np.diag(Ab[2, :])
>>> c = cholesky_banded(Ab)
>>> x = cho_solve_banded((c, False), np.ones(5))
>>> np.allclose(A @ x - np.ones(5), np.zeros(5))
True
scipy.linalg.polar
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.polar.html#scipy.linalg.polar
scipy.linalg.polar(a, side='right')
计算极分解。
返回极分解的因子 [1] u 和 p,使得 a = up
(如果 side 是“right”)或 a = pu
(如果 side 是“left”),其中 p 是半正定矩阵。根据 a 的形状,u 的行或列正交。当 a 是方阵时,u 是方酉矩阵。当 a 不是方阵时,计算“标准极分解” [2]。
参数:
a(m, n) array_like
要分解的数组。
side,可选
确定计算右极分解还是左极分解。如果 side 是“right”,那么 a = up
。如果 side 是“left”,那么 a = pu
。默认为“right”。
返回:
u(m, n) ndarray
如果 a 是方阵,则 u 是酉矩阵。如果 m > n,则 a 的列正交;如果 m < n,则 u 的行正交。
pndarray
p 是埃尔米特半正定矩阵。如果 a 非奇异,则 p 是正定的。p 的形状为 (n, n) 或 (m, m),具体取决于 side 是“right” 还是“left”。
参考文献
[1]
R. A. Horn 和 C. R. Johnson,《矩阵分析》,剑桥大学出版社,1985 年。
[2]
N. J. Higham,《矩阵函数:理论与计算》,SIAM,2008 年。
示例
>>> import numpy as np
>>> from scipy.linalg import polar
>>> a = np.array([[1, -1], [2, 4]])
>>> u, p = polar(a)
>>> u
array([[ 0.85749293, -0.51449576],
[ 0.51449576, 0.85749293]])
>>> p
array([[ 1.88648444, 1.2004901 ],
[ 1.2004901 , 3.94446746]])
一个非方阵示例,其中 m < n:
>>> b = np.array([[0.5, 1, 2], [1.5, 3, 4]])
>>> u, p = polar(b)
>>> u
array([[-0.21196618, -0.42393237, 0.88054056],
[ 0.39378971, 0.78757942, 0.4739708 ]])
>>> p
array([[ 0.48470147, 0.96940295, 1.15122648],
[ 0.96940295, 1.9388059 , 2.30245295],
[ 1.15122648, 2.30245295, 3.65696431]])
>>> u.dot(p) # Verify the decomposition.
array([[ 0.5, 1\. , 2\. ],
[ 1.5, 3\. , 4\. ]])
>>> u.dot(u.T) # The rows of u are orthonormal.
array([[ 1.00000000e+00, -2.07353665e-17],
[ -2.07353665e-17, 1.00000000e+00]])
另一个非方阵示例,其中 m > n:
>>> c = b.T
>>> u, p = polar(c)
>>> u
array([[-0.21196618, 0.39378971],
[-0.42393237, 0.78757942],
[ 0.88054056, 0.4739708 ]])
>>> p
array([[ 1.23116567, 1.93241587],
[ 1.93241587, 4.84930602]])
>>> u.dot(p) # Verify the decomposition.
array([[ 0.5, 1.5],
[ 1\. , 3\. ],
[ 2\. , 4\. ]])
>>> u.T.dot(u) # The columns of u are orthonormal.
array([[ 1.00000000e+00, -1.26363763e-16],
[ -1.26363763e-16, 1.00000000e+00]])
scipy.linalg.qr
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.qr.html#scipy.linalg.qr
scipy.linalg.qr(a, overwrite_a=False, lwork=None, mode='full', pivoting=False, check_finite=True)
计算矩阵的 QR 分解。
计算分解A = Q R
,其中 Q 是单位矩阵/正交矩阵,R 是上三角矩阵。
参数:
a(M, N) array_like
要分解的矩阵
overwrite_abool,可选
如果overwrite_a设置为 True,重复使用现有输入数据结构而不是创建新的数据结构,可能会提高性能。
lworkint,可选
工作数组大小,lwork >= a.shape[1]
。如果为 None 或-1,则计算最佳大小。
mode,可选
确定要返回的信息:返回 Q 和 R('full',默认),仅返回 R('r'),或者返回经济型大小计算的 Q 和 R('economic',详见备注)。最后一个选项'raw'(在 SciPy 0.11 中添加)使函数以 LAPACK 使用的内部格式返回两个矩阵(Q,TAU)。
pivotingbool,可选
是否应在用于排名显示 qr 分解的枢轴处理中包括枢轴。如果使用枢轴,则计算分解A P = Q R
,如上所述,但选择 P 使得 R 的对角线非递增。
check_finitebool,可选
是否检查输入矩阵仅包含有限数字。禁用可能会带来性能提升,但如果输入确实包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
Qfloat 或复数 ndarray
形状为(M, M)或者对于mode='economic'
为(M, K)的形状。如果mode='r'
,则不返回。如果mode='raw'
,则由元组(Q, TAU)
替代。
Rfloat 或复数 ndarray
形状为(M, N)或者对于mode in ['economic', 'raw']
为(K, N)。K = min(M, N)
。
Pint ndarray
对于pivoting=True
的形状为(N,)。如果pivoting=False
,则不返回。
引发:
LinAlgError
如果分解失败则引发
备注
这是 LAPACK 例程 dgeqrf、zgeqrf、dorgqr、zungqr、dgeqp3 和 zgeqp3 的接口。
如果mode=economic
,则 Q 和 R 的形状为(M, K)和(K, N),而不是(M,M)和(M,N),其中K=min(M,N)
。
示例
>>> import numpy as np
>>> from scipy import linalg
>>> rng = np.random.default_rng()
>>> a = rng.standard_normal((9, 6))
>>> q, r = linalg.qr(a)
>>> np.allclose(a, np.dot(q, r))
True
>>> q.shape, r.shape
((9, 9), (9, 6))
>>> r2 = linalg.qr(a, mode='r')
>>> np.allclose(r, r2)
True
>>> q3, r3 = linalg.qr(a, mode='economic')
>>> q3.shape, r3.shape
((9, 6), (6, 6))
>>> q4, r4, p4 = linalg.qr(a, pivoting=True)
>>> d = np.abs(np.diag(r4))
>>> np.all(d[1:] <= d[:-1])
True
>>> np.allclose(a[:, p4], np.dot(q4, r4))
True
>>> q4.shape, r4.shape, p4.shape
((9, 9), (9, 6), (6,))
>>> q5, r5, p5 = linalg.qr(a, mode='economic', pivoting=True)
>>> q5.shape, r5.shape, p5.shape
((9, 6), (6, 6), (6,))
scipy.linalg.qr_multiply
scipy.linalg.qr_multiply(a, c, mode='right', pivoting=False, conjugate=False, overwrite_a=False, overwrite_c=False)
计算 QR 分解并将 Q 与矩阵相乘。
计算分解A = Q R
,其中 Q 是单位/正交矩阵,R 是上三角矩阵。将 Q 与向量或矩阵 c 相乘。
参数:
a(M, N),array_like
输入数组
carray_like
要乘以q
的输入数组。
mode,可选
如果 mode 为‘left’,则返回Q @ c
,如果 mode 为‘right’,则返回c @ Q
。如果 mode 为‘left’,则 c 的形状必须适合矩阵乘法,min(a.shape) == c.shape[0]
;如果 mode 为‘right’,则a.shape[0] == c.shape[1]
。
pivotingbool,可选
是否应在 rank-revealing qr 分解中包含枢轴。有关 qr 的文档,请参阅。
conjugatebool,可选
是否应复合 Q。这可能比显式复合更快。
overwrite_abool,可选
数据是否在 a 中覆盖(可能会提高性能)
overwrite_cbool,可选
数据是否被覆盖(可能会提高性能)。如果使用此选项,则 c 必须足够大以保存结果,即如果 mode 为‘left’,则c.shape[0]
=a.shape[0]
。
Returns:
CQndarray
Q
和c
的乘积。
R(K, N),ndarray
结果 QR 分解的 R 数组,其中K = min(M, N)
。
P(N,) ndarray
整数枢轴数组。仅当pivoting=True
时返回。
Raises:
LinAlgError
如果 QR 分解失败,则引发。
Notes
这是 LAPACK 例程?GEQRF
、?ORMQR
、?UNMQR
和?GEQP3
的接口。
版本 0.11.0 中的新功能。
Examples
>>> import numpy as np
>>> from scipy.linalg import qr_multiply, qr
>>> A = np.array([[1, 3, 3], [2, 3, 2], [2, 3, 3], [1, 3, 2]])
>>> qc, r1, piv1 = qr_multiply(A, 2*np.eye(4), pivoting=1)
>>> qc
array([[-1., 1., -1.],
[-1., -1., 1.],
[-1., -1., -1.],
[-1., 1., 1.]])
>>> r1
array([[-6., -3., -5\. ],
[ 0., -1., -1.11022302e-16],
[ 0., 0., -1\. ]])
>>> piv1
array([1, 0, 2], dtype=int32)
>>> q2, r2, piv2 = qr(A, mode='economic', pivoting=1)
>>> np.allclose(2*q2 - qc, np.zeros((4, 3)))
True
scipy.linalg.qr_update
scipy.linalg.qr_update(Q, R, u, v, overwrite_qruv=False, check_finite=True)
排名 k 的 QR 更新
如果A = Q R
是A
的 QR 分解,则返回A + u v**T
的 QR 分解(对于实数A
)或A + u v**H
的 QR 分解(对于复数A
)。
参数:
Q(M, M) 或 (M, N) 类似数组
QR 分解后的酉/正交矩阵。
R(M, N) 或 (N, N) 类似数组
QR 分解后的上三角矩阵。
u(M,) 或 (M, k) 类似数组
左侧更新向量
v(N,) 或 (N, k) 类似数组
右侧更新向量
overwrite_qruv bool,可选
如果为 True,在执行更新时尽可能消耗 Q、R、u 和 v,否则根据需要进行复制。默认为 False。
check_finite bool,可选
是否检查输入矩阵仅包含有限数字。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。默认为 True。
返回:
Q1 ndarray
更新后的酉/正交因子
R1 ndarray
更新后的上三角因子
另见
qr
, qr_multiply
, qr_delete
, qr_insert
注释
此例程不保证R1的对角线条目是实数或正数。
新版本 0.16.0 中新增。
参考文献
[1]
高卢布(G. H.)与范伦(C. F.)·卢恩,《矩阵计算》,第三版(约翰·霍普金斯大学出版社,1996 年)。
[2]
丹尼尔(J. W.)、格拉格(W. B.)、考夫曼(L.)与斯图尔特(G. W.),《重正交化和稳定算法用于更新格拉姆-施密特 QR 分解》,数学与计算 30,772-795 页(1976 年)。
[3]
莱切尔(L.)与格拉格(W. B.),《用于更新 QR 分解的 FORTRAN 子程序》,ACM 数学软件事务 16,369-377 页(1990 年)。
示例
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[ 3., -2., -2.],
... [ 6., -9., -3.],
... [ -3., 10., 1.],
... [ 6., -7., 4.],
... [ 7., 8., -6.]])
>>> q, r = linalg.qr(a)
鉴于此 QR 分解,执行一个排名 1 的更新。
>>> u = np.array([7., -2., 4., 3., 5.])
>>> v = np.array([1., 3., -5.])
>>> q_up, r_up = linalg.qr_update(q, r, u, v, False)
>>> q_up
array([[ 0.54073807, 0.18645997, 0.81707661, -0.02136616, 0.06902409], # may vary (signs)
[ 0.21629523, -0.63257324, 0.06567893, 0.34125904, -0.65749222],
[ 0.05407381, 0.64757787, -0.12781284, -0.20031219, -0.72198188],
[ 0.48666426, -0.30466718, -0.27487277, -0.77079214, 0.0256951 ],
[ 0.64888568, 0.23001 , -0.4859845 , 0.49883891, 0.20253783]])
>>> r_up
array([[ 18.49324201, 24.11691794, -44.98940746], # may vary (signs)
[ 0\. , 31.95894662, -27.40998201],
[ 0\. , 0\. , -9.25451794],
[ 0\. , 0\. , 0\. ],
[ 0\. , 0\. , 0\. ]])
更新等效,但比以下更快。
>>> a_up = a + np.outer(u, v)
>>> q_direct, r_direct = linalg.qr(a_up)
检查我们是否有等价结果:
>>> np.allclose(np.dot(q_up, r_up), a_up)
True
而更新后的 Q 仍然是酉的:
>>> np.allclose(np.dot(q_up.T, q_up), np.eye(5))
True
还可以更新经济(减少、薄)分解:
>>> qe, re = linalg.qr(a, mode='economic')
>>> qe_up, re_up = linalg.qr_update(qe, re, u, v, False)
>>> qe_up
array([[ 0.54073807, 0.18645997, 0.81707661], # may vary (signs)
[ 0.21629523, -0.63257324, 0.06567893],
[ 0.05407381, 0.64757787, -0.12781284],
[ 0.48666426, -0.30466718, -0.27487277],
[ 0.64888568, 0.23001 , -0.4859845 ]])
>>> re_up
array([[ 18.49324201, 24.11691794, -44.98940746], # may vary (signs)
[ 0\. , 31.95894662, -27.40998201],
[ 0\. , 0\. , -9.25451794]])
>>> np.allclose(np.dot(qe_up, re_up), a_up)
True
>>> np.allclose(np.dot(qe_up.T, qe_up), np.eye(3))
True
类似上述,执行一个二阶更新。
>>> u2 = np.array([[ 7., -1,],
... [-2., 4.],
... [ 4., 2.],
... [ 3., -6.],
... [ 5., 3.]])
>>> v2 = np.array([[ 1., 2.],
... [ 3., 4.],
... [-5., 2]])
>>> q_up2, r_up2 = linalg.qr_update(q, r, u2, v2, False)
>>> q_up2
array([[-0.33626508, -0.03477253, 0.61956287, -0.64352987, -0.29618884], # may vary (signs)
[-0.50439762, 0.58319694, -0.43010077, -0.33395279, 0.33008064],
[-0.21016568, -0.63123106, 0.0582249 , -0.13675572, 0.73163206],
[ 0.12609941, 0.49694436, 0.64590024, 0.31191919, 0.47187344],
[-0.75659643, -0.11517748, 0.10284903, 0.5986227 , -0.21299983]])
>>> r_up2
array([[-23.79075451, -41.1084062 , 24.71548348], # may vary (signs)
[ 0\. , -33.83931057, 11.02226551],
[ 0\. , 0\. , 48.91476811],
[ 0\. , 0\. , 0\. ],
[ 0\. , 0\. , 0\. ]])
这个更新也是A + U V**T
的有效 QR 分解。
>>> a_up2 = a + np.dot(u2, v2.T)
>>> np.allclose(a_up2, np.dot(q_up2, r_up2))
True
>>> np.allclose(np.dot(q_up2.T, q_up2), np.eye(5))
True
scipy.linalg.qr_delete
scipy.linalg.qr_delete(Q, R, k, int p=1, which=u'row', overwrite_qr=False, check_finite=True)
行或列删除的 QR 下降
如果A = Q R
是A
的 QR 分解,则返回A
的 QR 分解,其中从行或列k
开始删除p
行或列。
参数:
Q(M, M)或(M, N) array_like
来自 QR 分解的酉/正交矩阵。
R(M, N)或(N, N) array_like
来自 QR 分解的上三角矩阵。
kint
要删除的第一行或列的索引。
pint,可选
要删除的行或列数,默认为 1。
which: {‘row’, ‘col’},可选
确定将删除行或列,默认为‘行’
overwrite_qrbool,可选
如果为 True,消耗 Q 和 R,用它们的下降版本覆盖它们的内容,并返回适当大小的视图。默认为 False。
check_finitebool,可选
是否检查输入矩阵仅包含有限数。禁用可能会带来性能提升,但如果输入确实包含无穷大或 NaN,可能会导致问题(崩溃,非终止)。默认为 True。
返回:
Q1ndarray
更新后的酉/正交因子
R1ndarray
更新后的上三角因子
另见
qr
, qr_multiply
, qr_insert
, qr_update
注意事项
此例程不保证R1
的对角线条目为正。
新版本 0.16.0 中加入。
参考资料
[1]
Golub, G. H. & Van Loan, C. F. Matrix Computations, 3rd Ed. (Johns Hopkins University Press, 1996).
[2]
Daniel, J. W., Gragg, W. B., Kaufman, L. & Stewart, G. W. Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization. Math. Comput. 30, 772-795 (1976).
[3]
Reichel, L. & Gragg, W. B. Algorithm 686: FORTRAN Subroutines for Updating the QR Decomposition. ACM Trans. Math. Softw. 16, 369-377 (1990).
示例
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[ 3., -2., -2.],
... [ 6., -9., -3.],
... [ -3., 10., 1.],
... [ 6., -7., 4.],
... [ 7., 8., -6.]])
>>> q, r = linalg.qr(a)
给定这个 QR 分解,当移除 2 行时更新 q 和 r。
>>> q1, r1 = linalg.qr_delete(q, r, 2, 2, 'row', False)
>>> q1
array([[ 0.30942637, 0.15347579, 0.93845645], # may vary (signs)
[ 0.61885275, 0.71680171, -0.32127338],
[ 0.72199487, -0.68017681, -0.12681844]])
>>> r1
array([[ 9.69535971, -0.4125685 , -6.80738023], # may vary (signs)
[ 0\. , -12.19958144, 1.62370412],
[ 0\. , 0\. , -0.15218213]])
此更新与以下方法等效,但速度更快。
>>> a1 = np.delete(a, slice(2,4), 0)
>>> a1
array([[ 3., -2., -2.],
[ 6., -9., -3.],
[ 7., 8., -6.]])
>>> q_direct, r_direct = linalg.qr(a1)
检查我们是否有等效的结果:
>>> np.dot(q1, r1)
array([[ 3., -2., -2.],
[ 6., -9., -3.],
[ 7., 8., -6.]])
>>> np.allclose(np.dot(q1, r1), a1)
True
更新后的 Q 仍然是酉的:
>>> np.allclose(np.dot(q1.T, q1), np.eye(3))
True
scipy.linalg.qr_insert
scipy.linalg.qr_insert(Q, R, u, k, which='row', rcond=None, overwrite_qru=False, check_finite=True)
QR 更新行或列插入
如果A = Q R
是A
的 QR 分解,则返回在行或列从 k 开始插入的A
的 QR 分解。
参数:
Q(M, M) array_like
A
的 QR 分解的单位/正交矩阵。
R(M, N) array_like
A
的 QR 分解的上三角矩阵。
u(N,), (p, N), (M,), or (M, p) array_like
要插入的行或列
kint
要插入u之前的索引。
which: {‘row’, ‘col’}, optional
决定是否插入行或列,默认为'row'
rcondfloat
Q
增广为u/||u||
的倒数条件数的下限。仅在更新经济模式(薄,(M,N) (N,N))分解时使用。如果为 None,则使用机器精度。默认为 None。
overwrite_qrubool, optional
如果为 True,则在执行更新时尽可能消耗 Q、R 和 u,否则根据需要制作副本。默认为 False。
check_finitebool, optional
是否检查输入矩阵只包含有限数字。禁用可能会带来性能提升,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃,非终止)。默认为 True。
返回:
Q1ndarray
更新后的单位/正交因子
R1ndarray
更新后的上三角因子
Raises:
LinAlgError
如果更新(M,N) (N,N)分解,并且带有 u/||u||增广的 Q 的倒数条件数小于 rcond。
另请参阅
qr
, qr_multiply
, qr_delete
, qr_update
注释
此例程不保证R1
的对角线条目为正。
新版本 0.16.0 中添加。
参考文献
[1]
Golub, G. H. & Van Loan, C. F. Matrix Computations, 3rd Ed. (Johns Hopkins University Press, 1996).
[2]
Daniel, J. W., Gragg, W. B., Kaufman, L. & Stewart, G. W. Reorthogonalization and stable algorithms for updating the Gram-Schmidt QR factorization. Math. Comput. 30, 772-795 (1976).
[3]
Reichel, L. & Gragg, W. B. Algorithm 686: FORTRAN Subroutines for Updating the QR Decomposition. ACM Trans. Math. Softw. 16, 369-377 (1990).
示例
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[ 3., -2., -2.],
... [ 6., -7., 4.],
... [ 7., 8., -6.]])
>>> q, r = linalg.qr(a)
给定此 QR 分解,当插入 2 行时更新 q 和 r。
>>> u = np.array([[ 6., -9., -3.],
... [ -3., 10., 1.]])
>>> q1, r1 = linalg.qr_insert(q, r, u, 2, 'row')
>>> q1
array([[-0.25445668, 0.02246245, 0.18146236, -0.72798806, 0.60979671], # may vary (signs)
[-0.50891336, 0.23226178, -0.82836478, -0.02837033, -0.00828114],
[-0.50891336, 0.35715302, 0.38937158, 0.58110733, 0.35235345],
[ 0.25445668, -0.52202743, -0.32165498, 0.36263239, 0.65404509],
[-0.59373225, -0.73856549, 0.16065817, -0.0063658 , -0.27595554]])
>>> r1
array([[-11.78982612, 6.44623587, 3.81685018], # may vary (signs)
[ 0\. , -16.01393278, 3.72202865],
[ 0\. , 0\. , -6.13010256],
[ 0\. , 0\. , 0\. ],
[ 0\. , 0\. , 0\. ]])
更新相当于但比以下更快。
>>> a1 = np.insert(a, 2, u, 0)
>>> a1
array([[ 3., -2., -2.],
[ 6., -7., 4.],
[ 6., -9., -3.],
[ -3., 10., 1.],
[ 7., 8., -6.]])
>>> q_direct, r_direct = linalg.qr(a1)
检查我们是否有相同的结果:
>>> np.dot(q1, r1)
array([[ 3., -2., -2.],
[ 6., -7., 4.],
[ 6., -9., -3.],
[ -3., 10., 1.],
[ 7., 8., -6.]])
>>> np.allclose(np.dot(q1, r1), a1)
True
并且更新后的 Q 仍然是单位的:
>>> np.allclose(np.dot(q1.T, q1), np.eye(5))
True
scipy.linalg.rq
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.rq.html#scipy.linalg.rq
scipy.linalg.rq(a, overwrite_a=False, lwork=None, mode='full', check_finite=True)
计算矩阵的 RQ 分解。
计算分解A = R Q
,其中 Q 是酉/正交的,R 是上三角形的。
参数:
a(M, N) 类似数组
要分解的矩阵
overwrite_a布尔型,可选
是否覆盖 a 中的数据(可能会提高性能)
lwork整型,可选
工作数组大小,lwork >= a.shape[1]。如果为 None 或-1,则计算一个最佳大小。
mode,可选
决定返回哪些信息:Q 和 R 都返回(‘full’,默认),只返回 R(‘r’),或者返回经济尺寸计算的 Q 和 R(‘economic’,参见注意事项)。
check_finite布尔型,可选
是否检查输入矩阵仅包含有限数。禁用可能会带来性能提升,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃、非终止)。
返回:
R浮点数或复数的 ndarray
形状为(M, N)或(M, K),对于mode='economic'
,K = min(M, N)
。
Q浮点数或复数的 ndarray
形状为(N, N)或(K, N),对于mode='economic'
。如果mode='r'
,则不返回。
抛出:
LinAlgError
如果分解失败。
注意事项
这是 LAPACK 例程 sgerqf,dgerqf,cgerqf,zgerqf,sorgrq,dorgrq,cungrq 和 zungrq 的接口。
如果mode=economic
,则 Q 和 R 的形状为(K, N)和(M, K),而不是(N,N)和(M,N),其中K=min(M,N)
。
示例
>>> import numpy as np
>>> from scipy import linalg
>>> rng = np.random.default_rng()
>>> a = rng.standard_normal((6, 9))
>>> r, q = linalg.rq(a)
>>> np.allclose(a, r @ q)
True
>>> r.shape, q.shape
((6, 9), (9, 9))
>>> r2 = linalg.rq(a, mode='r')
>>> np.allclose(r, r2)
True
>>> r3, q3 = linalg.rq(a, mode='economic')
>>> r3.shape, q3.shape
((6, 6), (6, 9))
scipy.linalg.qz
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.qz.html#scipy.linalg.qz
scipy.linalg.qz(A, B, output='real', lwork=None, sort=None, overwrite_a=False, overwrite_b=False, check_finite=True)
用于一对矩阵的广义特征值的 QZ 分解。
一对 n 乘 n 矩阵(A,B)的 QZ 或广义舒尔分解是:
(A,B) = (Q @ AA @ Z*, Q @ BB @ Z*)
如果 BB 是具有非负对角线的上三角形状,且 AA 是上三角形状,则 AA,BB 位于广义舒尔形式中;或者对于实 QZ 分解(output='real'
)块上三角形状,具有 1x1 和 2x2 块。在这种情况下,1x1 块对应于实广义特征值,而 2x2 块通过使 BB 的对应元素具有以下形式而‘标准化’:
[ a 0 ]
[ 0 b ]
并且 AA 和 BB 的对应的 2x2 块将具有一对复共轭的广义特征值。如果(output='complex'
)或 A 和 B 是复矩阵,则 Z’表示 Z 的共轭转置。Q 和 Z 是酉矩阵。
参数:
A(N, N) array_like
用于分解的二维数组
B(N, N) array_like
用于分解的二维数组
output,可选
构建实数或复数矩阵的 QZ 分解。默认为‘real’。
lworkint,可选
工作数组大小。如果为 None 或-1,则会自动计算。
sort,可选
注意:此输入目前已禁用。请使用 ordqz 代替。
指定是否应对上层特征值进行排序。可以传递一个可调用函数,给定一个特征值,返回一个布尔值,表示是否应将特征值排序到左上角(True)。对于实矩阵对,排序函数接受三个实参数(alphar, alphai, beta)。特征值 x = (alphar + alphai*1j)/beta
。对于复矩阵对或者 output=’complex’,排序函数接受两个复参数(alpha, beta)。特征值 x = (alpha/beta)
。也可以使用字符串参数:
- ‘lhp’ 左平面(x.real < 0.0)
- ‘rhp’ 右平面(x.real > 0.0)
- ‘iuc’ 单位圆内部(x*x.conjugate() < 1.0)
- ‘ouc’ 单位圆外部(x*x.conjugate() > 1.0)
默认为 None(不排序)。
overwrite_abool,可选
是否覆盖 a 中的数据(可能提高性能)
overwrite_bbool,可选
是否覆盖 b 中的数据(可能提高性能)
check_finitebool,可选
如果为 true,则检查A和B的元素是否为有限数。如果为 false,则不进行检查并将矩阵传递给底层算法。
返回:
AA(N, N) ndarray
一般化的 A 的舒尔形式。
BB(N, N) ndarray
一般化的 B 的舒尔形式。
Q(N, N) ndarray
左舒尔向量。
Z(N, N) ndarray
右舒尔向量。
参见
ordqz
注释
Q 相对于 Matlab 中等效函数是转置的。
新版本 0.11.0 中的新增内容。
示例
>>> import numpy as np
>>> from scipy.linalg import qz
>>> A = np.array([[1, 2, -1], [5, 5, 5], [2, 4, -8]])
>>> B = np.array([[1, 1, -3], [3, 1, -1], [5, 6, -2]])
计算分解。QZ 分解不唯一,因此根据所使用的基础库不同,以下输出中系数的符号可能会有所不同。
>>> AA, BB, Q, Z = qz(A, B)
>>> AA
array([[-1.36949157, -4.05459025, 7.44389431],
[ 0\. , 7.65653432, 5.13476017],
[ 0\. , -0.65978437, 2.4186015 ]]) # may vary
>>> BB
array([[ 1.71890633, -1.64723705, -0.72696385],
[ 0\. , 8.6965692 , -0\. ],
[ 0\. , 0\. , 2.27446233]]) # may vary
>>> Q
array([[-0.37048362, 0.1903278 , 0.90912992],
[-0.90073232, 0.16534124, -0.40167593],
[ 0.22676676, 0.96769706, -0.11017818]]) # may vary
>>> Z
array([[-0.67660785, 0.63528924, -0.37230283],
[ 0.70243299, 0.70853819, -0.06753907],
[ 0.22088393, -0.30721526, -0.92565062]]) # may vary
验证 QZ 分解。对于实数输出,在以下表达式中我们只需要Z
的转置。
>>> Q @ AA @ Z.T # Should be A
array([[ 1., 2., -1.],
[ 5., 5., 5.],
[ 2., 4., -8.]])
>>> Q @ BB @ Z.T # Should be B
array([[ 1., 1., -3.],
[ 3., 1., -1.],
[ 5., 6., -2.]])
重复分解,但使用output='complex'
。
>>> AA, BB, Q, Z = qz(A, B, output='complex')
为了输出简洁,我们使用np.set_printoptions()
来将 NumPy 数组的输出精度设置为 3,并将微小值显示为 0。
>>> np.set_printoptions(precision=3, suppress=True)
>>> AA
array([[-1.369+0.j , 2.248+4.237j, 4.861-5.022j],
[ 0\. +0.j , 7.037+2.922j, 0.794+4.932j],
[ 0\. +0.j , 0\. +0.j , 2.655-1.103j]]) # may vary
>>> BB
array([[ 1.719+0.j , -1.115+1.j , -0.763-0.646j],
[ 0\. +0.j , 7.24 +0.j , -3.144+3.322j],
[ 0\. +0.j , 0\. +0.j , 2.732+0.j ]]) # may vary
>>> Q
array([[ 0.326+0.175j, -0.273-0.029j, -0.886-0.052j],
[ 0.794+0.426j, -0.093+0.134j, 0.402-0.02j ],
[-0.2 -0.107j, -0.816+0.482j, 0.151-0.167j]]) # may vary
>>> Z
array([[ 0.596+0.32j , -0.31 +0.414j, 0.393-0.347j],
[-0.619-0.332j, -0.479+0.314j, 0.154-0.393j],
[-0.195-0.104j, 0.576+0.27j , 0.715+0.187j]]) # may vary
对于复数数组,在以下表达式中我们必须使用Z.conj().T
来验证分解。
>>> Q @ AA @ Z.conj().T # Should be A
array([[ 1.-0.j, 2.-0.j, -1.-0.j],
[ 5.+0.j, 5.+0.j, 5.-0.j],
[ 2.+0.j, 4.+0.j, -8.+0.j]])
>>> Q @ BB @ Z.conj().T # Should be B
array([[ 1.+0.j, 1.+0.j, -3.+0.j],
[ 3.-0.j, 1.-0.j, -1.+0.j],
[ 5.+0.j, 6.+0.j, -2.+0.j]])
scipy.linalg.ordqz
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.ordqz.html#scipy.linalg.ordqz
scipy.linalg.ordqz(A, B, sort='lhp', output='real', overwrite_a=False, overwrite_b=False, check_finite=True)
用于一对重新排序矩阵的 QZ 分解。
参数:
A(N, N) array_like
2-D 数组进行分解
B(N, N) array_like
2-D 数组进行分解
sort,可选
指定是否应对上部特征值进行排序。可以传递一个可调用函数,给定一个有序对 (alpha, beta)
表示特征值 x = (alpha/beta)
,返回一个布尔值,表示是否应将特征值排序到左上角(True)。对于实矩阵对,beta 是实数,而 alpha 可以是复数;对于复杂矩阵对,alpha 和 beta 都可以是复数。该可调用函数必须能够接受一个 NumPy 数组。另外,也可以使用字符串参数:
- ‘lhp’ 左半平面(x.real < 0.0)
- ‘rhp’ 右半平面(x.real > 0.0)
- ‘iuc’ 单位圆内(x*x.conjugate() < 1.0)
- ‘ouc’ 单位圆外(x*x.conjugate() > 1.0)
使用预定义的排序函数,无穷特征值(即 alpha != 0
且 beta = 0
)被认为既不位于左半平面也不位于右半平面,但被认为位于单位圆外。对于特征值 (alpha, beta) = (0, 0)
,预定义的排序函数都返回 False。
outputstr {'real','complex'},可选
构造实数或复数 QZ 分解的真实矩阵。默认为 'real'。
overwrite_abool, optional
如果为真,则覆盖 A 的内容。
overwrite_bbool, optional
如果为真,则覆盖 B 的内容。
check_finitebool, optional
如果为真,则检查 A 和 B 的元素是否为有限数。如果为假,则不进行检查并将矩阵传递给底层算法。
返回:
AA(N, N) ndarray
A 的广义舒尔形式。
BB(N, N) ndarray
B 的广义舒尔形式。
alpha(N,) ndarray
alpha = alphar + alphai * 1j。请参阅备注。
beta(N,) ndarray
请参阅备注。
Q(N, N) ndarray
左舒尔向量。
Z(N, N) ndarray
右舒尔向量。
另请参阅
备注
在退出时,(ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N
,将是广义特征值。 ALPHAR(j) + ALPHAI(j)*i
和 BETA(j),j=1,...,N
是复杂舒尔形式(S,T)的对角线,如果实广义舒尔形式(A,B)的 2×2 对角块进一步通过复杂酉变换化为三角形式,则结果如此。如果 ALPHAI(j) 为零,则第 j 个特征值为实数;如果为正,则第 j
个和 (j+1)
个特征值为复共轭对,其中 ALPHAI(j+1)
为负数。
自版本 0.17.0 起新增。
示例
>>> import numpy as np
>>> from scipy.linalg import ordqz
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> B = np.array([[0, 6, 0, 0], [5, 0, 2, 1], [5, 2, 6, 6], [4, 7, 7, 7]])
>>> AA, BB, alpha, beta, Q, Z = ordqz(A, B, sort='lhp')
由于我们已对左半平面特征值进行了排序,负值首先出现
>>> (alpha/beta).real < 0
array([ True, True, False, False], dtype=bool)
scipy.linalg.schur
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.schur.html#scipy.linalg.schur
scipy.linalg.schur(a, output='real', lwork=None, overwrite_a=False, sort=None, check_finite=True)
计算矩阵的舒尔分解。
舒尔分解为:
A = Z T Z^H
其中 Z 为酉,T 为上三角,或对于实舒尔分解(output=’real’),准上三角。在准三角形式中,描述复值特征值对的 2x2 块可能从对角线突出。
参数:
a(M, M) array_like
矩阵分解
output,可选
构造实数或复数舒尔分解(对于实矩阵)。
lwork整数,可选
工作数组大小。如果为 None 或-1,则会自动计算。
overwrite_a布尔值,可选
是否覆盖数据在(可能提高性能)。
sort,可选
指定是否应对上特征值进行排序。可以传递一个可调用对象,给定一个特征值,返回一个布尔值,表示是否应将该特征值排序到左上角(True)。另外,也可以使用字符串参数:
'lhp' Left-hand plane (x.real < 0.0)
'rhp' Right-hand plane (x.real > 0.0)
'iuc' Inside the unit circle (x*x.conjugate() <= 1.0)
'ouc' Outside the unit circle (x*x.conjugate() > 1.0)
默认为 None(不排序)。
check_finite布尔值,可选
是否检查输入矩阵仅包含有限数。禁用可能会带来性能提升,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃、非终止)。
返回:
T(M, M) 数组
A 的舒尔形式。对于实数舒尔分解,它是实数值的。
Z(M, M) 数组
一个 A 的酉舒尔变换矩阵。对于实数舒尔分解,它是实数值的。
sdim整数
如果只有在请求排序时,第三个返回值才会包含满足排序条件的特征值数量。
引发:
线性代数错误
三种条件下引发的错误:
-
算法由于 QR 算法未能计算所有特征值而失败。
-
如果请求特征值排序,由于未能分离特征值而导致无法重新排序特征值,通常是因为条件不佳。
-
如果请求特征值排序,由于舍入误差导致主特征值不再满足排序条件。
另见
rsf2csf
将实数舒尔形式转换为复数舒尔形式
示例
>>> import numpy as np
>>> from scipy.linalg import schur, eigvals
>>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
>>> T, Z = schur(A)
>>> T
array([[ 2.65896708, 1.42440458, -1.92933439],
[ 0\. , -0.32948354, -0.49063704],
[ 0\. , 1.31178921, -0.32948354]])
>>> Z
array([[0.72711591, -0.60156188, 0.33079564],
[0.52839428, 0.79801892, 0.28976765],
[0.43829436, 0.03590414, -0.89811411]])
>>> T2, Z2 = schur(A, output='complex')
>>> T2
array([[ 2.65896708, -1.22839825+1.32378589j, 0.42590089+1.51937378j], # may vary
[ 0\. , -0.32948354+0.80225456j, -0.59877807+0.56192146j],
[ 0\. , 0\. , -0.32948354-0.80225456j]])
>>> eigvals(T2)
array([2.65896708, -0.32948354+0.80225456j, -0.32948354-0.80225456j])
一个任意的自定义特征值排序条件,具有正虚部,仅由一个特征值满足
>>> T3, Z3, sdim = schur(A, output='complex', sort=lambda x: x.imag > 0)
>>> sdim
1
scipy.linalg.rsf2csf
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.rsf2csf.html#scipy.linalg.rsf2csf
scipy.linalg.rsf2csf(T, Z, check_finite=True)
将实 Schur 形式转换为复 Schur 形式。
将准对角的实值 Schur 形式转换为上三角形的复值 Schur 形式。
参数:
T(M, M) array_like
原始数组的实 Schur 形式
Z(M, M) array_like
Schur 变换矩阵
check_finitebool,可选
是否检查输入数组仅包含有限数。禁用可能会提高性能,但如果输入确实包含无穷大或 NaN,则可能会导致问题(崩溃、非终止)。
返回值:
T(M, M) ndarray
原始数组的复 Schur 形式
Z(M, M) ndarray
对应于复形式的 Schur 变换矩阵
参见
schur
数组的 Schur 分解
示例
>>> import numpy as np
>>> from scipy.linalg import schur, rsf2csf
>>> A = np.array([[0, 2, 2], [0, 1, 2], [1, 0, 1]])
>>> T, Z = schur(A)
>>> T
array([[ 2.65896708, 1.42440458, -1.92933439],
[ 0\. , -0.32948354, -0.49063704],
[ 0\. , 1.31178921, -0.32948354]])
>>> Z
array([[0.72711591, -0.60156188, 0.33079564],
[0.52839428, 0.79801892, 0.28976765],
[0.43829436, 0.03590414, -0.89811411]])
>>> T2 , Z2 = rsf2csf(T, Z)
>>> T2
array([[2.65896708+0.j, -1.64592781+0.743164187j, -1.21516887+1.00660462j],
[0.+0.j , -0.32948354+8.02254558e-01j, -0.82115218-2.77555756e-17j],
[0.+0.j , 0.+0.j, -0.32948354-0.802254558j]])
>>> Z2
array([[0.72711591+0.j, 0.28220393-0.31385693j, 0.51319638-0.17258824j],
[0.52839428+0.j, 0.24720268+0.41635578j, -0.68079517-0.15118243j],
[0.43829436+0.j, -0.76618703+0.01873251j, -0.03063006+0.46857912j]])
scipy.linalg.hessenberg
scipy.linalg.hessenberg(a, calc_q=False, overwrite_a=False, check_finite=True)
计算矩阵的 Hessenberg 形式。
Hessenberg 分解为:
A = Q H Q^H
其中 Q 是单位 ary/正交的,H 除了第一个次对角线以下的元素外都为零。
参数:
a(M, M) array_like
要转换为 Hessenberg 形式的矩阵。
calc_qbool, 可选
是否计算变换矩阵。默认为 False。
overwrite_abool, 可选
是否覆盖 a;可能提高性能。默认为 False。
check_finitebool, 可选
是否检查输入矩阵仅包含有限数字。禁用可能提高性能,但如果输入包含无穷大或 NaN,可能导致问题(崩溃、非终止)。
返回:
H(M, M) ndarray
a 的 Hessenberg 形式。
Q(M, M) ndarray
单位 ary/正交相似变换矩阵 A = Q H Q^H
。仅在 calc_q=True
时返回。
示例
>>> import numpy as np
>>> from scipy.linalg import hessenberg
>>> A = np.array([[2, 5, 8, 7], [5, 2, 2, 8], [7, 5, 6, 6], [5, 4, 4, 8]])
>>> H, Q = hessenberg(A, calc_q=True)
>>> H
array([[ 2\. , -11.65843866, 1.42005301, 0.25349066],
[ -9.94987437, 14.53535354, -5.31022304, 2.43081618],
[ 0\. , -1.83299243, 0.38969961, -0.51527034],
[ 0\. , 0\. , -3.83189513, 1.07494686]])
>>> np.allclose(Q @ H @ Q.conj().T - A, np.zeros((4, 4)))
True
scipy.linalg.cdf2rdf
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.cdf2rdf.html#scipy.linalg.cdf2rdf
scipy.linalg.cdf2rdf(w, v)
将复数特征值w
和特征向量v
转换为实块对角形式的实特征值wr
及相关的实特征向量vr
,使得:
vr @ wr = X @ vr
保持不变,其中X
是w
和v
是其特征值和特征向量的原始数组。
1.1.0 版本新增。
参数:
w(…, M) array_like
复数或实特征值,数组或数组堆栈
如果交错排列共轭对,将会产生错误结果。因此,[1+1j, 1, 1-1j]
将给出正确结果,但[1+1j, 2+1j, 1-1j, 2-1j]
则不会。
v(…, M, M) array_like
复数或实特征向量,方阵或方阵堆栈。
返回:
wr(…, M, M) ndarray
特征值的实对角块形式
vr(…, M, M) ndarray
与wr
相关的实特征向量
参见
eig
函数
对于非对称数组的特征值和右特征向量
rsf2csf
函数
将实舒尔形式转换为复舒尔形式
注释
w
、v
必须是某些实矩阵X
的特征结构,例如通过w, v = scipy.linalg.eig(X)
或w, v = numpy.linalg.eig(X)
获得,其中X
也可以表示为堆叠的数组。
1.1.0 版本新增。
示例
>>> import numpy as np
>>> X = np.array([[1, 2, 3], [0, 4, 5], [0, -5, 4]])
>>> X
array([[ 1, 2, 3],
[ 0, 4, 5],
[ 0, -5, 4]])
>>> from scipy import linalg
>>> w, v = linalg.eig(X)
>>> w
array([ 1.+0.j, 4.+5.j, 4.-5.j])
>>> v
array([[ 1.00000+0.j , -0.01906-0.40016j, -0.01906+0.40016j],
[ 0.00000+0.j , 0.00000-0.64788j, 0.00000+0.64788j],
[ 0.00000+0.j , 0.64788+0.j , 0.64788-0.j ]])
>>> wr, vr = linalg.cdf2rdf(w, v)
>>> wr
array([[ 1., 0., 0.],
[ 0., 4., 5.],
[ 0., -5., 4.]])
>>> vr
array([[ 1\. , 0.40016, -0.01906],
[ 0\. , 0.64788, 0\. ],
[ 0\. , 0\. , 0.64788]])
>>> vr @ wr
array([[ 1\. , 1.69593, 1.9246 ],
[ 0\. , 2.59153, 3.23942],
[ 0\. , -3.23942, 2.59153]])
>>> X @ vr
array([[ 1\. , 1.69593, 1.9246 ],
[ 0\. , 2.59153, 3.23942],
[ 0\. , -3.23942, 2.59153]])
scipy.linalg.cossin
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.cossin.html#scipy.linalg.cossin
scipy.linalg.cossin(X, p=None, q=None, separate=False, swap_sign=False, compute_u=True, compute_vh=True)
计算正交/酉矩阵的余弦-正弦(CS)分解。
X 是一个(m, m)
正交/酉矩阵,分块如下,其中左上角块的形状为(p, q)
:
┌ ┐
│ I 0 0 │ 0 0 0 │
┌ ┐ ┌ ┐│ 0 C 0 │ 0 -S 0 │┌ ┐*
│ X11 │ X12 │ │ U1 │ ││ 0 0 0 │ 0 0 -I ││ V1 │ │
│ ────┼──── │ = │────┼────││─────────┼─────────││────┼────│
│ X21 │ X22 │ │ │ U2 ││ 0 0 0 │ I 0 0 ││ │ V2 │
└ ┘ └ ┘│ 0 S 0 │ 0 C 0 │└ ┘
│ 0 0 I │ 0 0 0 │
└ ┘
U1
, U2
, V1
, V2
是维度分别为(p,p)
,(m-p,m-p)
,(q,q)
和 (m-q,m-q)
的方正交/酉矩阵,C
和 S
是满足 C² + S² = I
的(r,r)
非负对角矩阵,其中 r = min(p, m-p, q, m-q)
。
此外,单位矩阵的秩分别为min(p, q) - r
,min(p, m - q) - r
,min(m - p, q) - r
和min(m - p, m - q) - r
。
X 可以通过其自身和块规格 p
, q
或其子块的可迭代对象提供。参见下面的示例。
参数:
X类数组,可迭代对象
要分解的复数酉或实正交矩阵,或子块 X11
, X12
, X21
, X22
的可迭代对象,当省略 p
, q
时。
p整数,可选
左上角块 X11
的行数,仅在给定 X 作为数组时使用。
q整数,可选
左上角块 X11
的列数,仅在给定 X 作为数组时使用。
separate布尔值,可选
如果为True
,则返回低级组件而不是矩阵因子,即 (u1,u2)
, theta
, (v1h,v2h)
而不是 u
, cs
, vh
。
swap_sign布尔值,可选
如果为True
,则-S
,-I
块将位于左下角,否则(默认情况下)它们将位于右上角。
compute_u布尔值,可选
如果为False
,u
将不会被计算,并返回一个空数组。
compute_vh布尔值,可选
如果为False
,vh
将不会被计算,并返回一个空数组。
返回:
u数组
当compute_u=True
时,包含由块对角线正交/酉矩阵组成的块 U1
(p
x p
) 和 U2
(m-p
x m-p
)。如果separate=True
,则包含元组(U1, U2)
。
cs数组
具有上述结构的余弦-正弦因子。
如果separate=True
,则包含角度以弧度表示的theta
数组。
vh数组
当pycompute_vh=True`, contains the block diagonal orthogonal/unitary matrix consisting of the blocks ``V1H
(q
x q
) 和 V2H
(m-q
x m-q
) 正交/酉矩阵。如果separate=True
,则包含元组(V1H, V2H)
。
参考文献
[1]
Brian D. Sutton. 计算完整的 CS 分解。Numer. Algorithms, 50(1):33-65, 2009.
示例
>>> import numpy as np
>>> from scipy.linalg import cossin
>>> from scipy.stats import unitary_group
>>> x = unitary_group.rvs(4)
>>> u, cs, vdh = cossin(x, p=2, q=2)
>>> np.allclose(x, u @ cs @ vdh)
True
也可以通过子块输入,无需 p
和 q
。还让我们跳过 u
的计算。
>>> ue, cs, vdh = cossin((x[:2, :2], x[:2, 2:], x[2:, :2], x[2:, 2:]),
... compute_u=False)
>>> print(ue)
[]
>>> np.allclose(x, u @ cs @ vdh)
True
scipy.linalg.expm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.expm.html#scipy.linalg.expm
scipy.linalg.expm(A)
计算数组的矩阵指数。
参数:
Andarray
输入的最后两个维度是方形的(..., n, n)
。
返回:
eAndarray
结果矩阵指数与A
的形状相同
注意
实现了在[1]中给出的算法,这实质上是一种带有基于数组数据决定的可变阶数的 Pade 逼近。
对于大小为n
的输入,在最坏情况下,内存使用量是8*(n**2)
的数量级。如果输入数据不是单精度和双精度的实数和复数数据类型,则将其复制到一个新的数组。
对于n >= 400
的情况,精确的 1-范数计算成本与 1-范数估计持平,并且从那一点开始,使用[2]中给出的估计方案来决定逼近阶数。
参考文献
[1]
Awad H. Al-Mohy 和 Nicholas J. Higham(2009),"矩阵指数的新缩放和平方算法",SIAM J. Matrix Anal. Appl. 31(3):970-989,DOI:10.1137/09074721X
[2]
Nicholas J. Higham 和 Francoise Tisseur(2000),"用于矩阵 1-范数估计的块算法,及其在 1-范数伪谱中的应用",SIAM J. Matrix Anal. Appl. 21(4):1185-1201,DOI:10.1137/S0895479899356080
示例
>>> import numpy as np
>>> from scipy.linalg import expm, sinm, cosm
公式 exp(0) = 1 的矩阵版本:
>>> expm(np.zeros((3, 2, 2)))
array([[[1., 0.],
[0., 1.]],
[[1., 0.],
[0., 1.]],
[[1., 0.],
[0., 1.]]])
欧拉恒等式(exp(itheta) = cos(theta) + isin(theta))应用于矩阵:
>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
scipy.linalg.logm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.logm.html#scipy.linalg.logm
scipy.linalg.logm(A, disp=True)
计算矩阵对数。
矩阵对数是 expm 的逆:expm(logm(A)) == A
参数:
A(N, N) 类似数组
要评估其对数的矩阵
disp 布尔值,可选项
如果估计的结果误差较大,则打印警告,而不是返回估计的误差。 (默认:True)
返回:
logm(N, N) ndarray
A 的矩阵对数
errest 浮点数
(如果 disp == False)
估计误差的 1-范数,||err||_1 / ||A||_1
参考文献
[1]
Awad H. Al-Mohy 和 Nicholas J. Higham (2012) “矩阵对数的改进逆缩放和平方算法。”《SIAM 科学计算杂志》,34 (4). C152-C169. ISSN 1095-7197
[2]
Nicholas J. Higham (2008) “矩阵函数:理论与计算” ISBN 978-0-898716-46-7
[3]
Nicholas J. Higham 和 Lijing lin (2011) “矩阵的分数幂的 Schur-Pade 算法。”《SIAM 矩阵分析和应用杂志》,32 (3). pp. 1056-1078. ISSN 0895-4798
示例
>>> import numpy as np
>>> from scipy.linalg import logm, expm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> b = logm(a)
>>> b
array([[-1.02571087, 2.05142174],
[ 0.68380725, 1.02571087]])
>>> expm(b) # Verify expm(logm(a)) returns a
array([[ 1., 3.],
[ 1., 4.]])
scipy.linalg.cosm
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.cosm.html#scipy.linalg.cosm
scipy.linalg.cosm(A)
计算矩阵余弦。
这个 这个例程使用 expm 来计算矩阵指数。
参数:
A(N, N) array_like
输入数组
返回值:
cosm(N, N) ndarray
A 的矩阵余弦
例子
>>> import numpy as np
>>> from scipy.linalg import expm, sinm, cosm
应用于矩阵的欧拉恒等式(exp(itheta) = cos(theta) + isin(theta)):
>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
scipy.linalg.sinm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.sinm.html#scipy.linalg.sinm
scipy.linalg.sinm(A)
计算矩阵正弦。
此例程使用expm
来计算矩阵的指数。
参数:
A(N, N) array_like
输入数组。
返回:
sinm(N, N) ndarray
A 的矩阵正弦
示例
>>> import numpy as np
>>> from scipy.linalg import expm, sinm, cosm
Euler’s identity (exp(i*theta) = cos(theta) + i*sin(theta)
) 应用于矩阵:
>>> a = np.array([[1.0, 2.0], [-1.0, 3.0]])
>>> expm(1j*a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
>>> cosm(a) + 1j*sinm(a)
array([[ 0.42645930+1.89217551j, -2.13721484-0.97811252j],
[ 1.06860742+0.48905626j, -1.71075555+0.91406299j]])
scipy.linalg.tanm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.tanm.html#scipy.linalg.tanm
scipy.linalg.tanm(A)
计算矩阵的切线。
此例程使用 expm 计算矩阵指数。
参数:
A(N, N) 数组样式
输入数组。
返回:
tanm(N, N) ndarray
A 的矩阵切线
示例
>>> import numpy as np
>>> from scipy.linalg import tanm, sinm, cosm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> t = tanm(a)
>>> t
array([[ -2.00876993, -8.41880636],
[ -2.80626879, -10.42757629]])
验证 tanm(a) = sinm(a).dot(inv(cosm(a)))
>>> s = sinm(a)
>>> c = cosm(a)
>>> s.dot(np.linalg.inv(c))
array([[ -2.00876993, -8.41880636],
[ -2.80626879, -10.42757629]])
scipy.linalg.coshm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.coshm.html#scipy.linalg.coshm
scipy.linalg.coshm(A)
计算双曲矩阵余弦。
此例程使用 expm 来计算矩阵指数。
参数:
A(N, N) 数组类似物
输入数组。
返回:
coshm(N, N) 数组形状
A 的双曲矩阵余弦
示例
>>> import numpy as np
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> c = coshm(a)
>>> c
array([[ 11.24592233, 38.76236492],
[ 12.92078831, 50.00828725]])
验证 tanhm(a) = sinhm(a).dot(inv(coshm(a)))
>>> t = tanhm(a)
>>> s = sinhm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[ 2.72004641e-15, 4.55191440e-15],
[ 0.00000000e+00, -5.55111512e-16]])
scipy.linalg.sinhm
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.sinhm.html#scipy.linalg.sinhm
scipy.linalg.sinhm(A)
计算双曲矩阵正弦。
此例程使用 expm 计算矩阵指数。
参数:
A(N, N) array_like
输入数组。
返回:
sinhm(N, N) ndarray
A 的双曲矩阵正弦
示例
>>> import numpy as np
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> s = sinhm(a)
>>> s
array([[ 10.57300653, 39.28826594],
[ 13.09608865, 49.86127247]])
验证 tanhm(a) = sinhm(a).dot(inv(coshm(a)))
>>> t = tanhm(a)
>>> c = coshm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[ 2.72004641e-15, 4.55191440e-15],
[ 0.00000000e+00, -5.55111512e-16]])
scipy.linalg.tanhm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.tanhm.html#scipy.linalg.tanhm
scipy.linalg.tanhm(A)
计算双曲矩阵切线。
该例程使用 expm 计算矩阵指数。
参数:
A(N, N) array_like
输入数组
返回:
tanhm(N, N) ndarray
A 的双曲矩阵切线
示例
>>> import numpy as np
>>> from scipy.linalg import tanhm, sinhm, coshm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> t = tanhm(a)
>>> t
array([[ 0.3428582 , 0.51987926],
[ 0.17329309, 0.86273746]])
验证 tanhm(a) = sinhm(a).dot(inv(coshm(a)))
>>> s = sinhm(a)
>>> c = coshm(a)
>>> t - s.dot(np.linalg.inv(c))
array([[ 2.72004641e-15, 4.55191440e-15],
[ 0.00000000e+00, -5.55111512e-16]])
scipy.linalg.signm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.signm.html#scipy.linalg.signm
scipy.linalg.signm(A, disp=True)
矩阵签名函数。
标量 sign(x)对矩阵的扩展。
参数:
A(N, N) 数组型
评估签名函数的矩阵
disp布尔值,可选
打印警告,如果结果中估计的误差较大而不是返回估计的错误。(默认:True)
返回:
signm(N, N) 数组型
签名函数在A处的值
errest浮点数
(如果 disp == False)
1-范数的估计误差,||err||_1 / ||A||_1
示例
>>> from scipy.linalg import signm, eigvals
>>> a = [[1,2,3], [1,2,1], [1,1,1]]
>>> eigvals(a)
array([ 4.12488542+0.j, -0.76155718+0.j, 0.63667176+0.j])
>>> eigvals(signm(a))
array([-1.+0.j, 1.+0.j, 1.+0.j])
scipy.linalg.sqrtm
原文:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.sqrtm.html#scipy.linalg.sqrtm
scipy.linalg.sqrtm(A, disp=True, blocksize=64)
矩阵平方根。
参数:
A(N, N) array_like
要评估其平方根的矩阵
disp布尔值,可选
如果结果中的误差估计较大,则打印警告,而不是返回估计的误差。(默认:True)
blocksize整数,可选
如果块大小与输入数组的大小不同,则使用块算法。(默认:64)
返回:
sqrtm(N, N) ndarray
A处的 sqrt 函数值。数据类型为 float 或 complex。精度(数据大小)基于输入A的精度。当数据类型为 float 时,精度与A相同。当数据类型为 complex 时,精度是A的两倍。每种数据类型的精度可能会被各自的精度范围所限制。
errest浮点数
(如果 disp == False)
估计误差的 Frobenius 范数,||err||_F / ||A||_F
参考文献
[1]
Edvin Deadman, Nicholas J. Higham, Rui Ralha (2013) “Blocked Schur Algorithms for Computing the Matrix Square Root, Lecture Notes in Computer Science, 7782. pp. 171-182.
示例
>>> import numpy as np
>>> from scipy.linalg import sqrtm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> r = sqrtm(a)
>>> r
array([[ 0.75592895, 1.13389342],
[ 0.37796447, 1.88982237]])
>>> r.dot(r)
array([[ 1., 3.],
[ 1., 4.]])
scipy.linalg.funm
原文链接:
docs.scipy.org/doc/scipy-1.12.0/reference/generated/scipy.linalg.funm.html#scipy.linalg.funm
scipy.linalg.funm(A, func, disp=True)
评估由可调用对象指定的矩阵函数。
返回函数 f
在 A 处的矩阵值。函数 f
是将标量函数 func 推广到矩阵的扩展。
参数:
A(N, N) array_like
用于评估函数的矩阵
funccallable
评估标量函数 f 的可调用对象。必须是矢量化的(例如,使用 vectorize)。
dispbool, 可选
如果结果中的误差估计较大,则打印警告而不是返回估计的误差。(默认:True)
返回:
funm(N, N) ndarray
在 A 处评估的由 func 指定的矩阵函数的值
errestfloat
(如果 disp == False)
估计误差的 1-范数,||err||_1 / ||A||_1
注意
该函数实现了基于舒尔分解的一般算法(算法 9.1.1 在 [1] 中)。
如果已知输入矩阵可对角化,则依赖于特征分解可能更快。例如,如果你的矩阵是埃尔米特矩阵,你可以执行
>>> from scipy.linalg import eigh
>>> def funm_herm(a, func, check_finite=False):
... w, v = eigh(a, check_finite=check_finite)
... ## if you further know that your matrix is positive semidefinite,
... ## you can optionally guard against precision errors by doing
... # w = np.maximum(w, 0)
... w = func(w)
... return (v * w).dot(v.conj().T)
参考文献
[1]
Gene H. Golub, Charles F. van Loan, 《Matrix Computations》第四版。
示例
>>> import numpy as np
>>> from scipy.linalg import funm
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> funm(a, lambda x: x*x)
array([[ 4., 15.],
[ 5., 19.]])
>>> a.dot(a)
array([[ 4., 15.],
[ 5., 19.]])
scipy.linalg.expm_frechet
scipy.linalg.expm_frechet(A, E, method=None, compute_expm=True, check_finite=True)
A 的矩阵指数在 E 方向上的 Frechet 导数。
参数:
A(N, N) 类似数组
矩阵的矩阵指数。
E(N, N) 类似数组
用于计算 Frechet 导数的矩阵方向。
method字符串,可选
算法的选择。应该是以下之一:
-
SPS(默认)
-
blockEnlarge
compute_expm布尔值,可选
是否同时计算expm_A和expm_frechet_AE。默认为 True。
check_finite布尔值,可选
是否检查输入矩阵是否仅包含有限数值。禁用可能会提高性能,但如果输入包含无穷大或 NaN,则可能导致问题(崩溃、非终止)。
返回:
expm_A数组
A 的矩阵指数。
expm_frechet_AE数组
A 的矩阵指数在 E 方向上的 Frechet 导数。
对于compute_expm = False
,只返回expm_frechet_AE。
另请参见:
计算矩阵的指数。
注:
本节描述了可以通过method参数选择的可用实现。默认方法是SPS。
方法blockEnlarge是一个朴素的算法。
方法SPS是 Scaling-Pade-Squaring [1]。这是一个复杂的实现,其执行时间只需朴素实现的 3/8。渐近性质相同。
从版本 0.13.0 开始。
参考资料:
Awad H. Al-Mohy 和 Nicholas J. Higham(2009 年)计算矩阵指数的 Frechet 导数,及其在条件数估计中的应用。SIAM Journal On Matrix Analysis and Applications.,30(4)。pp. 1639-1657。ISSN 1095-7162
示例:
>>> import numpy as np
>>> from scipy import linalg
>>> rng = np.random.default_rng()
>>> A = rng.standard_normal((3, 3))
>>> E = rng.standard_normal((3, 3))
>>> expm_A, expm_frechet_AE = linalg.expm_frechet(A, E)
>>> expm_A.shape, expm_frechet_AE.shape
((3, 3), (3, 3))
创建一个包含[[A, E], [0, A]]的 6x6 矩阵:
>>> M = np.zeros((6, 6))
>>> M[:3, :3] = A
>>> M[:3, 3:] = E
>>> M[3:, 3:] = A
>>> expm_M = linalg.expm(M)
>>> np.allclose(expm_A, expm_M[:3, :3])
True
>>> np.allclose(expm_frechet_AE, expm_M[:3, 3:])
True
scipy.linalg.expm_cond
scipy.linalg.expm_cond(A, check_finite=True)
矩阵指数在 Frobenius 范数中的相对条件数。
参数:
A2 维数组类型
形状为(N, N)的方形输入矩阵。
check_finite布尔型,可选
是否检查输入矩阵只包含有限数字。禁用可能会提高性能,但如果输入包含无穷大或 NaN,可能会导致问题(崩溃,非终止)。
返回:
kappa浮点型
矩阵指数在 Frobenius 范数中的相对条件数。
另请参见
计算矩阵的指数。
计算矩阵指数的 Frechet 导数。
注意事项
已发布 1 范数中条件数的更快估计,但尚未在 SciPy 中实现。
自版本 0.14.0 开始新增。
示例
>>> import numpy as np
>>> from scipy.linalg import expm_cond
>>> A = np.array([[-0.3, 0.2, 0.6], [0.6, 0.3, -0.1], [-0.7, 1.2, 0.9]])
>>> k = expm_cond(A)
>>> k
1.7787805864469866
scipy.linalg.fractional_matrix_power
scipy.linalg.fractional_matrix_power(A, t)
计算矩阵的分数幂。
按照[1]中第六部分的讨论进行。
参数:
- A(N, N) 数组类型
评估其分数幂的矩阵。
- tfloat
分数幂。
返回:
- X(N, N) 数组类型
矩阵的分数幂。
参考文献
[1]
Nicholas J. Higham 和 Lijing lin(2011 年)“矩阵的分数幂的舒尔-帕德算法。”《SIAM 矩阵分析和应用杂志》,32(3)。第 1056-1078 页。ISSN 0895-4798
示例
>>> import numpy as np
>>> from scipy.linalg import fractional_matrix_power
>>> a = np.array([[1.0, 3.0], [1.0, 4.0]])
>>> b = fractional_matrix_power(a, 0.5)
>>> b
array([[ 0.75592895, 1.13389342],
[ 0.37796447, 1.88982237]])
>>> np.dot(b, b) # Verify square root
array([[ 1., 3.],
[ 1., 4.]])
scipy.linalg.solve_sylvester
scipy.linalg.solve_sylvester(a, b, q)
计算 Sylvester 方程 (AX + XB = Q) 的解(X)。
参数:
a(M, M) 数组
Sylvester 方程的首部矩阵
b(N, N) 数组
Sylvester 方程的尾部矩阵
q(M, N) 数组
右手边
返回:
x(M, N) 数组
Sylvester 方程的解。
引发:
LinAlgError
如果找不到解决方案
注意事项
通过巴特尔斯-斯图尔特算法计算 Sylvester 矩阵方程的解。首先对 A 和 B 矩阵进行 Schur 分解。然后利用得到的矩阵构造一个替代的 Sylvester 方程 (RY + YS^T = F
),其中 R 和 S 矩阵处于准三角形形式(或当 R、S 或 F 是复数时,为三角形形式)。简化的方程然后直接使用 LAPACK 中的 *TRSYL
解决。
自版本 0.11.0 起新增
示例
给定 a, b 和 q 求解 x:
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[-3, -2, 0], [-1, -1, 3], [3, -5, -1]])
>>> b = np.array([[1]])
>>> q = np.array([[1],[2],[3]])
>>> x = linalg.solve_sylvester(a, b, q)
>>> x
array([[ 0.0625],
[-0.5625],
[ 0.6875]])
>>> np.allclose(a.dot(x) + x.dot(b), q)
True
scipy.linalg.solve_continuous_are
scipy.linalg.solve_continuous_are(a, b, q, r, e=None, s=None, balanced=True)
解连续时间代数 Riccati 方程(CARE)。
CARE 的定义为
[X A + A^H X - X B R^{-1} B^H X + Q = 0]
解的存在条件限制为:
- A 的所有特征值在右半平面上,应该是可控的。
- 关联的哈密顿笔(见注释),其特征值应足够远离虚轴。
此外,如果e
或s
不精确为None
,则 CARE 的广义版本
[E^HXA + A^HXE - (E^HXB + S) R^{-1} (B^HXE + S^H) + Q = 0]
被解决。当省略时,假设e
为单位矩阵,s
与a
和b
兼容且为零矩阵的大小相同。
参数:
a(M, M) array_like
方阵
b(M, N) array_like
输入
q(M, M) array_like
输入
r(N, N) array_like
非奇异方阵
e(M, M) array_like, 可选
非奇异方阵
s(M, N) array_like, 可选
输入
balancedbool, 可选
指示数据是否进行平衡步骤的布尔值,默认设置为 True。
返回:
x(M, M) ndarray
连续时间代数 Riccati 方程的解。
异常:
LinAlgError
对于无法分离出笔的稳定子空间的情况。请参阅 Notes 部分和详细信息的参考资料。
另见
解决离散时间代数 Riccati 方程
注释
方程式通过形成扩展的哈密顿矩阵笔来解决,如[1]中描述的,由块矩阵给出,[H - \lambda J]
[ A 0 B ] [ E 0 0 ]
[-Q -A^H -S ] - \lambda * [ 0 E^H 0 ]
[ S^H B^H R ] [ 0 0 0 ]
并使用 QZ 分解方法。
在此算法中,失败条件与产品(U_2 U_1^{-1})的对称性和(U_1)的条件数相关。这里,(U)是一个 2m×m 矩阵,包含了稳定子空间的特征向量,具有 2-m 行并分成两个 m 行矩阵。详见[1]和[2]获取更多详细信息。
为了提高 QZ 分解的准确性,笔在进行平衡步骤时,根据[3]中给出的配方,平衡绝对值的和(在删除对角元素后)。
从版本 0.11.0 开始新增。
参考文献
[1] (1,2)
P. van Dooren,《用于解 Riccati 方程的广义特征值方法》,SIAM 科学与统计计算杂志,Vol.2(2),DOI:10.1137/0902010
[2]
A.J. Laub,“用于解代数 Riccati 方程的 Schur 方法”,麻省理工学院。信息与决策系统实验室。LIDS-R ; 859。在线查看:hdl.handle.net/1721.1/1301
[3]
P. Benner,“哈密顿矩阵的辛平衡”,2001 年,SIAM J. Sci. Comput.,2001 年,Vol.22(5),DOI:10.1137/S1064827500367993
示例
给定 a, b, q, 和 r,解出 x:
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[4, 3], [-4.5, -3.5]])
>>> b = np.array([[1], [-1]])
>>> q = np.array([[9, 6], [6, 4.]])
>>> r = 1
>>> x = linalg.solve_continuous_are(a, b, q, r)
>>> x
array([[ 21.72792206, 14.48528137],
[ 14.48528137, 9.65685425]])
>>> np.allclose(a.T.dot(x) + x.dot(a)-x.dot(b).dot(b.T).dot(x), -q)
True
scipy.linalg.solve_discrete_are
scipy.linalg.solve_discrete_are(a, b, q, r, e=None, s=None, balanced=True)
解离散时间代数 Riccati 方程(DARE)。
DARE 定义为
[A^HXA - X - (A^HXB) (R + BHXB) (B^HXA) + Q = 0]
存在解的限制条件是:
- 所有 (A) 的特征值都在单位圆外,应该是可控的。
- 相关的辛特征对(见注释),其特征值应远离单位圆。
此外,如果 e
和 s
都不精确为 None
,则求解广义版本的 DARE
[A^HXA - E^HXE - (A^HXB+S) (R+BHXB) (BHXA+SH) + Q = 0]
被解决。当省略时,假定 e
为单位矩阵, s
为零矩阵。
参数:
a(M, M) 数组样式
方阵
b(M, N) 数组样式
输入
q(M, M) 数组样式
输入
r(N, N) 数组样式
方阵
e(M, M) 数组样式,可选
非奇异方阵
s(M, N) 数组样式,可选
输入
balanced布尔值
布尔值,指示是否在数据上执行平衡步骤。默认设置为 True。
返回:
x(M, M) ndarray
离散代数 Riccati 方程的解。
引发:
LinAlgError
对于无法隔离铅笔的稳定子空间的情况,请参见注释部分和详细的参考文献。
另请参阅
solve_continuous_are
解连续代数 Riccati 方程
注释
通过形成扩展辛矩阵对,求解方程 (H - \lambda J),如[1]所述,其中 (H - \lambda J)由块矩阵给出
[ A 0 B ] [ E 0 B ]
[ -Q E^H -S ] - \lambda * [ 0 A^H 0 ]
[ S^H 0 R ] [ 0 -B^H 0 ]
使用 QZ 分解方法。
在该算法中,失败条件与 (U_2 U_1^{-1}) 的对称性和 (U_1) 的条件数相关。这里,(U) 是一个 2m-by-m 矩阵,包含了稳定子空间的特征向量,具有 2-m 行,并分成两个 m 行的矩阵。详见[1]和[2]。
为了提高 QZ 分解的精度,铅笔经历了一个平衡步骤,其中绝对值的和(去除对角线条目后)按照[3]给出的配方平衡。如果数据有小的数值噪声,平衡可能会放大它们的影响,需要进行一些清理。
新功能在版本 0.11.0 中引入。
参考文献
[1] (1,2)
P. van Dooren,“用于解决 Riccati 方程的广义特征值方法”,SIAM 科学与统计计算杂志,Vol.2(2),DOI:10.1137/0902010
[2]
A.J. Laub, “用于解决代数 Riccati 方程的 Schur 方法”, 麻省理工学院. 信息与决策系统实验室. LIDS-R ; 859. 在线提供 : hdl.handle.net/1721.1/1301
[3]
P. Benner, “Hamiltonian 矩阵的辛平衡”, 2001, SIAM J. Sci. Comput., 2001, Vol.22(5), DOI:10.1137/S1064827500367993
例子
给定 a, b, q, 和 r 求解 x:
>>> import numpy as np
>>> from scipy import linalg as la
>>> a = np.array([[0, 1], [0, -1]])
>>> b = np.array([[1, 0], [2, 1]])
>>> q = np.array([[-4, -4], [-4, 7]])
>>> r = np.array([[9, 3], [3, 1]])
>>> x = la.solve_discrete_are(a, b, q, r)
>>> x
array([[-4., -4.],
[-4., 7.]])
>>> R = la.solve(r + b.T.dot(x).dot(b), b.T.dot(x).dot(a))
>>> np.allclose(a.T.dot(x).dot(a) - x - a.T.dot(x).dot(b).dot(R), -q)
True
scipy.linalg.solve_continuous_lyapunov
scipy.linalg.solve_continuous_lyapunov(a, q)
解决连续李亚普诺夫方程 (AX + XA^H = Q)。
使用巴特尔斯-斯图尔特算法找到 (X)。
参数:
aarray_like
方阵
qarray_like
右手边方阵
返回:
xndarray
连续李亚普诺夫方程的解
另请参阅
solve_discrete_lyapunov
计算离散时间李亚普诺夫方程的解
solve_sylvester
计算斯普尔斯特方程的解
注意
连续时间李亚普诺夫方程是斯普尔斯特方程的特殊形式,因此该解算器依赖于 LAPACK 例程 ?TRSYL。
新版本新增于 0.11.0。
示例
给定 a 和 q 解出 x:
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[-3, -2, 0], [-1, -1, 0], [0, -5, -1]])
>>> b = np.array([2, 4, -1])
>>> q = np.eye(3)
>>> x = linalg.solve_continuous_lyapunov(a, q)
>>> x
array([[ -0.75 , 0.875 , -3.75 ],
[ 0.875 , -1.375 , 5.3125],
[ -3.75 , 5.3125, -27.0625]])
>>> np.allclose(a.dot(x) + x.dot(a.T), q)
True
scipy.linalg.solve_discrete_lyapunov
scipy.linalg.solve_discrete_lyapunov(a, q, method=None)
解决离散 Lyapunov 方程 (AXA^H - X + Q = 0)。
参数:
a, q(M, M) array_like
对应于上述方程的 A 和 Q 的方阵。必须具有相同的形状。
方法,可选
求解器的类型。
如果未给出,则选择为 direct 如果 M
小于 10,否则为 bilinear
。
返回:
xndarray
离散 Lyapunov 方程的解
另见
solve_continuous_lyapunov
计算连续时间 Lyapunov 方程的解
注释
本节描述了可以通过 ‘method’ 参数选择的可用求解器。如果 M
小于 10,则默认方法为 direct,否则为 bilinear
。
方法 direct 使用直接的分析解来解离散 Lyapunov 方程。该算法在例如[1]中给出。然而,它要求线性解一个维度为 (M²) 的系统,因此即使对于中等大小的矩阵,性能也会迅速下降。
方法 bilinear 使用双线性变换将离散 Lyapunov 方程转换为连续 Lyapunov 方程 ((BX+XB'=-C)),其中 (B=(A-I)(A+I)^{-1}) 并且 (C=2(A' + I)^{-1} Q (A + I)^{-1})。连续方程可以有效地求解,因为它是 Sylvester 方程的特例。变换算法来自 Popov(1964),如[2]中描述。
自版本 0.11.0 新增。
参考文献
[1]
Hamilton, James D. Time Series Analysis, Princeton: Princeton University Press, 1994. 265. Print. doc1.lbfl.li/aca/FLMF037168.pdf
[2]
Gajic, Z., and M.T.J. Qureshi. 2008. Lyapunov Matrix Equation in System Stability and Control. Dover Books on Engineering Series. Dover Publications.
示例
给定 a 和 q 求解 x:
>>> import numpy as np
>>> from scipy import linalg
>>> a = np.array([[0.2, 0.5],[0.7, -0.9]])
>>> q = np.eye(2)
>>> x = linalg.solve_discrete_lyapunov(a, q)
>>> x
array([[ 0.70872893, 1.43518822],
[ 1.43518822, -2.4266315 ]])
>>> np.allclose(a.dot(x).dot(a.T)-x, -q)
True
scipy.linalg.clarkson_woodruff_transform
scipy.linalg.clarkson_woodruff_transform(input_matrix, sketch_size, seed=None)
应用 Clarkson-Woodruff 变换/草图到输入矩阵。
给定大小为 (n, d)
的输入矩阵 A
,计算大小为 (sketch_size, d)
的矩阵 A'
,以便
[|Ax| \approx |A'x|]
通过 Clarkson-Woodruff 变换,通常称为 CountSketch 矩阵,以高概率。
参数:
input_matrixarray_like
输入矩阵,形状为 (n, d)
。
sketch_sizeint
草图的行数。
seed{None, int, numpy.random.Generator
, numpy.random.RandomState
}, 可选
如果 seed 是 None(或 np.random),则使用 numpy.random.RandomState
单例。如果 seed 是一个整数,则使用新的带有 seed 种子的 RandomState
实例。如果 seed 已经是 Generator
或 RandomState
实例,则使用该实例。
返回:
A’array_like
对输入矩阵 A
的草图,大小为 (sketch_size, d)
。
注意事项
为了说明以下的结论
[|Ax| \approx |A'x|]
精确,观察以下的结果,它是从定理 14 的证明中适应的 [2] 通过马尔科夫不等式。如果我们有一个 sketch_size=k
的草图大小,它至少是
[k \geq \frac{2}{\epsilon²\delta}]
针对任意固定向量 x
,
[|Ax| = (1\pm\epsilon)|A'x|]
至少以概率1 - delta
。
此实现利用稀疏性:计算草图所需时间与 A.nnz
成正比。数据 A
以 scipy.sparse.csc_matrix
格式给出时,提供了稀疏输入的最快计算时间。
>>> import numpy as np
>>> from scipy import linalg
>>> from scipy import sparse
>>> rng = np.random.default_rng()
>>> n_rows, n_columns, density, sketch_n_rows = 15000, 100, 0.01, 200
>>> A = sparse.rand(n_rows, n_columns, density=density, format='csc')
>>> B = sparse.rand(n_rows, n_columns, density=density, format='csr')
>>> C = sparse.rand(n_rows, n_columns, density=density, format='coo')
>>> D = rng.standard_normal((n_rows, n_columns))
>>> SA = linalg.clarkson_woodruff_transform(A, sketch_n_rows) # fastest
>>> SB = linalg.clarkson_woodruff_transform(B, sketch_n_rows) # fast
>>> SC = linalg.clarkson_woodruff_transform(C, sketch_n_rows) # slower
>>> SD = linalg.clarkson_woodruff_transform(D, sketch_n_rows) # slowest
也就是说,在稠密输入上,这种方法表现良好,只是相对来说速度较慢。
参考资料
[1]
Kenneth L. Clarkson 和 David P. Woodruff。在 STOC, 2013 中的低秩逼近与输入稀疏时间回归。
[2]
David P. Woodruff。作为数值线性代数工具的草图化。在 Foundations and Trends in Theoretical Computer Science, 2014 中。
示例
创建一个大的密集矩阵 A
作为例子:
>>> import numpy as np
>>> from scipy import linalg
>>> n_rows, n_columns = 15000, 100
>>> rng = np.random.default_rng()
>>> A = rng.standard_normal((n_rows, n_columns))
应用变换来创建一个新的矩阵,其中有 200 行:
>>> sketch_n_rows = 200
>>> sketch = linalg.clarkson_woodruff_transform(A, sketch_n_rows, seed=rng)
>>> sketch.shape
(200, 100)
现在以高概率,真实范数的绝对值接近于草图范数。
>>> linalg.norm(A)
1224.2812927123198
>>> linalg.norm(sketch)
1226.518328407333
类似地,应用我们的草图保留了线性回归的解 (\min |Ax - b|)。
>>> b = rng.standard_normal(n_rows)
>>> x = linalg.lstsq(A, b)[0]
>>> Ab = np.hstack((A, b.reshape(-1, 1)))
>>> SAb = linalg.clarkson_woodruff_transform(Ab, sketch_n_rows, seed=rng)
>>> SA, Sb = SAb[:, :-1], SAb[:, -1]
>>> x_sketched = linalg.lstsq(SA, Sb)[0]
就像矩阵范数示例一样,linalg.norm(A @ x - b)
与高概率接近于 linalg.norm(A @ x_sketched - b)
。
>>> linalg.norm(A @ x - b)
122.83242365433877
>>> linalg.norm(A @ x_sketched - b)
166.58473879945151
scipy.linalg.block_diag
scipy.linalg.block_diag(*arrs)
从提供的数组创建块对角线矩阵。
给定输入的 A, B 和 C,输出将在对角线上排列这些数组。
[[A, 0, 0],
[0, B, 0],
[0, 0, C]]
参数:
A, B, C, … array_like,最多为二维
输入数组。长度为 n 的一维数组或类数组序列被视为形状为 (1,n)
的二维数组。
返回:
D ndarray
数组与 A, B, C, … 对角线上的元素。D 与 A 具有相同的数据类型。
注意事项
如果所有输入数组都是方阵,则输出称为块对角线矩阵。
空序列(即大小为零的类数组)不会被忽略。值得注意的是,[]
和 [[]]
都被视为形状为 (1,0)
的矩阵。
示例
>>> import numpy as np
>>> from scipy.linalg import block_diag
>>> A = [[1, 0],
... [0, 1]]
>>> B = [[3, 4, 5],
... [6, 7, 8]]
>>> C = [[7]]
>>> P = np.zeros((2, 0), dtype='int32')
>>> block_diag(A, B, C)
array([[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 3, 4, 5, 0],
[0, 0, 6, 7, 8, 0],
[0, 0, 0, 0, 0, 7]])
>>> block_diag(A, P, B, C)
array([[1, 0, 0, 0, 0, 0],
[0, 1, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 0, 0, 0, 0],
[0, 0, 3, 4, 5, 0],
[0, 0, 6, 7, 8, 0],
[0, 0, 0, 0, 0, 7]])
>>> block_diag(1.0, [2, 3], [[4, 5], [6, 7]])
array([[ 1., 0., 0., 0., 0.],
[ 0., 2., 3., 0., 0.],
[ 0., 0., 0., 4., 5.],
[ 0., 0., 0., 6., 7.]])
【推荐】国内首个AI IDE,深度理解中文开发场景,立即下载体验Trae
【推荐】编程新体验,更懂你的AI,立即体验豆包MarsCode编程助手
【推荐】抖音旗下AI助手豆包,你的智能百科全书,全免费不限次数
【推荐】轻量又高性能的 SSH 工具 IShell:AI 加持,快人一步
· 被坑几百块钱后,我竟然真的恢复了删除的微信聊天记录!
· 【自荐】一款简洁、开源的在线白板工具 Drawnix
· 没有Manus邀请码?试试免邀请码的MGX或者开源的OpenManus吧
· 园子的第一款AI主题卫衣上架——"HELLO! HOW CAN I ASSIST YOU TODAY
· 无需6万激活码!GitHub神秘组织3小时极速复刻Manus,手把手教你使用OpenManus搭建本