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)

Findet die kleinsten Quadrate Lösung eines großen, dünnbesetzten linearen Gleichungssystems.

Die Funktion löst Ax = b oder min ||Ax - b||^2 oder min ||Ax - b||^2 + d^2 ||x - x0||^2.

Die Matrix A kann quadratisch oder rechteckig sein (überbestimmt oder unterbestimmt) und kann beliebigen Rang haben.

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

Parameter:

A{sparse array, ndarray, LinearOperator}

Repräsentation einer m-mal-n-Matrix. Alternativ kann A ein linearer Operator sein, der Ax und A^T x erzeugen kann, z. B. über scipy.sparse.linalg.LinearOperator.

barray_like, shape (m,)

Vektor der rechten Seite b.

dampfloat

Dämpfungskoeffizient. Standard ist 0.

atol, btolfloat, optional

Abbruch-Toleranzen. lsqr setzt die Iterationen fort, bis ein bestimmter Fehler in Rückwärtsrichtung kleiner als eine Größe ist, die von atol und btol abhängt. Sei r = b - Ax der Residuenvektor für die aktuelle ungefähre Lösung x. Wenn Ax = b konsistent zu sein scheint, terminiert lsqr, wenn norm(r) <= atol * norm(A) * norm(x) + btol * norm(b). Andernfalls terminiert lsqr, wenn norm(A^H r) <= atol * norm(A) * norm(r). Wenn beide Toleranzen 1,0e-6 (Standard) betragen, sollte die finale norm(r) auf etwa 6 Stellen genau sein. (Die finale x wird normalerweise weniger korrekte Stellen haben, abhängig von cond(A) und der Größe von LAMBDA.) Wenn atol oder btol None ist, wird ein Standardwert von 1,0e-6 verwendet. Idealerweise sollten sie Schätzungen des relativen Fehlers in den Einträgen von A bzw. b sein. Wenn beispielsweise die Einträge von A 7 korrekte Ziffern haben, setzen Sie atol = 1e-7. Dies verhindert, dass der Algorithmus unnötige Arbeit über die Unsicherheit der Eingabedaten hinaus leistet.

conlimfloat, optional

Eine weitere Abbruch-Toleranz. lsqr bricht ab, wenn eine Schätzung von cond(A) conlim überschreitet. Für konsistente Systeme Ax = b könnte conlim bis zu 1,0e+12 (z. B.) betragen. Für Probleme mit kleinsten Quadraten sollte conlim kleiner als 1,0e+8 sein. Maximale Präzision kann durch Setzen von atol = btol = conlim = null erreicht werden, aber die Anzahl der Iterationen kann dann exzessiv sein. Standard ist 1e8.

iter_limint, optional

Explizite Begrenzung der Anzahl der Iterationen (zur Sicherheit).

showbool, optional

Anzeige eines Iterationsprotokolls. Standard ist False.

calc_varbool, optional

