Python SciPy linalg.lsmr用法及代码示例

本文简要介绍 python 语言中 scipy.sparse.linalg.lsmr 的用法。

用法:

scipy.sparse.linalg.lsmr(A, b, damp=0.0, atol=1e-06, btol=1e-06, conlim=100000000.0, maxiter=None, show=False, x0=None)#

最小二乘问题的迭代求解器。

lsmr 求解线性方程组 Ax = b 。如果系统不一致，则解决最小二乘问题min ||b - Ax||_2。 A 是维度为 m-by-n 的矩形矩阵，其中允许所有情况：m = n、m > n 或 m < n。 b是长度为m的向量。矩阵 A 可以是稠密的，也可以是稀疏的(通常是稀疏的)。

参数 ：：

A： {稀疏矩阵，ndarray，LinearOperator}

线性系统中的矩阵 A。或者，A 可以是一个线性运算符，它可以使用例如 scipy.sparse.linalg.LinearOperator 生成 Ax 和 A^H x。

b： 数组, 形状 (m,)

线性系统中的向量b。

damp： 浮点数

正则化最小二乘的阻尼因子。 lsmr 解决正则化最小二乘问题：

min ||(b) - (  A   )x||
    ||(0)   (damp*I) ||_2

其中潮湿是一个标量。如果damp 为None 或0，则系统在没有正则化的情况下求解。默认值为 0。

atol, btol： 浮点数，可选

停止公差。lsmr继续迭代，直到某个后向误差估计小于某个数量，具体取决于 atol 和 btol。让r = b - Ax是当前近似解的残差向量x.如果Ax = b似乎是一致的，lsmr终止时norm(r) <= atol * norm(A) * norm(x) + btol * norm(b).否则，lsmr终止时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： 浮点数，可选

lsmr终止，如果估计cond(A)超过康林.对于兼容系统Ax = b，conlim 可能大到 1.0e+12(比如说)。对于最小二乘问题，康林应小于 1.0e+8。如果康林为无，默认值为 1e+8。可以通过设置获得最大精度atol = btol = conlim = 0，但迭代次数可能过多。默认为 1e8。

maxiter： 整数，可选

lsmr如果迭代次数达到则终止马克西特.默认是maxiter = min(m, n).对于ill-conditioned 系统，较大的值马克西特可能需要。默认为假。

show： 布尔型，可选

如果 show=True 打印迭代日志。默认为假。

x0： 数组，形状(n，)，可选

x 的初始猜测，如果没有使用零。默认为无。

返回 ：：

x： 浮点数

Least-square 解决方案返回。

istop： int

istop 给出了停止的原因：

istop   = 0 means x=0 is a solution.  If x0 was given, then x=x0 is a
            solution.
        = 1 means x is an approximate solution to A@x = B,
            according to atol and btol.
        = 2 means x approximately solves the least-squares problem
            according to atol.
        = 3 means COND(A) seems to be greater than CONLIM.
        = 4 is the same as 1 with atol = btol = eps (machine
            precision)
        = 5 is the same as 2 with atol = eps.
        = 6 is the same as 3 with CONLIM = 1/eps.
        = 7 means ITN reached maxiter before the other stopping
            conditions were satisfied.

itn： int

使用的迭代次数。

normr： 浮点数

norm(b-Ax)

normar： 浮点数

norm(A^H (b - Ax))

norma： 浮点数

norm(A)

conda： 浮点数

A的条件编号。

normx： 浮点数

norm(x)

注意：

参考：

[1]

D.C.-L。 Fong 和 M. A. Saunders，“LSMR：稀疏最小二乘问题的迭代算法”，SIAM J. Sci。计算机，卷。 33，第 2950-2971 页，2011 年。arXiv:1006.0758

[2]

LSMR 软件，https://web.stanford.edu/group/SOL/software/lsmr/

例子：

>>> import numpy as np
>>> from scipy.sparse import csc_matrix
>>> from scipy.sparse.linalg import lsmr
>>> 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 = lsmr(A, b)[:4]
>>> istop
0
>>> x
array([0., 0.])

停止代码停止=0返回表示找到了一个由零组成的向量作为解。返回的解决方案x确实包含[0., 0.]。下一个示例有一个重要的解决方案：

>>> b = np.array([1., 0., -1.], dtype=float)
>>> x, istop, itn, normr = lsmr(A, b)[:4]
>>> istop
1
>>> x
array([ 1., -1.])
>>> itn
1
>>> normr
4.440892098500627e-16

如所示停止=1,lsmr找到了符合公差限制的解决方案。给定的解决方案[1., -1.]显然解方程。其余返回值包括有关迭代次数的信息 (它=1) 以及求解方程左右两边的剩余差。最后一个示例演示了方程无解时的行为：

>>> b = np.array([1., 0.01, -1.], dtype=float)
>>> x, istop, itn, normr = lsmr(A, b)[:4]
>>> istop
2
>>> x
array([ 1.00333333, -0.99666667])
>>> A.dot(x)-b
array([ 0.00333333, -0.00333333,  0.00333333])
>>> normr
0.005773502691896255

istop 表明系统是不一致的，因此 x 是相应最小二乘问题的近似解。 norr 包含找到的最小距离。