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)[source]#
找到大型稀疏线性方程组的最小二乘解。
该函数求解
Ax = b
或min ||Ax - b||^2
或min ||Ax - b||^2 + d^2 ||x - x0||^2
。矩阵 A 可以是方形或矩形(超定或欠定),并且可以具有任何秩。
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×n 矩阵的表示。或者,
A
可以是一个线性算子,它可以使用例如scipy.sparse.linalg.LinearOperator
生成Ax
和A^T x
。- barray_like, 形状 (m,)
右侧向量
b
。- dampfloat
阻尼系数。默认为 0。
- atol, btolfloat, optional
停止容差。
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 的大小。)如果 atol 或 btol 为 None,则将使用默认值 1.0e-6。理想情况下,它们应该是A
和b
的条目中相对误差的估计值。例如,如果A
的条目具有 7 个正确的数字,则设置atol = 1e-7
。这可以防止算法执行超出输入数据不确定性的不必要工作。- conlimfloat, optional
另一个停止容差。如果对
cond(A)
的估计超过 conlim,则 lsqr 终止。对于兼容系统Ax = b
,conlim 可能高达 1.0e+12(例如)。对于最小二乘问题,conlim 应小于 1.0e+8。可以通过设置atol = btol = conlim = zero
来获得最大精度,但迭代次数可能过多。默认值为 1e8。- iter_limint, optional
显式限制迭代次数(为了安全起见)。
- showbool, optional
显示迭代日志。默认为 False。
- calc_varbool, optional
是否估计
(A'A + damp^2*I)^{-1}
的对角线。- x0array_like, 形状 (n,), optional
x 的初始猜测,如果为 None,则使用零。默认为 None。
1.0.0 版本中新增。
- 返回:
- xfloat 的 ndarray
最终解。
- istopint
给出终止的原因。1 表示 x 是 Ax = b 的近似解。2 表示 x 近似解决了最小二乘问题。
- itnint
终止时的迭代次数。
- r1normfloat
norm(r)
,其中r = b - Ax
。- r2normfloat
sqrt( norm(r)^2 + damp^2 * norm(x - x0)^2 )
。如果damp == 0
,则等于 r1norm。- anormfloat
估计
Abar = [[A]; [damp*I]]
的 Frobenius 范数。- acondfloat
估计
cond(Abar)
。- arnormfloat
估计
norm(A'@r - damp^2*(x - x0))
。- xnormfloat
norm(x)
- varfloat 的 ndarray
如果
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 中,解不受行缩放的影响。如果 A 的一行与其他行相比非常小或大,则 ( A b ) 的相应行应向上或向下缩放。
在问题 1 和 2 中,在列缩放后很容易恢复解 x。除非已知更好的信息,否则 A 的非零列应进行缩放,使其都具有相同的欧几里德范数(例如,1.0)。
在问题 3 中,如果 damp 非零,则没有自由度来重新缩放。但是,只有在注意 A 的缩放之后,才应分配 damp 的值。
参数 damp 旨在通过防止真实解变得非常大来帮助正则化病态系统。参数 acond 也提供了正则化帮助,可用于在计算解变得非常大之前终止迭代。
如果已知一些初始估计值
x0
并且如果damp == 0
,则可以按如下方式进行计算残差向量
r0 = b - A@x0
。使用 LSQR 求解系统
A@dx = r0
。添加校正 dx 以获得最终解
x = x0 + dx
。
这要求在调用 LSQR 之前和之后都可以使用
x0
。为了判断好处,假设 LSQR 需要 k1 次迭代来求解 A@x = b,需要 k2 次迭代来求解 A@dx = r0。如果 x0 是“好的”,则 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. Saunders (1995)。“使用 LSQR 和 CRAIG 求解稀疏矩形系统”,BIT 35, 588-604。
示例
>>> import numpy as np >>> from scipy.sparse import csc_array >>> from scipy.sparse.linalg import lsqr >>> A = csc_array([[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.])
返回的停止代码
istop=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
如
istop=1
所示,lsqr
找到了一个满足容差限制的解。给定的解[1., -1.]
显然解决了该方程。其余返回值包括有关迭代次数 (itn=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 包含找到的最小残差的范数。