Ich versuche, den Algorithmus zur Anpassung der kleinsten Quadrate an Python zu implementieren, nachdem ich ihn bereits auf Matlab geschrieben habe. Allerdings habe ich Probleme, die richtige Transformationsmatrix zu erhalten, und das Problem scheint beim Lösungsschritt zu passieren. (Bearbeiten: Meine Transformationsmatrix ist unglaublich genau mit Matlab, aber komplett mit Python.)
Ich habe mir online zahlreiche Quellen angeschaut, und alle weisen darauf hin, dass man zur Übersetzung von Matlabs "mlldivide" "np.linalg.solve" verwenden muss, wenn die Matrix quadratisch und nichtsingulär ist, und "np.linalg.lstsq" ' Andernfalls. Doch meine Ergebnisse stimmen nicht überein.
Was ist das Problem? Wenn es mit der Implementierung der Funktionen zu tun hat, was ist die korrekte Übersetzung von mldivide in numpy?
Ich habe beide Versionen des Codes unten angehängt. Sie sind im Wesentlichen die exakt gleiche Implementierung, mit Ausnahme des Teils lösen.
Matlab-Code:
%Vor%Python-Code:
%Vor%Bearbeiten (gemäß einem Vorschlag):
%Vor% Interessant ist, dass Sie mit np.linalg.lstsq und np.linalg.solve ziemlich unterschiedliche Ergebnisse erhalten.
Beide sollten eine Lösung für die Gleichung Ax = B bieten. Diese ergeben jedoch zwei recht unterschiedliche Arrays:
%Vor%Beide sollten eine genaue (bis auf die Berechnungsgenauigkeit) Lösung für die Gruppe der linearen Gleichungen geben, und für eine nicht singuläre Matrix gibt es genau eine Lösung.
Etwas muss dann falsch sein. Lassen Sie uns sehen, ob diese beiden Kandidaten Lösungen für die ursprüngliche Gleichung sein könnten:
%Vor% Sie scheinen die gleiche Lösung zu geben, die im Wesentlichen die gleiche ist wie B_star , wie es sein sollte. Dies führt uns zur Erklärung. Mit einfacher linearer Algebra konnten wir voraussagen, dass A. (x1-x2) sollte sehr nahe bei Null liegen:
Und das ist es tatsächlich. Es scheint also, dass es eine nicht-triviale Lösung für die Gleichung Ax = 0 gibt, wobei die Lösung x = x1-x2 ist, was bedeutet, dass die Matrix singulär ist und somit eine unendliche Anzahl verschiedener Lösungen für Ax = B existiert.
Das Problem ist also nicht in NumPy oder Matlab, es ist in der Matrix selbst.
Im Falle dieser Matrix ist die Situation jedoch etwas schwierig. A_star scheint durch die obige Definition singulär zu sein (Ax = 0 für x & lt; & gt; 0). Auf der anderen Seite ist seine Determinante nicht Null und nicht singulär.
In diesem Fall ist A_star ein Beispiel für eine Matrix, die numerisch instabil und nicht singulär ist. Die Methode solve löst sie, indem sie die einfache Multiplikation durch Inverse verwendet. Dies ist in diesem Fall eine schlechte Wahl, da das Inverse sehr große und sehr kleine Werte enthält. Dies macht die Multiplikation anfällig für Rundungsfehler. Dies kann anhand der Zustandsnummer der Matrix gesehen werden:
Dies ist eine sehr hohe Bedingungsnummer, und die Matrix ist schlecht konditioniert.
In diesem Fall ist die Verwendung einer Umkehrung zur Lösung des Problems ein schlechter Ansatz. Der Ansatz der kleinsten Quadrate liefert bessere Ergebnisse, wie Sie vielleicht sehen. Die bessere Lösung besteht jedoch darin, die Eingabewerte so zu skalieren, dass x und x ^ 2 im gleichen Bereich liegen. Eine sehr gute Skalierung besteht darin, alles zwischen -1 und 1 zu skalieren.
Eine Sache, die Sie in Erwägung ziehen könnten, ist die Verwendung der Indexierungsfunktionen von NumPy. Zum Beispiel:
%Vor%entspricht:
%Vor% und Sie können Ihr Array A auf diese Weise konstruieren:
Sie müssen keine Listen mit Listen oder Schleifen erstellen.