Wie löst NumPy kleinste Quadrate für unterbestimmte Systeme?

14

Nehmen wir an, wir haben X (2, 5)
und Y (2,)

Das funktioniert: np.linalg.lstsq(X, y)

Wir erwarten, dass dies nur dann funktioniert, wenn X die Form (N, 5) hat, wobei N> = 5 ist. Aber warum und wie?

Wir bekommen wie erwartet 5 Gewichte zurück, aber wie wird dieses Problem gelöst?

Ist es nicht so, als hätten wir 2 Gleichungen und 5 Unbekannte?
Wie könnte Numpy das lösen?
Es muss so etwas wie eine Interpolation sein, um künstlichere Gleichungen zu erstellen? ..

George Pligoropoulos
quelle
3
Warum sollte es nicht funktionieren? Ein unbestimmtes System hat viele Lösungen.
Matthew Gunn
Haben Sie vielleicht einen Link zu relevanter Theorie?
George Pligoropoulos

Antworten:

18

Nach meinem Verständnis basiert numpy.linalg.lstsq auf der LAPACK- Routine dgelsd .

Das Problem ist zu lösen:

minimize(overx)Axb2

Dies hat natürlich keine eindeutige Lösung für eine Matrix A, deren Rang kleiner als die Länge des Vektors b . dgelsdBietet im Falle eines unbestimmten Systems eine Lösung z , die Folgendes beinhaltet :

  • Az=b
  • z2x2 für allex , dieAx=b erfüllen. (dhz ist die minimale Normlösung für das unbestimmte System.

Beispiel, wenn System ist x+y=1 , gibt numpy.linalg.lstsqx=.5,y=.5 .

Wie funktioniert dgelsd?

Die Routine dgelsdberechnet die Singularwertzerlegung (SVD) von A.

Ich skizziere nur die Idee hinter der Verwendung einer SVD zur Lösung eines linearen Systems. Die Singulärwertzerlegung ist eine Faktorisierung UΣV=A wobeiU undV orthogonale Matrizen sind undΣ eine Diagonalmatrix ist, bei der die Diagonaleinträge als Singularwerte bezeichnet werden.

Der effektive Rang der Matrix A ist die Anzahl der singulären Werte, die effektiv nicht Null sind (dh sich ausreichend von Null in Bezug auf die Maschinengenauigkeit usw. unterscheiden). Sei S eine Diagonalmatrix der von Null verschiedenen Singularwerte. Die SVD ist also:

A=U[S000]V

Die Pseudoinverse von A ist gegeben durch:

A=V[S1000]U

Betrachten Sie die Lösung x=Ab . Dann:

Axb=U[S000]VV[S1000]Ubb=U[I000]Ubb

Grundsätzlich gibt es hier zwei Fälle:

  1. Die Anzahl der singulären Werte ungleich Null (dh die Größe der Matrix I ) ist kleiner als die Länge von b . Die Lösung hier wird nicht genau sein; Wir werden das lineare System im Sinne der kleinsten Quadrate lösen.
  2. Axb=0

Dieser letzte Teil ist etwas knifflig. Sie müssen die Matrixdimensionen verfolgen und verwenden, dass U eine orthogonale Matrix ist.

Äquivalenz von Pseudo-Inverse

Wenn A linear unabhängige Reihen hat (z. B. wir haben eine Fettmatrix), dann gilt:

A=A(AA)1

Für ein unbestimmtes System können Sie zeigen, dass das Pseudo-Inverse die minimale Normlösung ergibt.

Wenn A linear unabhängige Spalten hat (z. B. haben wir eine dünne Matrix), dann gilt:

A=(AA)1A

Matthew Gunn
quelle
dgelsd benutzt SVD aber R lm benutzt QR?
Haitao Du
@ hxd1011R lmverwendet standardmäßig die QR-Faktorisierung, Sie können jedoch Alternativen angeben.
Sycorax sagt Reinstate Monica