本文简要介绍 python 语言中 scipy.sparse.linalg.lsqr
的用法。
用法:
scipy.sparse.linalg.lsqr(A, b, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, iter_lim=None, show=False, calc_var=False, x0=None)#
找到大型稀疏线性方程组的最小二乘解。
该函数求解
Ax = b
或min ||Ax - b||^2
或min ||Ax - b||^2 + d^2 ||x - x0||^2
。矩阵A可以是正方形或矩形(over-determined或未确定),并且可以具有任何秩。
1. Unsymmetric equations -- solve Ax = b 2. Linear least squares -- solve Ax = b in the least-squares sense 3. Damped least squares -- solve ( A )*x = ( b ) ( damp*I ) ( damp*x0 ) in the least-squares sense
- A: {稀疏矩阵,ndarray,LinearOperator}
m-by-n 矩阵的表示。或者,
A
可以是一个线性运算符,它可以使用例如scipy.sparse.linalg.LinearOperator
生成Ax
和A^T x
。- b: 数组, 形状 (m,)
右侧向量
b
。- damp: 浮点数
阻尼系数。默认值为 0。
- atol, btol: 浮点数,可选
停止公差。
lsqr
继续迭代,直到某个后向误差估计小于某个数量,具体取决于 atol 和 btol。让r = b - Ax
是当前近似解的残差向量x
.如果Ax = b
似乎是一致的,lsqr
终止时norm(r) <= atol * norm(A) * norm(x) + btol * norm(b)
.否则,lsqr
终止时norm(A^H r) <= atol * norm(A) * norm(r)
.如果两个公差均为 1.0e-6(默认),则最终norm(r)
应该精确到大约 6 位数。 (决赛x
通常会有更少的正确数字,具体取决于cond(A)
和 LAMBDA 的大小。)如果环礁或者btol为无,将使用默认值 1.0e-6。理想情况下,它们应该是对条目中相对误差的估计A
和b
分别。例如,如果条目A
有 7 个正确的数字,设置atol = 1e-7
.这可以防止算法在输入数据的不确定性之外做不必要的工作。- conlim: 浮点数,可选
另一个停止公差。 lsqr 终止,如果估计
cond(A)
超过康林.对于兼容系统Ax = b
,康林可能大到 1.0e+12(比如说)。对于最小二乘问题,conlim 应小于 1.0e+8。通过设置可以获得最大精度atol = btol = conlim = zero
,但迭代次数可能过多。默认为 1e8。- iter_lim: 整数,可选
明确限制迭代次数(为了安全)。
- show: 布尔型,可选
显示迭代日志。默认为假。
- calc_var: 布尔型,可选
是否估计
(A'A + damp^2*I)^{-1}
的对角线。- x0: 数组,形状(n,),可选
如果没有使用零,则对 x 的初始猜测。默认为无。
- x: 浮点数
最终的解决方案。
- istop: int
给出终止的原因。 1 表示 x 是 Ax = b 的近似解。 2 表示 x 近似解决了最小二乘问题。
- itn: int
终止时的迭代次数。
- r1norm: 浮点数
norm(r)
,其中r = b - Ax
。- r2norm: 浮点数
sqrt( norm(r)^2 + damp^2 * norm(x - x0)^2 )
.等于r1范数如果damp == 0
.- anorm: 浮点数
Abar = [[A]; [damp*I]]
的 Frobenius 范数估计。- acond: 浮点数
cond(Abar)
的估计值。- arnorm: 浮点数
norm(A'@r - damp^2*(x - x0))
的估计值。- xnorm: 浮点数
norm(x)
- var: 浮点数
如果
calc_var
为 True,则估计(A'A)^{-1}
(如果damp == 0
)或更一般的(A'A + damp^2*I)^{-1}
的所有对角线。如果 A 具有完整的列排名或damp > 0
,则这是明确定义的。 (如果rank(A) < n
和damp = 0.
不确定 var 是什么意思)
参数 ::
返回 ::
注意:
LSQR 使用迭代方法来近似解。达到一定精度所需的迭代次数很大程度上取决于问题的规模。因此,应尽可能避免 A 的行或列的不良缩放。
例如,在问题 1 中,row-scaling 不会改变解决方案。如果 A 的一行与 A 的其他行相比非常小或大,则 (A b ) 的相应行应按比例放大或缩小。
在问题 1 和 2 中,解 x 在 column-scaling 之后很容易恢复。除非知道更好的信息,否则 A 的非零列应该被缩放,以便它们都具有相同的欧几里得范数(例如,1.0)。
在问题 3 中,如果湿气不为零,则没有重新缩放的自由。但是,只有在注意 A 的缩放之后才应该分配阻尼值。
参数阻尼旨在通过防止真正的解决方案非常大来帮助规范ill-conditioned 系统。参数 acond 提供了另一个对正则化的帮助,它可用于在计算的解变得非常大之前终止迭代。
如果某些初始估计
x0
已知并且如果damp == 0
,则可以进行如下操作:Compute a residual vector
r0 = b - A@x0
.Use LSQR to solve the system
A@dx = r0
.Add the correction dx to obtain a final solution
x = x0 + dx
.
这就要求
x0
在调用 LSQR 之前和之后可用。为了判断收益,假设 LSQR 需要 k1 次迭代来解决 A@x = b 和 k2 次迭代来解决A@dx= r0。如果 x0 是 “good”,则 norm(r0) 将小于 norm(b)。如果每个系统使用相同的停止公差 atol 和 btol,k1 和 k2 将相似,但最终解 x0 + dx 应该更准确。减少总功的唯一方法是对第二个系统使用更大的停止容差。如果某个值 btol 适用于 A@x = b,则较大的值 btol*norm(b)/norm(r0) 应该适用于A@dx= r0.预处理是另一种减少迭代次数的方法。如果可以有效地求解相关系统
M@x = b
,其中 M 以某种有用的方式逼近 A(例如,M - A 的秩较低或其元素相对于 A 的元素较小),则 LSQR 可能会更快地在系统上收敛A@M(inverse)@z = b
,之后可以通过求解 M@x = z 来恢复 x。如果 A 是对称的,则不应使用 LSQR!
替代方法是对称共轭梯度法 (cg) 和/或 SYMMLQ。 SYMMLQ 是对称 cg 的一种实现,适用于任何对称 A,并且比 LSQR 收敛得更快。如果 A 是正定的,则对称 cg 的其他实现每次迭代所需的工作量比 SYMMLQ 略少(但需要相同数量的迭代)。
参考:
[1]C. C. Paige 和 M. A. Saunders (1982a)。 “LSQR:稀疏线性方程和稀疏最小二乘的算法”,ACM TOMS 8(1),43-71。
[2]C. C. Paige 和 M. A. Saunders (1982b)。 “算法 583。LSQR:稀疏线性方程和最小二乘问题”,ACM TOMS 8(2),195-209。
[3]M. A. 桑德斯 (1995)。 “使用 LSQR 和 CRAIG 解决稀疏矩形系统”,BIT 35, 588-604。
例子:
>>> import numpy as np >>> from scipy.sparse import csc_matrix >>> from scipy.sparse.linalg import lsqr >>> A = csc_matrix([[1., 0.], [1., 1.], [0., 1.]], dtype=float)
第一个示例有简单的解决方案
[0, 0]
>>> b = np.array([0., 0., 0.], dtype=float) >>> x, istop, itn, normr = lsqr(A, b)[:4] >>> istop 0 >>> x array([ 0., 0.])
停止代码停止=0返回表示找到了一个由零组成的向量作为解。返回的解决方案x确实包含
[0., 0.]
。下一个示例有一个重要的解决方案:>>> b = np.array([1., 0., -1.], dtype=float) >>> x, istop, itn, r1norm = lsqr(A, b)[:4] >>> istop 1 >>> x array([ 1., -1.]) >>> itn 1 >>> r1norm 4.440892098500627e-16
如所示停止=1,
lsqr
找到了符合公差限制的解决方案。给定的解决方案[1., -1.]
显然解方程。其余返回值包括有关迭代次数的信息 (它=1) 以及求解方程左右两边的剩余差。最后一个示例演示了方程无解时的行为:>>> b = np.array([1., 0.01, -1.], dtype=float) >>> x, istop, itn, r1norm = lsqr(A, b)[:4] >>> istop 2 >>> x array([ 1.00333333, -0.99666667]) >>> A.dot(x)-b array([ 0.00333333, -0.00333333, 0.00333333]) >>> r1norm 0.005773502691896255
istop 表明系统是不一致的,因此 x 是相应最小二乘问题的近似解。 r1norm 包含找到的最小残差的范数。
相关用法
- Python SciPy linalg.lsmr用法及代码示例
- Python SciPy linalg.lstsq用法及代码示例
- Python SciPy linalg.lu_factor用法及代码示例
- Python SciPy linalg.lu_solve用法及代码示例
- Python SciPy linalg.lu用法及代码示例
- Python SciPy linalg.logm用法及代码示例
- Python SciPy linalg.ldl用法及代码示例
- Python SciPy linalg.leslie用法及代码示例
- Python SciPy linalg.lobpcg用法及代码示例
- Python SciPy linalg.lgmres用法及代码示例
- Python SciPy linalg.eigvalsh_tridiagonal用法及代码示例
- Python SciPy linalg.cdf2rdf用法及代码示例
- Python SciPy linalg.LaplacianNd用法及代码示例
- Python SciPy linalg.solve_circulant用法及代码示例
- Python SciPy linalg.polar用法及代码示例
- Python SciPy linalg.clarkson_woodruff_transform用法及代码示例
- Python SciPy linalg.rsf2csf用法及代码示例
- Python SciPy linalg.hessenberg用法及代码示例
- Python SciPy linalg.tril用法及代码示例
- Python SciPy linalg.triu用法及代码示例
- Python SciPy linalg.svd用法及代码示例
- Python SciPy linalg.ishermitian用法及代码示例
- Python SciPy linalg.invhilbert用法及代码示例
- Python SciPy linalg.factorized用法及代码示例
- Python SciPy linalg.SuperLU用法及代码示例
注:本文由纯净天空筛选整理自scipy.org大神的英文原创作品 scipy.sparse.linalg.lsqr。非经特殊声明,原始代码版权归原作者所有,本译文未经允许或授权,请勿转载或复制。