Wie findet man das Minimum einer nichtlinearen, multivariaten Funktion mit der Newton-Methode (Code, nicht lineare Algebra)?

Question

Ich versuche, eine Parameterschätzung durchzuführen und möchte Parameterschätzungen wählen, die den quadratischen Fehler in einer vorhergesagten Gleichung minimieren über etwa 30 Variablen . Wäre die Gleichung linear, würde ich einfach die 30 partiellen Ableitungen berechnen, sie alle auf Null setzen und einen Löser für lineare Gleichungen verwenden. Aber leider die Gleichung ist nichtlinear und ihre Ableitungen sind es auch.

Wäre die Gleichung über eine einzige Variable, würde ich einfach Folgendes verwenden Newton-Verfahren (auch bekannt als Newton-Raphson). Das Internet ist reich an Beispielen und Code zur Umsetzung der Newton-Methode für Funktionen mit einer einzigen Variablen .

Da ich etwa 30 Variablen habe, wie kann ich eine numerische Lösung für dieses Problem mit der Newton-Methode programmieren? ? Ich habe die Gleichung in geschlossener Form und kann die erste und zweite Ableitung berechnen, aber ich weiß nicht so recht, wie es weitergehen soll. Ich habe im Internet eine große Anzahl von Behandlungen gefunden, aber sie gehen schnell in die schwere Matrixnotation über. Ich habe gefunden etwas mäßig Hilfreiches auf Wikipedia, aber ich habe Probleme, es in Code zu übersetzen.

Ich befürchte, dass ich bei der Matrixalgebra und den Matrixinvertierungen nicht weiterkomme. Ich kann eine Matrix mit einem Löser für lineare Gleichungen invertieren, aber ich mache mir Sorgen um die richtigen Zeilen und Spalten, die Vermeidung von Transpositionsfehlern und so weiter.

Um ganz konkret zu sein:

• Ich möchte mit Tabellen arbeiten, die Variablen auf ihre Werte abbilden. Ich kann eine Funktion für eine solche Tabelle schreiben, die den quadratischen Fehler bei einer solchen Tabelle als Argument zurückgibt. Ich kann auch Funktionen erstellen, die eine partielle Ableitung in Bezug auf eine bestimmte Variable zurückgeben.

• Ich habe eine vernünftige Ausgangsschätzung für die Werte in der Tabelle, daher mache ich mir keine Sorgen um die Konvergenz.

• Ich bin mir nicht sicher, wie ich die Schleife schreiben soll, die eine Schätzung (Wertetabelle für jede Variable), die Funktion und eine Tabelle mit partiell abgeleiteten Funktionen verwendet, um eine neue Schätzung zu erstellen.

Bei letzterem bräuchte ich Hilfe. Jede direkte Hilfe oder Hinweise auf gute Quellen werden sehr geschätzt.

---

Edit: Da ich die erste und zweite Ableitung in geschlossener Form habe, möchte ich sie nutzen und langsamer konvergierende Methoden wie Simplex-Suchen vermeiden.

Answer 1

Vielleicht denken Sie, dass Sie eine gute Lösung haben, aber für mich ist es am einfachsten, das Problem zunächst im Fall von einer Variablen zu verstehen und es dann auf den Matrixfall zu übertragen.

Wenn Sie im Fall der 1-Variablen die erste Ableitung durch die zweite Ableitung teilen, erhalten Sie die (negative) Schrittweite bis zum nächsten Versuchspunkt, z. B. -V/A.

Im Fall von N-Variablen ist die erste Ableitung ein Vektor und die zweite Ableitung eine Matrix (die Hessian). Man multipliziert den Ableitungsvektor mit dem Kehrwert der zweiten Ableitung, und das Ergebnis ist der negative Schrittvektor zum nächsten Versuchspunkt, z. B. -V*(1/A)

Ich gehe davon aus, dass Sie die hessische Matrix der 2. Ableitung erhalten können. Sie benötigen eine Routine, um sie zu invertieren. Davon gibt es viele in verschiedenen linearen Algebra-Paketen, und sie sind recht schnell.

(Für Leser, die mit dieser Idee nicht vertraut sind: Nehmen wir an, die beiden Variablen sind x und y, und die Oberfläche ist v(x,y). Dann ist die erste Ableitung der Vektor:

 V = [ dv/dx, dv/dy ]

und die zweite Ableitung ist die Matrix:

 A = [dV/dx]
    [dV/dy]

oder:

 A = [ d(dv/dx)/dx, d(dv/dy)/dx]
    [ d(dv/dx)/dy, d(dv/dy)/dy]

oder:

 A = [d^2v/dx^2, d^2v/dydx]
    [d^2v/dxdy, d^2v/dy^2]

die symmetrisch ist).

Wenn die Oberfläche parabolisch ist (konstante 2. Ableitung), wird die Antwort in einem Schritt erreicht. Wenn die 2. Ableitung jedoch nicht konstant ist, kann es zu Schwingungen kommen. Wenn man jeden Schritt halbiert (oder einen Bruchteil davon), sollte das Ergebnis stabil sein.

Wenn N == 1 ist, werden Sie sehen, dass es das Gleiche tut wie im Fall der 1-Variablen.

Viel Glück!

Hinzugefügt: Sie wollten einen Code:

 double X[N];
// Set X to initial estimate
while(!done){
    double V[N];    // 1st derivative "velocity" vector
    double A[N*N];  // 2nd derivative "acceleration" matrix
    double A1[N*N]; // inverse of A
    double S[N];    // step vector
    CalculateFirstDerivative(V, X);
    CalculateSecondDerivative(A, X);
    // A1 = 1/A
    GetMatrixInverse(A, A1);
    // S = V*(1/A)
    VectorTimesMatrix(V, A1, S);
    // if S is small enough, stop
    // X -= S
    VectorMinusVector(X, S, X);
}

Answer 2

Meiner Meinung nach sollte ein stochastischer Optimierer verwendet werden, z. B. ein Partikelschwarmverfahren.