Ob die Diagonalen von (A'A + damp^2*I)^{-1} geschätzt werden sollen.

x0array_like, shape (n,), optional

Anfangsschätzung von x, bei None werden Nullen verwendet. Standard ist None.

Hinzugefügt in Version 1.0.0.

Rückgabe:

xndarray von float

Die finale Lösung.

istopint

Gibt den Grund für die Beendigung an. 1 bedeutet, dass x eine ungefähre Lösung von Ax = b ist. 2 bedeutet, dass x das Problem der kleinsten Quadrate ungefähr löst.

itnint

Iteraionsnummer bei Beendigung.

r1normfloat

norm(r), wobei r = b - Ax.

r2normfloat

sqrt( norm(r)^2  +  damp^2 * norm(x - x0)^2 ). Gleich r1norm, wenn damp == 0.

anormfloat

Schätzung der Frobenius-Norm von Abar = [[A]; [damp*I]].

acondfloat

Schätzung von cond(Abar).

arnormfloat

Schätzung von norm(A'@r - damp^2*(x - x0)).

xnormfloat

norm(x)

varndarray von float

Wenn calc_var True ist, werden alle Diagonalen von (A'A)^{-1} (wenn damp == 0) oder allgemeiner von (A'A + damp^2*I)^{-1} geschätzt. Dies ist gut definiert, wenn A vollen Spaltenrang hat oder damp > 0. (Nicht sicher, was var bedeutet, wenn rank(A) < n und damp = 0.)

Hinweise

LSQR verwendet eine iterative Methode, um die Lösung zu approximieren. Die Anzahl der Iterationen, die zur Erreichung einer bestimmten Genauigkeit erforderlich sind, hängt stark von der Skalierung des Problems ab. Eine schlechte Skalierung der Zeilen oder Spalten von A sollte daher, wo immer möglich, vermieden werden.

Beispielsweise bleibt die Lösung im ersten Beispiel von der Zeilenskalierung unverändert. Wenn eine Zeile von A sehr klein oder groß im Vergleich zu den anderen Zeilen von A ist, sollte die entsprechende Zeile von ( A b ) hoch- oder runterskaliert werden.

In den Problemen 1 und 2 lässt sich die Lösung x nach der Spaltenskalierung leicht wiederherstellen. Sofern keine besseren Informationen vorliegen, sollten die von Null verschiedenen Spalten von A so skaliert werden, dass sie alle die gleiche euklidische Norm haben (z. B. 1,0).

In Problem 3 gibt es keine Freiheit zur Neuskalierung, wenn damp ungleich Null ist. Der Wert von damp sollte jedoch erst zugewiesen werden, nachdem die Skalierung von A berücksichtigt wurde.

Der Parameter damp soll zur Regularisierung schlecht konditionierter Systeme beitragen, indem er verhindert, dass die wahre Lösung sehr groß wird. Eine weitere Hilfe zur Regularisierung bietet der Parameter acond, der verwendet werden kann, um Iterationen zu beenden, bevor die berechnete Lösung sehr groß wird.

Wenn eine anfängliche Schätzung x0 bekannt ist und damp == 0 ist, könnte man wie folgt vorgehen:

1. Berechnen Sie einen Residuenvektor r0 = b - A@x0.

2. Verwenden Sie LSQR, um das System A@dx = r0 zu lösen.

3. Addieren Sie die Korrektur dx, um eine endgültige Lösung x = x0 + dx zu erhalten.

Dies erfordert, dass x0 vor und nach dem Aufruf von LSQR verfügbar ist. Um die Vorteile zu beurteilen, nehmen wir an, LSQR benötigt k1 Iterationen, um A@x = b zu lösen, und k2 Iterationen, um A@dx = r0 zu lösen. Wenn x0 "gut" ist, wird norm(r0) kleiner als norm(b) sein. Wenn die gleichen Abbruch-Toleranzen atol und btol für jedes System verwendet werden, werden k1 und k2 ähnlich sein, aber die endgültige Lösung x0 + dx sollte genauer sein. Der einzige Weg, die Gesamtarbeit zu reduzieren, ist die Verwendung einer größeren Abbruch-Toleranz für das zweite System. Wenn ein Wert btol für A@x = b geeignet ist, sollte der größere Wert btol*norm(b)/norm(r0) für A@dx = r0 geeignet sein.

Vorkonditionierung ist eine weitere Möglichkeit, die Anzahl der Iterationen zu reduzieren. Wenn es möglich ist, ein verwandtes System M@x = b effizient zu lösen, wobei M A auf hilfreiche Weise approximiert (z. B. M - A hat niedrigen Rang oder seine Elemente sind klein im Verhältnis zu denen von A), kann LSQR auf dem System A@M(inverse)@z = b schneller konvergieren, wonach x durch Lösen von M@x = z wiederhergestellt werden kann.

Wenn A symmetrisch ist, sollte LSQR nicht verwendet werden!

Alternativen sind die symmetrische Konjugierten-Gradienten-Methode (cg) und/oder SYMMLQ. SYMMLQ ist eine Implementierung von symmetrischem cg, die für jedes symmetrische A gilt und schneller konvergiert als LSQR. Wenn A positiv definit ist, gibt es andere Implementierungen von symmetrischem cg, die pro Iteration etwas weniger Arbeit erfordern als SYMMLQ (aber die gleiche Anzahl von Iterationen benötigen).

Referenzen

[1]

C. C. Paige und 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 und 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.

Beispiele

Probieren Sie es in Ihrem Browser aus!

>>> 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)

Das erste Beispiel hat die triviale Lösung [0, 0].

>>> b = np.array([0., 0., 0.], dtype=float)
>>> x, istop, itn, normr = lsqr(A, b)[:4]
>>> istop
0
>>> x
array([ 0.,  0.])

Der Abbruchcode istop=0 zeigt an, dass ein Nullvektor als Lösung gefunden wurde. Die zurückgegebene Lösung x enthält tatsächlich [0., 0.]. Das nächste Beispiel hat eine nicht-triviale Lösung.

>>> 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

Wie durch istop=1 angezeigt, hat lsqr eine Lösung gefunden, die den Toleranzgrenzen entspricht. Die gegebene Lösung [1., -1.] löst offensichtlich die Gleichung. Die restlichen Rückgabewerte enthalten Informationen über die Anzahl der Iterationen (itn=1) und die verbleibende Differenz zwischen linker und rechter Seite der gelösten Gleichung. Das letzte Beispiel demonstriert das Verhalten für den Fall, dass keine Lösung für die Gleichung existiert.

>>> 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 zeigt an, dass das System inkonsistent ist und somit x eher eine ungefähre Lösung des entsprechenden Problems der kleinsten Quadrate ist. r1norm enthält die Norm des minimalen gefundenen Residuums.

ZurückIn Tab öffnen