Ctrl+K
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 可以是方形或矩形（超定或欠定），并且可以具有任何秩。

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，线性算子}

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 似乎是一致的，则当 norm(r) <= atol * norm(A) * norm(x) + btol * norm(b) 时，lsqr 终止。 否则，当 norm(A^H r) <= atol * norm(A) * norm(r) 时，lsqr 终止。 如果两个容差均为 1.0e-6（默认值），则最终的 norm(r) 的精度应约为 6 位。（最终的 x 通常具有较少的正确位数，具体取决于 cond(A) 和 LAMBDA 的大小。） 如果 atol 或 btol 为 None，则将使用 1.0e-6 的默认值。 理想情况下，它们应该分别是 A 和 b 的条目中的相对误差的估计值。 例如，如果 A 的条目有 7 个正确的数字，则设置 atol = 1e-7。 这可防止算法执行超出输入数据不确定性的不必要的工作。

conlimfloat，可选

另一个停止容差。 如果 cond(A) 的估计值超过 conlim，则 lsqr 终止。 对于兼容系统 Ax = b，conlim 可以高达 1.0e+12（假设）。 对于最小二乘问题，conlim 应小于 1.0e+8。 可以通过设置 atol = btol = conlim = 零 来获得最大精度，但迭代次数可能会过多。 默认值为 1e8。

iter_limint，可选

对迭代次数的显式限制（为了安全）。

showbool，可选

显示迭代日志。 默认为 False。

calc_varbool，可选

是否估计 (A'A + damp^2*I)^{-1} 的对角线。

x0array_like，形状 (n,)，可选

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 的其他行相比非常小或非常大，则应向上或向下缩放 ( A b ) 的相应行。

在问题 1 和 2 中，列缩放后很容易恢复解 x。除非已知更好的信息，否则应缩放 A 的非零列，使它们都具有相同的欧几里得范数（例如，1.0）。

在问题 3 中，如果 damp 为非零值，则没有重新缩放的自由度。但是，只有在关注 A 的缩放比例之后，才应分配 damp 的值。

参数 damp 旨在通过防止真实解变得非常大来帮助正则化病态系统。参数 acond 也提供了另一种正则化辅助，可用于在计算出的解变得非常大之前终止迭代。

如果已知某些初始估计 x0 并且如果 damp == 0，则可以按如下方式进行

  1. 计算残差向量 r0 = b - A@x0。

  2. 使用 LSQR 求解系统 A@dx = r0。

  3. 添加校正 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: An algorithm for sparse linear equations and sparse least squares”, ACM TOMS 8(1), 43-71。

[2]

C. C. Paige 和 M. A. Saunders (1982b)。“Algorithm 583. LSQR: Sparse linear equations and least squares problems”, ACM TOMS 8(2), 195-209。

[3]

M. A. Saunders (1995)。“Solution of sparse rectangular systems using LSQR and 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 包含找到的最小残差的范数。

本页内容