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)[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，可選

另一個停止容忍度。如果 cond(A) 的估計值超過 conlim，則 lsqr 終止。對於相容系統 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 版本中新增。

返回:
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。

這要求 x0 在呼叫 LSQR 之前和之後都可用。為了判斷好處，假設 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 包含找到的最小殘差的範數。

在本頁中