Gauß-Newton-Verfahren

Grundidee und Einsatz

Das Gauß-Newton-Verfahren ist ein numerisches Verfahren zur Lösung nichtlinearer Minimierungsprobleme nach der Methode der kleinsten Quadrate. Es ist mit dem Newton-Verfahren verwandt, benötigt aber keine Berechnung der zweiten Ableitung. Das ist besonders bei großen Problemen mit mehreren zehntausend Parametern wichtig, weil die zweite Ableitung dort ein begrenzender Faktor sein kann.

Gesucht wird ein Parametervektor, der die Summe quadrierter Fehler möglichst klein macht. Dazu wird die nichtlineare Funktion zunächst lokal durch eine lineare Funktion angenähert. Diese Linearisierung wird anschließend mit der Methode der kleinsten Quadrate optimiert. Aus dem Ergebnis entsteht ein neuer Iterationspunkt.

Das Optimierungsproblem

Das Verfahren löst Probleme der Form

minₓ∈ℝⁿ { ½ ∑ᵢ₌₁ᵐ (fᵢ(x))² },

wobei fᵢ: ℝⁿ → ℝ stetig differenzierbare Funktionen und m ≥ n gilt. Mit der euklidischen Norm kann dasselbe Problem als

minₓ∈ℝⁿ { ½ ‖f(x)‖² }

geschrieben werden. Dabei ist f = (f₁, …, fₘ): ℝⁿ → ℝᵐ.

Das nichtlineare Gleichungsproblem f(x) = 0 ist unter der Voraussetzung, dass f eine Nullstelle besitzt, äquivalent zur Minimierung von ½‖f(x)‖². Das Verfahren wird außerdem häufig für nichtlineare Ausgleichsprobleme eingesetzt. In diesem Fall beschreiben die Funktionen fᵢ(x) die Abweichung einer Modellfunktion g(x) von bekannten Beobachtungen oder Messwerten yᵢ:

fᵢ(x) = g(x) − yᵢ.

Ist g eine lineare Abbildung, erhält man den Standardfall der Methode der kleinsten Quadrate mit linearer Modellfunktion.

Linearisierung und Iterationsschritt

Im Punkt x⁰ ∈ ℝⁿ wird f durch die Taylorentwicklung erster Ordnung angenähert:

f̃(x) = f(x⁰) + ∇f(x⁰)ᵀ(x − x⁰).

Die Matrix ∇f(x⁰)ᵀ wird als Jacobi-Matrix J bezeichnet. Durch Einsetzen der Linearisierung entsteht das lineare kleinste-Quadrate-Problem

minₓ∈ℝⁿ { ½‖J(x − x⁰) + f(x⁰)‖² }.

Der Gradient dieser Zielfunktion ist

Jᵀ(J(x − x⁰) + f(x⁰)).

Setzt man ihn gleich null, erhält man die Normalgleichungen

JᵀJ(x − x⁰) = −Jᵀf(x⁰).

Falls JᵀJ invertierbar ist, lautet die Lösung

x = x⁰ − (JᵀJ)⁻¹Jᵀf(x⁰).

Daraus folgt der Gauß-Newton-Iterationsschritt

xᵏ⁺¹ = xᵏ − αᵏ((J|ₓᵏ)ᵀJ|ₓᵏ)⁻¹(J|ₓᵏ)ᵀf(xᵏ).

Hier bedeutet J|ₓᵏ, dass die Jacobi-Matrix am Punkt xᵏ ausgewertet wird. αᵏ ≥ 0 ist die Schrittweite. In jedem Schritt wird also die aktuelle nichtlineare Funktion linearisiert, ein lineares kleinste-Quadrate-Problem gelöst und der neue Punkt berechnet.

Für das entstehende lineare Gleichungssystem gibt es je nach Größe und Struktur des Problems verschiedene Verfahren:

• Kleine Probleme mit n < 1000 und m < 10000 werden am besten mit der QR-Zerlegung gelöst.
• Für große Probleme eignet sich die Cholesky-Zerlegung, weil JᵀJ von Konstruktion aus symmetrisch ist. Für dünnbesetzte Matrizen JᵀJ gibt es speziell angepasste Cholesky-Varianten.
• Allgemein kann das CG-Verfahren verwendet werden; dabei ist üblicherweise eine Vorkonditionierung notwendig.

Konvergenz und Stabilisierung

Der Update-Vektor hat die Form

d = −DJᵀf(x), wobei D = (JᵀJ)⁻¹.

Hat J vollen Rang, dann sind JᵀJ und D positiv definit. Da Jᵀf(x) der Gradient des quadratischen Problems ½‖f(x)‖² ist, ist d eine Abstiegsrichtung. Es gilt

(Jᵀf(x))ᵀd < 0.

