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 包含找到的最小距離。