Python SciPy linalg.lsqr用法及代碼示例

本文簡要介紹 python 語言中 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可以是正方形或矩形(over-determined或未確定)，並且可以具有任何秩。

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-by-n 矩陣的表示。或者，A 可以是一個線性運算符，它可以使用例如 scipy.sparse.linalg.LinearOperator 生成 Ax 和 A^T x。

b： 數組, 形狀 (m,)

右側向量b。

damp： 浮點數

阻尼係數。默認值為 0。

atol, btol： 浮點數，可選

停止公差。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 的大小。)如果環礁或者btol為無，將使用默認值 1.0e-6。理想情況下，它們應該是對條目中相對誤差的估計A和b分別。例如，如果條目A有 7 個正確的數字，設置atol = 1e-7.這可以防止算法在輸入數據的不確定性之外做不必要的工作。

conlim： 浮點數，可選

另一個停止公差。 lsqr 終止，如果估計cond(A)超過康林.對於兼容係統Ax = b,康林可能大到 1.0e+12(比如說)。對於最小二乘問題，conlim 應小於 1.0e+8。通過設置可以獲得最大精度atol = btol = conlim = zero，但迭代次數可能過多。默認為 1e8。

iter_lim： 整數，可選

明確限製迭代次數(為了安全)。

show： 布爾型，可選

顯示迭代日誌。默認為假。

calc_var： 布爾型，可選

是否估計 (A'A + damp^2*I)^{-1} 的對角線。

x0： 數組，形狀(n，)，可選

如果沒有使用零，則對 x 的初始猜測。默認為無。

返回 ：：

x： 浮點數

最終的解決方案。

istop： int

給出終止的原因。 1 表示 x 是 Ax = b 的近似解。 2 表示 x 近似解決了最小二乘問題。

itn： int

終止時的迭代次數。

r1norm： 浮點數

norm(r) ，其中 r = b - Ax 。

r2norm： 浮點數

sqrt( norm(r)^2  +  damp^2 * norm(x - x0)^2 ).等於r1範數如果damp == 0.

anorm： 浮點數

Abar = [[A]; [damp*I]] 的 Frobenius 範數估計。

acond： 浮點數

cond(Abar) 的估計值。

arnorm： 浮點數

norm(A'@r - damp^2*(x - x0)) 的估計值。

xnorm： 浮點數

norm(x)

var： 浮點數

如果 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 中，row-scaling 不會改變解決方案。如果 A 的一行與 A 的其他行相比非常小或大，則 (A b ) 的相應行應按比例放大或縮小。

在問題 1 和 2 中，解 x 在 column-scaling 之後很容易恢複。除非知道更好的信息，否則 A 的非零列應該被縮放，以便它們都具有相同的歐幾裏得範數(例如，1.0)。

在問題 3 中，如果濕氣不為零，則沒有重新縮放的自由。但是，隻有在注意 A 的縮放之後才應該分配阻尼值。

參數阻尼旨在通過防止真正的解決方案非常大來幫助規範ill-conditioned 係統。參數 acond 提供了另一個對正則化的幫助，它可用於在計算的解變得非常大之前終止迭代。

如果某些初始估計 x0 已知並且如果 damp == 0 ，則可以進行如下操作：

1. Compute a residual vector r0 = b - A@x0.

2. Use LSQR to solve the system A@dx = r0.

3. Add the correction dx to obtain a final solution x = x0 + dx.

這就要求x0在調用 LSQR 之前和之後可用。為了判斷收益，假設 LSQR 需要 k1 次迭代來解決 A@x = b 和 k2 次迭代來解決A@dx= r0。如果 x0 是 “good”，則 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. 桑德斯 (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.])

停止代碼停止=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

如所示停止=1,lsqr找到了符合公差限製的解決方案。給定的解決方案[1., -1.]顯然解方程。其餘返回值包括有關迭代次數的信息 (它=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 包含找到的最小殘差的範數。