Bei geeigneter Wahl der Schrittweite αᵏ folgt daraus die Konvergenz zu einem stationären Punkt. Das Verfahren kann daher im Wesentlichen als skaliertes Gradientenverfahren mit der positiv definiten Skalierungsmatrix D verstanden werden.

Eine allgemeine Aussage über die Konvergenzgeschwindigkeit ist nicht möglich. Liegt der Startpunkt x⁰ sehr weit vom Optimum entfernt oder ist JᵀJ schlecht konditioniert, kann die Konvergenz unter Umständen nur sublinear sein. In vielen praktischen Anwendungen ist sie jedoch deutlich schneller und kann in bestimmten Fällen sogar die quadratische Konvergenz des Newton-Verfahrens erreichen.

Die Taylorentwicklung zweiter Ordnung der Zielfunktion wird beschrieben durch

f̃(x) ≈ f(x⁰) + ∇f(x⁰)ᵀ(x − x⁰) + ½(x − x⁰)ᵀH(x⁰)(x − x⁰),

wobei H(x⁰) die Hesse-Matrix im Punkt x⁰ ist. Ist H(x) klein, etwa weil f(x) fast linear ist oder die Komponentenfunktionen fᵢ(x) in der Nähe des Optimums sehr klein sind, kann der quadratische Term vernachlässigt werden; dann konvergiert das Gauß-Newton-Verfahren superlinear. Gilt im optimalen Punkt x* außerdem H(x*) = 0, ist der entsprechende Newton-Schritt identisch mit dem Gauß-Newton-Schritt, und die Konvergenz ist quadratisch.

Bei schlecht konditioniertem oder singulärem JᵀJ kann der Iterationsschritt durch eine Diagonalmatrix Δᵏ stabilisiert werden:

xᵏ⁺¹ = xᵏ − αᵏ((J|ₓᵏ)ᵀJ|ₓᵏ + Δᵏ)⁻¹(J|ₓᵏ)ᵀf(xᵏ).

Δᵏ wird so gewählt, dass J|ₓᵏᵀJ|ₓᵏ + Δᵏ positiv definit ist. Mit Δᵏ = βᵏI und β ≥ 0, also einem skalaren Vielfachen der Identitätsmatrix, erhält man den Levenberg-Marquardt-Algorithmus.

Beispiel: Rosenbrock-Funktion

Die Rosenbrock-Funktion ist ein typischer Test für Optimierungsverfahren, weil sie ein schmales und flaches Tal besitzt, in dem iterative Methoden nur kleine Schritte machen können. Sie lautet

g: ℝ² → ℝ: x ↦ (a − x₁)² + b(x₂ − x₁²)².

Üblicherweise werden a = 1 und b = 100 gewählt. Dann liegt das globale Optimum bei x* = (1, 1) mit g(x*) = 0.

Damit die Funktion als Summe von Quadraten dargestellt werden kann, setzt man

f₁(x) = √2(a − x₁),

f₂(x) = √(2b)(x₂ − x₁²).

Damit lautet das Gauß-Newton-Problem minₓ∈ℝ² ½‖f(x)‖² mit

f(x) = (√2(a − x₁), √(2b)(x₂ − x₁²))ᵀ.

Die Jacobi-Matrix ist

J = [−√2, 0; −2x₁√(2b), √(2b)],

und

JᵀJ = [8bx₁² + 2, −4bx₁; −4bx₁, 2b].

Da J vollen Rang hat, ist JᵀJ positiv definit und seine Inverse existiert.

Für die Schrittweite wird eine einfache Liniensuche verwendet: Man startet mit αᵏ = 1, berechnet den Kandidaten x̃ = xᵏ + αᵏd und akzeptiert ihn, wenn g(x̃) < g(xᵏ) gilt. Andernfalls wird αᵏ halbiert und der Kandidat erneut berechnet. Diese Suche garantiert wegen der Abstiegsrichtung d einen kleineren Funktionswert, gegebenenfalls mit sehr kleiner Schrittweite.

Beim Startpunkt x⁰ = (0, −0.1) erreicht das Verfahren das Optimum in wenigen Iterationen. Die Funktionswerte sinken von 2 über 1.8291, 1.6306, 1.6131, 1.3000, 1.0300 und 0.2150 bis zu 1.1212e−27 beim Punkt (1.0, 1.0) in Iteration 7.

Zum Vergleich erreicht das Gradientenverfahren mit derselben Liniensuche selbst nach 500 Iterationen das Optimum nicht. Es endet bei ungefähr (0.8513, 0.7233) mit dem Funktionswert 0.0223. Das Beispiel zeigt den möglichen Geschwindigkeitsvorteil des Gauß-Newton-Verfahrens bei dieser nichtlinearen kleinste-Quadrate-Aufgabe.