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，可选
停止容差。
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，可选
另一个停止容差。lsqr 在
cond(A)的估计值超过 conlim 时终止。对于兼容系统Ax = b，conlim 可以高达 1.0e+12（例如）。对于最小二乘问题，conlim 应该小于 1.0e+8。可以通过设置atol = btol = conlim = zero来获得最大精度，但迭代次数可能过多。默认值为 1e8。- iter_limint，可选
对迭代次数的显式限制（为了安全起见）。
- showbool，可选
显示迭代日志。默认值为 False。
- calc_varbool，可选
是否估计
(A'A + damp^2*I)^{-1}的对角线。- x0array_like，形状 (n,)，可选
x 的初始猜测，如果为 None，则使用零。默认值为 None。
在版本 1.0.0 中添加。
- 返回值：
- xndarray of float
最终解。
- 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)- varndarray of float
如果
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 的其他行相比非常小或非常大，则 ( A b ) 的对应行应该进行缩放。
在问题 1 和 2 中，解 x 在列缩放后很容易恢复。除非知道更好的信息，否则应将 A 的非零列进行缩放，以使它们都具有相同的欧几里得范数（例如，1.0）。
在问题 3 中，如果 damp 不为零，则没有自由缩放的空间。但是，damp 的值应该在注意 A 的缩放之后再分配。
参数 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_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.])
返回的停止代码 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 包含找到的最小残差的范数。