Numerische Mathematik 1 – Nicht-lineare Gleichungen und Optimierung

Nullstellen von Funktionen

Bisektionsverfahren

Sei \(f\colon [a, b] \rightarrow \real \) eine stetige Funktion mit \(f(a) f(b) \le 0\).
Nach dem Zwischenwertsatz besitzt \(f\) mindestens eine Nullstelle in \([a, b]\).

Halbiert man das Intervall und wertet \(f\) an der Intervallmitte \(c = \frac {a + b}{2}\) aus, so kann man mithilfe des Vorzeichens entscheiden, in welchem Teilintervall eine Nullstelle liegen muss:

\begin{align*} f(a)f(c) \le 0 \quad \Rightarrow & \quad \text {es gibt eine Nullstelle in } [a, c] \\ f(a)f(c) \ge 0 \quad \Rightarrow & \quad \text {es gibt eine Nullstelle in } [c, b]. \end{align*} Nun wählt man das entsprechende Teilintervall aus und iteriert das Verfahren, bis die Länge des Intervalls die gewünschte Genauigkeit erreicht hat.

Dieses Verfahren zur Nullstellenbestimmung heißt Bisektionsverfahren (Zweiteilung).
Die Schnelligkeit der Konvergenz eines Iterationsverfahrens kann mit der Konvergenzordnung \(p \ge 1\) beurteilt werden. Ist der \(\ell \)-te Fehler \(e_\ell := x_\ell - x_\ast \) mit der \(\ell \)-ten Näherung \(x_\ell \) und der exakten Lösung \(x_\ast \), so gilt \(|e_{\ell +1}| \le c \cdot |e_\ell |^p\) mit \(p \ge 1\) und \(c < 1\) für \(p = 1\) (lineare Konvergenz). Die Konvergenzordnung \(p\) hat wesentlichen Einfluss auf die Konvergenzgeschwindigkeit, zum Beispiel bei \(p = 3\) (kubische Konvergenz) und \(e_0 = \frac {1}{10}\) ist schon \(|e_3| \le 10^{-14}\). Im Falle des Bisektionsverfahrens ist \(p = 1\) und \(c = \frac {1}{2}\), die Bisektion ist also recht langsam (ungefähr drei Schritte für eine Dezimalstelle).

Sekanten-Verfahren

Zwei hinreichend gute Näherungen \(x_{\ell -1}\) und \(x_\ell \) einer Nullstelle \(x_\ast \) von \(f\) können im Allgemeinen durch Bestimmung der Nullstelle der interpolierenden Gerade

\begin{align*} x_{\ell +1} = \frac {x_{\ell -1} f(x_\ell ) - x_\ell f(x_{\ell -1})} {f(x_\ell ) - f(x_{\ell -1})} = x_\ell - \frac {(x_\ell - x_{\ell -1}) f(x_\ell )} {f(x_\ell ) - f(x_{\ell -1})} \end{align*} verbessert werden.
Die wiederholte Anwendung dieser linearen Approximation heißt Sekanten-Verfahren.

Ist die approximierte Nullstelle \(x_\ast \) einfach (d. h. \(f’(x_\ast ) \not = 0\)), so besitzt das Verfahren die Konvergenzordnung \(r = \frac {1 + \sqrt {5}}{2}\) (Goldener Schnitt). Genauer gilt für den Fehler \(e_\ell := |x_\ell - x_\ast |\)

\begin{align*} \lim _{\ell \to \infty } \frac {e_{\ell +1}}{e_\ell ^r} = \left |\frac {f”(x_\ast )}{2 f’(x_\ast )}\right |^{1/r}. \end{align*} Da pro Iterationsschritt nur eine Funktionsauswertung erforderlich ist, ist das Sekanten-Verfahren etwas effizienter als die Newton-Iteration (\(r^2 > 2\)).

Inverse Interpolation

Aus Näherungen \(x_{\ell -n}, \dotsc , x_\ell \) für eine Nullstelle \(x_\ast \) einer Funktion \(f\) kann man eine Approximation \(x_{\ell +1} \approx x_\ast \) durch inverse Interpolation der Funktionswerte \(f_k = f(x_k)\) mit einem Polynom \(p\) vom Grad \(\le n\) gewinnen. Sind die Werte \(f_k\) paarweise verschieden, so ist

\begin{align*} x_{\ell +1} = p(0) = \sum _{k=\ell -n}^\ell \left (x_k \prod _{j\not =k} \frac {f_j}{f_j - f_k}\right ). \end{align*} Die Iteration des Verfahrens ist für glatte Funktionen bei einfachen Nullstellen lokal konvergent und sehr effizient. Allerdings ist der Iterationsschritt nicht immer durchführbar: Die möglichen Ausnahmefälle müssen mithilfe eines anderen Verfahrens (z. B. Bisektion) überbrückt werden.

Häufig angewendet wird die quadratische inverse Interpolation (\(n = 2\)).
Für die Daten \((x_0, f_0)\), \((x_1, f_1)\) und \((x_2, f_2)\) hat die Approximation dann die Form

\begin{align*} x_\ast \approx x_0 \frac {f_1 f_2}{(f_1 - f_0)(f_2 - f_0)} + x_1 \frac {f_0 f_2}{(f_0 - f_1)(f_2 - f_1)} + x_2 \frac {f_0 f_1}{(f_0 - f_2)(f_1 - f_2)}. \end{align*} Beginnt man mit einem Intervall \([a, b]\), an dessen Endpunkten die Funktion \(f\) ihr Vorzeichen wechselt, so kann man die inverse Interpolation sehr effektiv mit dem Bisektionsverfahren verbinden. Neben einer Folge von Approximationen \(x_0, x_1, \dotsc \) werden Intervalle \(I_\ell \) mit Vorzeichenwechsel von \(f\) gespeichert, deren einer Eckpunkt \(x_\ell \) ist. Ein Iterationsschritt \(x_\ell \rightarrow x_{\ell +1}\) verläuft wie folgt:

Ist die inverse quadratische Interpolation mit \(x_{\ell -2}\), \(x_{\ell -1}\) und \(x_\ell \) durchführbar und liefert einen Punkt \(x_{\ell +1}\) im Inneren von \(I_\ell \), so wird die Approximation akzeptiert.
Andernfalls (oder falls im letzten Schritt mit Interpolation keine signifikante Verbesserung erzielt wurde) wird \(x_{\ell +1}\) mit Bisektion aus den Endpunkten von \(I_\ell \) bestimmt.
Durch Einbeziehung des neuen Punktes wird das Intervall aktualisiert. Bei einem Bisektionsschritt werden zusätzlich \(x_{\ell -1}\) und \(x_\ell \) durch die neuen Intervallendpunkte ersetzt.

Newton-Verfahren

Mit dem Newton-Verfahren kann eine Nullstelle \(x_\ast \) einer Funktion \(f\) numerisch bestimmt werden. Die Folge \(x_0, x_1, \dotsc \) der Approximationen wird durch Linearisierung gewonnen. Die Näherung \(x_{\ell +1}\) ist der Schnittpunkt der Tangente im Punkt \((x_\ell , f(x_\ell ))\) mit der \(x\)-Achse:

\begin{align*} x_{\ell +1} = x_\ell - \frac {f(x_\ell )}{f’(x_\ell )}. \end{align*} Für eine einfache Nullstelle \(x_\ast \) konvergiert die Newton-Iteration lokal quadratisch, d. h.
\(|x_{\ell +1} - x_\ast | \le c |x_\ell - x_\ast |^2\) für Startpunkte \(x_0\) in einer hinreichend kleinen Umgebung von \(x_\ast \).

Die Voraussetzung, dass der Startwert \(x_0\) in einer hinreichend kleinen Umgebung von \(x_\ast \) liegt, ist notwendig, d. h. das Newton-Verfahren konvergiert i. A. nicht für Startwerte außerhalb der Umgebung.

Das Heron-Verfahren (auch babylonisches Wurzelziehen genannt) \(x \leftarrow (x + a/x) / 2\) zur Berechnung der Wurzel einer positiven Zahl \(a\) stellt sich bei genauerer Betrachtung als das Newton-Verfahren für \(f(x) = x^2 - a\) heraus:

\begin{align*} x \leftarrow \frac {1}{2} \left (x + \frac {a}{x}\right ) = x - \frac {x^2 - a}{2x}. \end{align*} Die Konvergenz ist z. B. für \(a = 2\), \(x_0 = 1\) äußerst schnell. Bei jedem Schritt verdoppelt sich die Anzahl der korrekten Stellen.

Färbt man bei der komplexen Newton-Iteration Bereiche der komplexen Zahlenebene je nach Konvergenzgeschwindigkeit unterschiedlich ein, so erhält man teilweise merkwürdig aussehende Fraktale (Newton-Fraktale). Ein solches ergibt sich bspw. für die Gleichung \(z^3 - 1 = 0\), \(z \in \complex \).

Müllers Verfahren

Mit Müllers Verfahren können sowohl reelle als auch komplexe Nullstellen einer Funktion \(f\) approximiert werden. Dabei wird eine Folge \(z_0, z_1, \dotsc \) von Näherungen für eine Nullstelle \(z_\ast \) mithilfe von quadratischer Interpolation generiert. Die Approximation \(z_{\ell +1}\) ist die am nächsten bei \(z_\ell \) gelegene Nullstelle der Parabel, die die Punkte \((z_\ell , f(z_\ell ))\), \((z_{\ell -1}, f(z_{\ell -1}))\), \((z_{\ell -2}, f(z_{\ell -2}))\) interpoliert:

\begin{align*} z_{\ell +1} = z_\ell - \frac {2 f(z_\ell )}{\beta _\ell \pm \sqrt {\beta _\ell ^2 - 4 f(z_\ell ) \alpha _\ell }}, \end{align*}

\begin{align*} \alpha _\ell = \Delta (z_\ell , z_{\ell -1}, z_{\ell -2}) f, \quad \beta _\ell = \Delta (z_\ell , z_{\ell -1}) f + \alpha _\ell (z_\ell - z_{\ell -1}). \end{align*} Dabei ist das Vorzeichen so gewählt, dass der Betrag des Nenners am größten wird.

Im Fall von zusammenfallenden Punkten sind die Dividierten Differenzen mithilfe entsprechender Ableitungen definiert. Allerdings steht dies nicht im Einklang mit dem ableitungsfreien Charakter des Verfahrens. In der Praxis treten jedoch solche und andere Ausnahmefälle (\(\alpha _\ell = 0\), \(z_{\ell +1} = \infty \)) sehr selten auf.

Für glatte Funktionen konvergiert Müllers Verfahren lokal fast mit Ordnung \(2\).

Schranken für Nullstellen von Polynonmen

Die Beträge der Nullstellen eines Polynoms \(p(x) = a_n x^n + \dotsb + a_1 x + a_0\) können abgeschätzt werden:

\begin{align*} |x| \le \max \left \{1, \sum _{i=0}^{n-1} \frac {|a_i|}{|a_n|}\right \}. \end{align*}

Sturmsche Kette

Die Polynomfolge \(p_n, p_{n-1}, \dotsc , p_m\) bildet eine Sturmsche Kette, falls

alle reellen Nullstellen von \(p_n\) einfach sind,
\(p_m\) sein Vorzeichen nicht ändert,
\(p_{n-1}\) an allen reellen Nullstellen von \(p_n\) ein anderes Vorzeichen als \(p_n’\) hat und
\(p_{k+1}(x) p_{k-1}(x) < 0\) für alle reellen Nullstellen \(x\) von \(p_k\), \(k = n - 1, \dotsc , m + 1\) ist.

Eine Sturmsche Kette kann zur Nullstellen-Bestimmung von Polynomen verwendet werden: Bezeichnet \(s(x)\) die Anzahl der Vorzeichenwechsel der Folge \(p_n(x), \dotsc , p_m(x)\), dann ist die Anzahl der reellen Nullstellen von \(p_n\) im Intervall \(\left [a, b\right )\) gleich der Differenz \(s(b) - s(a)\).

Diese Eigenschaft bildet die Grundlage für ein Bisektionsverfahren: Man beginnt mit einem Intervall, dass alle reelle Nullstellen von \(p_n\) enthält. Durch fortgesetzte Unterteilung, gemäß der Anzahl der Vorzeichenwechsel der Kette an den Teilintervall-Endpunkten, können so alle reellen Nullstellen von \(p_n\) bestimmt werden.

Zu einem beliebigen Polynom \(q_n\) kann man eine Sturmsche Kette \(p_n, \dotsc , p_m\) mithilfe des Euklidischen Algorithmus konstruieren. Man bestimmt dazu zunächst den größten gemeinsamen Teiler \(q_m\) von \(q_n\) und \(q_{n-1} := -q_n’\) durch sukzessive Polynomdivision

\begin{align*} q_{k+1} = r_k q_k - q_{k-1}, \quad k = n - 1, \dotsc , m \end{align*} mit \(q_{m-1} = 0\) und setzt anschließend \(p_k := q_k / q_m\).

Die charakteristischen Polynome der symmetrischen Tridiagonalmatrix

\begin{align*} \begin{pmatrix} a_1 & b_1 & & & 0 \\ b_1 & a_2 & b_2 \\ & b_2 & \ddots & \ddots \\ & & \ddots \\ & & & & b_{n-1} \\ 0 & & & b_{n-1} & a_n \end {pmatrix} \end{align*} erfüllen die Rekursion

\begin{align*} p_{n+1}(\lambda ) = (a_{n+1} - \lambda ) p_n(\lambda ) - b_n^2 p_{n-1}(\lambda ), \quad n \in \natural \end{align*} mit \(p_0(\lambda ) = 1\) und \(p_1(\lambda ) = a_1 - \lambda \). Sind alle Nebendiagonalelemente \(b_k\) ungleich Null, so sind die Nullstellen wie bei orthogonalen Polynomen geschachtelt und die Folge \(p_n, \dotsc , p_0\) bildet eine Sturmsche Kette. Mithilfe Sturmscher Ketten können somit die Eigenwerte beliebiger symmetrischer tridiagonaler Matrizen bestimmt werden.

Nullstellenbestimmung mit MATLAB

Nullstellen einer reellen Funktion f können in M ATLAB mit x = fzero(f, x0); bestimmt werden. f ist dabei als Funktionshandle, Funktionsname und Inline-Funktion gegeben. x0 ist ein Punkt, in dessen Umgebung eine Nullstelle gefunden werden soll. Statt eines Punktes kann auch ein Intervall [a, b] übergeben werden. Der verwendete Algorithmus basiert auf einer Kombination von Bisektion und inverser quadratischer Interpolation. Wird kein Intervall übergeben, so wird zunächst ein Intervall mit einem Vorzeichenwechsel der Funktion an den Endpunkten bestimmt. Kann kein solches Intervall bestimmt werden, wird ein Fehler ausgegeben und NaN zurückgegeben.

Komplexe Nullstellen können nicht gefunden werden.
Für Polynome p (\(p(z) = p_1 z^n + \dotsb + p_n z + p_{n+1}\)) schafft hier der Befehl z = roots(c); Abhilfe. Hier werden alle Nullstellen (reell wie komplex) bestimmt.

Nicht-lineare Systeme

Nicht-lineares Gleichungssystem

Ein nicht-lineares Gleichungssystem hat die Form

\begin{align*} f(x) = 0 \quad \Leftrightarrow \quad \left \{\begin{array}{rcl} f_1(x_1, \dotsc , x_n) & = & 0 \\ & \vdots \\ f_m(x_1, \dotsc , x_n) & = & 0 \end {array}\right . \end{align*} mit Unbekannten \(x_i\) und gegebenen Funktionen \(f_j\) (\(i = 1, \dotsc , n\), \(j = 1, \dotsc , m\)).

Im Gegensatz zu linearen Gleichungssystemen (LGS) können keine generellen Aussagen über die Lösbarkeit eines nicht-linearen Gleichungssystems gemacht werden. Im Allgemeinen existieren jedoch für \(m = n\) nur endlich viele Lösungen. Für \(m > n\) ist das System normalerweise überbestimmt. Es existiert keine Lösung und man spricht von einem nicht-linearen Ausgleichsproblem. Für \(m < n\) ist das System i. A. unterbestimmt, d. h. Unbekannte \(x_j\) können frei gewählt werden.

Banachscher Fixpunktsatz

Sei \(g\colon D \rightarrow \real ^n\) mit \(D \subset \real ^n\), \(D \not = \emptyset \) gegeben. Ist \(D\) abgeschlossen (\(D = \overline {D}\)), wird \(D\) von \(g\) in sich selbst abgebildet (\(g(D) \subset D\)) und ist \(g\) eine Kontraktion (\(\norm {g(x) - g(y)} \le c \norm {x - y}\) für alle \(x, y \in D\), wobei \(0 \le c < 1\)), so besitzt \(g\) einen eindeutigen Fixpunkt \(x_\ast = g(x_\ast ) \in D\).

Ausgehend von einem beliebigen Punkt \(x_0 \in D\) kann \(x_\ast \) durch die Folge

\begin{align*} x_0, \quad x_1 = g(x_0), \quad x_2 = g(x_1), \quad \dotsc \end{align*} approximiert werden. Für den Fehler gilt dabei

\begin{align*} \norm {x_\ast - x_k} \le \frac {c^k}{1 - c} \norm {x_1 - x_0}, \end{align*} d. h. die Iterationsfolge konvergiert für jeden Startwert linear.

Der Fixpunktsatz gilt auch allgemein in vollständigen metrischen Räumen. Dabei wird die Norm \(\norm {x - y}\) einfach durch die Abstandsfunktion \(d(x, y)\) ersetzt.

Zum Nachweis der Kontraktionseigenschaft von differenzierbaren Abbildungen \(g\) benutzt man oft den Mittelwertsatzes (Satz von Lagrange). Für \(n = 1\) lautet dieser

\begin{align*} g(y) - g(x) = g’(z) \cdot (y - x) \quad \text {fÃijr ein}\quad z \in \overline {x,y}. \end{align*} Lässt sich nun die Ableitung nach oben abschätzen, d. h. \(g’(z) \le \max _{z \in D} |g’(z)| = c\) mit \(c < 1\), so ist die Kontraktionseigenschaft nachgewiesen.

Für \(n \ge 2\) bedient man sich des verallgemeinerten Mittelwertsatzes:

\begin{align*} g(y) - g(x) = \int _0^1 g’(x + t(y - x))(y - x) \dt . \end{align*} Hier muss die Norm der Jacobi-Matrix nach oben abgeschätzt werden:
\(\norm {g’(x + t(y - x))} \le \max _{z \in D} \norm {g’(z)} = c\). Für \(c < 1\) liegt wieder eine Kontraktion vor.

Eine typische Anwendung ist ein leicht gestörtes System \(Ax + \varepsilon f(x) = b\) mit einer quadratischen Matrix \(A\). Die Funktion \(f(x)\) stellt die Störung dar und besitzt eine komplizierte Abhängigkeit von \(x\). Für hinreichend kleine \(\varepsilon \) dominiert aber lineare Anteil und das System kann mit der Iteration \(x \leftarrow g(x) := A^{-1} (b - \varepsilon f(x))\) gelöst werden. Man nimmt dabei an, dass \(A\) invertierbar und \(f\) Lipschitz-stetig mit Konstante \(c_f\) ist.

Zur Überprüfung der Voraussetzungen des Banachschen Fixpunktsatzes wählt man
\(D := \overline {U_r(p)} = \{y \in \real ^n \;|\; \norm {y - p} \le r\}\) mit \(p = A^{-1} b\), denn für kleines \(\varepsilon \) liegt die Lösung nahe bei der Lösung von \(Ax = b\). Für \(x \in D\) beliebig gilt dann \(\norm {g(x) - p} = \varepsilon \norm {A^{-1} f(x)} \le \varepsilon \norm {A^{-1}} \max _{y \in D} \norm {f(y)}\). Wählt man also \(\varepsilon \le \frac {r}{\norm {A^{-1}} \max _{y \in D} \norm {f(y)}}\), so gilt \(g(x) \in D\), d. h. \(g\) bildet \(D\) in sich selbst ab. Für die Kontraktionskonstante gilt \(\norm {g(x) - g(y)} = \varepsilon \norm {A^{-1} (f(x) - f(y))} \le \varepsilon \norm {A^{-1}} c_f \norm {x - y} = c \norm {x - y}\) mit \(c := \varepsilon \norm {A^{-1}} c_f\). Es gilt \(c < 1\), falls \(\varepsilon < \frac {1}{\norm {A^{-1}} c_f}\) gilt. Für ein hinreichend kleines \(\varepsilon > 0\) sind beide Bedingungen erfüllt und der Fixpunktsatz von Banach lässt sich anwenden.

Multivariates Newton-Verfahren

Ein Iterationsschritt \(x \rightarrow y\) der Newton-Iteration zur Bestimmung einer Nullstelle \(x_\ast \) eines nicht-linearen Gleichungssystems \(f_k(x_1, \dotsc , x_n) = 0\) für \(k = 1, \dotsc , n\) hat die Form

\begin{align*} \Delta := -f’(x)^{-1} f(x), \qquad y := x + \Delta . \end{align*} Dabei wird die Jacobi-Matrix \(f’(x)\) nicht invertiert, sondern das Inkrement \(\Delta \) wird als Lösung des LGS \(f’(x) \Delta = -f(x)\) bestimmt.

Für \(\det f’(x_\ast ) \not = 0\) konvergiert die Newton-Iteration lokal quadratisch, d. h.

\begin{align*} \norm {y - x_\ast } \le c \norm {x - x_\ast }^2 \quad \text {fÃijr}\quad x \approx x_\ast . \end{align*}

Kantorovich-Kriterium

Das Kantorovich-Kriterium gibt eine hinreichende Bedingung für die Konvergenz des Newton-Verfahrens für ein System nicht-linearer Gleichungen \(f_k(x_1, \dotsc , x_n) = 0\), \(k = 1, \dotsc , n\).
Gilt für einen Startvektor \(y\)

\begin{align*} \norm {f’(y)^{-1} f(y)} & \le \frac {r}{2} \\ \norm {f’(y)^{-1} (f’(x) - f’(\widetilde {x}))} & \le \frac {1}{r} \norm {x - \widetilde {x}} \end{align*} für alle \(x, \widetilde {x} \in U_r(y)\), dann existiert eine Lösung \(x_\ast \) des Systems in \(\overline {U_r(y)}\) und das Newton-Verfahren mit Startvektor \(y\) konvergiert gegen \(x_\ast \).

Fortsetzungsmethode

Bei einem von einem Parameter \(t\) abhängigen nicht-linearen Gleichungssystem
\(f_k(x_1, \dotsc , x_n, t) = 0\), \(k = 1, \dotsc , n\) kann eine Lösung \(x(t)\) für kleines \(\Delta t\) als Näherung für \(x(t + \Delta t)\) verwendet werden. Ist die Jacobi-Matrix von \(f\) bzgl. \(x\) bei \(x(t)\) invertierbar, so erhält man durch die lineare Taylor-Entwicklung

\begin{align*} x(t + \Delta t) \approx x(t) - f_x(x(t), t)^{-1} f_t(x(t), t) \Delta t \end{align*} eine verbesserte Approximation (Fortsetzungsmethode).

Die Fortsetzungsmethode kann insbesondere in Kombination mit iterativen Verfahren benutzt werden, um das nicht-lineare Gleichungssystem für eine Parameterfolge \(t_0 < t_1 < \dotsb \) zu lösen. Die Lösungen \(x(t_k)\) dienen dabei jeweils als Startwerte zur Berechnung von \(x(t_{k+1})\).

Gedämpftes Newton-Verfahren

Mit dem Newton-Verfahren wird eine Lösung \(x_\ast \) eines nicht-linearen Gleichungssystems
\(f_k(x_1, \dotsc , x_n) = 0\), \(k = 1, \dotsc , n\) ausgehend von einer hinreichend guten Startnäherung \(x\) approximiert. Beim gedämpften Newton-Verfahren will man in jedem Fall eine Verkleinerung der Norm des Funktionswerts erreichen. Daher hat ein Iterationsschritt \(x \rightarrow y\) die folgende Form:

Das Inkrement \(\Delta x\) wird durch Lösung des LGS \(f’(x) \Delta x = f(x)\) berechnet.
Man bestimmt einen Dämpfungsparameter \(\lambda \in \{1, 1/2, 1/4, \dotsc \}\), sodass \(\norm {f(y)}\) für \(y = x - \lambda \Delta x\) signifikant kleiner als \(\norm {f(x)}\) ist.

Als Test zur Bestimmung von \(\lambda \) dient der Vergleich

\begin{align*} \norm {\Delta y} \le (1 - \lambda /2) \norm {\Delta x}, \quad f’(x) \Delta y = f(y) \end{align*} mit einer geeignet gewählten Norm \(\norm {\cdot }\). Dabei kann eine bereits bestimmte LR- oder QR-Zerlegung der Jacobi-Matrix \(f’(x)\) zur schnelleren Berechnung von \(\Delta y\) benutzt werden. Durch die Multiplikation der Funktionswerte mit \(f’(x)^{-1}\) wird der Vergleich affin invariant. Insbesondere ist damit auch das gedämpfte Newton-Verfahren skalierungsunabhängig, was in einem Vergleich der Form \(\norm {f(y)} \le (1 - \rho ) \norm {f(x)}\) nicht gewährleistet wäre.

Bei der Implementierung empfiehlt es sich, \(\lambda \) nicht abrupt zu ändern. Man beginnt den Vergleich mit dem zuletzt gewählten Dämpfungsparameter. Entsprechend dem Result der Abfrage wird \(\lambda \) halbiert oder verdoppelt, wobei eine Verdoppelung dann frühestens im nächsten Iterationsschritt umgesetzt wird. Ist die Jacobi-Matrix für die approximierte Lösung \(x_\ast \) regulär, so ist \(\lambda = 1\) für Approximationen nahe genug bei \(x_\ast \). Die quadratische Konvergenz wird somit durch die Dämpfung nicht beeinträchtigt.

Gauß-Newton-Verfahren

Die Lösung \(x_\ast \in \real ^n\) eines nicht-linearen Ausgleichsproblems

\begin{align*} \norm {f(x_1, \dotsc , x_n)}_2^2 = \sum _{k=1}^m |f_k(x)|^2 \rightarrow \min \qquad (m > n) \end{align*} kann mit der durch

\begin{align*} & \norm {f(x) + f’(x) \Delta x}_2 \rightarrow \min \\ & x \leftarrow x + \Delta x \end{align*} definierten Gauß-Newton-Iteration bestimmt werden. In jedem Iterationsschritt wird dabei ein lineares Ausgleichsproblem mit der \(m \times n\)-Matrix \(f’(x)\) gelöst.

Minimierung ohne Nebenbedingungen

Goldene Suche

Ein lokales Minimum einer stetigen Funktion \(f\) kann mithilfe eines Unterteilungsalgorithmus (Goldene Suche, analog der Bisektion bei Nullstellen) bestimmt werden. Man geht dazu von drei Punkten \(a\), \(b\) und \(c\) mit

\begin{align*} f(a) \ge f(b) \le f(c) \end{align*} aus. Es muss mindestens ein lokales Minimum von \(f\) im Intervall \((a, c)\) liegen. Zur Verkleinerung des Intervalls wird \(f\) nun an einem weiterem Punkt \(x\) im größeren der Teilintervalle \((a, b)\) und \((b, c)\) ausgewertet. Dann wird einer der Eckpunkte durch \(x\) ersetzt, sodass für das neue Tripel \(\{a’, b’, c’\}\) wiederum \(f(a’) \ge f(b’) \le f(c’)\) gilt. Die Prozedur wird solange wiederholt, bis eine vorgegebene Genauigkeit erreicht ist.

Die optimale Unterteilung der Intervalle erfolgt im Verhältnis (Goldener Schnitt)

\begin{align*} r : (1 - r) \quad \text {mit}\quad r := \frac {\sqrt {5} - 1}{2} \approx 0.61803. \end{align*} Der Parameter \(r\) ist die positive Lösung der Gleichung \(r^2 = 1 - r\). Gilt \((c’ - a’) = r(c - a)\), so wird durch dieses Teilverhältnis eine konstante Reduktion der Intervalllänge pro Schritt unabhängig von dem Vergleich zwischen \(f(x)\) und \(f(b)\) erreicht.
Für beliebige Startpunkte \(a \le b \le c\) ist die Intervalllänge nach \(n\) Schritten \(\le r^{n-1} (c - a)\).

Quadratische Suche

Aus Approximationen \(x_{\ell -2}\), \(x_{\ell -1}\), \(x_\ell \) für das Minimum \(x_\ast \) einer Funktion \(f\) mit
\(\Delta (x_\ell , x_{\ell -1}, x_{\ell -2}) f > 0\) kann eine verbesserte Approximation durch Minimierung der interpolierenden Parabel bestimmt werden (quadratische Suche):

\begin{align*} x_{\ell +1} := \frac {1}{2} \left (x_\ell + x_{\ell -1} - \frac {\Delta (x_\ell , x_{\ell -1}) f} {\Delta (x_\ell , x_{\ell -1}, x_{\ell -2}) f}\right ). \end{align*} Fallen Punkte \(x_k\) zusammen, so sind die auftretenden Divideren Differenzen mithilfe der entsprechenden Ableitungswerte zu berechnen.

Ist \(f”(x_\ast ) > 0\), dann ist die quadratische Suche lokal konvergent. Insbesondere ist
\(\Delta (x_\ell , x_{\ell -1}, x_{\ell -2}) f > 0\) für \(x_k\) nahe bei \(x_\ast \), sodass die Methode immer durchführbar ist.

Steilster Abstieg

Die Methode des steilsten Abstiegs dient zur Minimierung multivariater Funktionen
\(f\colon \real ^n \rightarrow \real \). Zur Durchführung eines Iterationsschrittes wird zunächst der negative Gradient

\begin{align*} d := -\grad f(x) \end{align*} als lokal beste Abstiegsrichtung berechnet. Dann wird \(y\) als eine Minimalstelle von \(f\) in Richtung von \(d\) bestimmt:

\begin{align*} f(y) = \min _{t \ge 0} f(x + td). \end{align*} Die Suchrichtung ist dabei orthogonal zu der Niveaumenge durch \(x\) und berührt eine Niveaumenge zu einem kleineren Funktionswert in \(y\).

Die Konvergenz der durch die Methode des steilsten Abstiegs erzeugten Folge \(x_0, x_1, \dotsc \) kann unter sehr allgemeinen Voraussetzungen gezeigt werden. Hinreichend ist, dass \(f\) nach unten beschränkt ist und \(\grad f\) in einer Umgebung \(U\) der Menge \(\{x \in \real ^n \;|\; f(x) \le f(x_0)\}\) Lipschitz-stetig ist, d. h.

\begin{align*} \norm {\grad f(x) - \grad f(\widetilde {x})} \le L \norm {x - \widetilde {x}}, \quad x, \widetilde {x} \in U. \end{align*} Dann gilt

\begin{align*} \sum _{\ell =0}^\infty \norm {\grad f(x_\ell )}^2 < \infty . \end{align*} Dies impliziert insbesondere, dass jeder Häufungspunkt der Folge \(x_0, x_1, \dotsc \) ein kritischer Punkt von \(f\) ist. Dass es sich um ein lokales Minimum handelt, ist statistisch gesehen fast sicher, kann jedoch nicht zwingend gefolgert werden.

Im Algorithmus braucht die eindimensionale Minimierung nur näherungsweise durchgeführt werden. Die Suchrichtung \(d\) muss nicht als der negative Gradient gewählt und eine globale Minimalstelle \(y\) nicht bestimmt werden. Entscheidend für Konvergenz ist lediglich, dass in jedem Iterationsschritt eine Reduktion des Funktionswerts proportional zu \(\norm {\grad f(x)}^2\) erreicht wird.

Für eine symmetrische Matrix \(A\), einen Vektor \(b\) und einen Skalar \(c\) wird durch

\begin{align*} q(x) = \frac {1}{2} x^t A x + b^t x + c \end{align*} eine quadratische Funktion definiert.

Ist \(A\) symmetrisch und positiv definit, so kann man für die quadratische Funktion
\(f(x) = \frac {1}{2} x^t A x - b^t x\) ein Iterationsschritt \(x \rightarrow y\) der Methode des steilsten Abstiegs explizit angeben. Es ist \(y = x + td\) mit \(d = -\grad f(x) = b - Ax\) und \(t = \frac {d^t d}{d^t A d}\), denn das Minimum von \(f(x + td) = \frac {1}{2} (x + td)^t A (x + td) - b^t (x + td) = \frac {1}{2} d^t A d t^2 + (x^t A d - b^t d) t + c\) kann durch Nullsetzen der Ableitung nach \(t\) bestimmt werden: \(0 = d^t A d t - (b - Ax)^t d = d^t A d t - d^t d\).

Es kann zu unerwünschten Oszillationen kommen, wenn \(A\) Eigenwerte stark unterschiedlicher Größenordnung besitzt. Für

\begin{align*} A = \begin{pmatrix}1 & 0\\0 & 100\end {pmatrix}, \quad b = \begin{pmatrix}0\\0\end {pmatrix}, \quad x = \begin{pmatrix}100\\1\end {pmatrix} \end{align*} verringert sich zum Beispiel in jedem Iterationsschritt der Abstand zum Minimum im Ursprung nur um weniger als \(1\) Prozent.

Kantorovich-Ungleichung

Sei \(x \rightarrow y\) ein Schritt bei der Minimierung der quadratischen Funktion

\begin{align*} f(x) = \frac {1}{2} x^t A x - b^t x \end{align*} mit symmetrischer, positiv definiter Matrix \(A\) durch die Methode des steilsten Abstiegs.
Dann gilt die Kantorovich-Ungleichung:

\begin{align*} \norm {y - x_\ast }_A \le \frac {\kappa - 1}{\kappa + 1} \norm {x - x_\ast }_A. \end{align*} Dabei bezeichnet \(x_\ast := A^{-1} b\) die Lösung, die Kondition \(\kappa := \lambda _{\max } / \lambda _{\min }\) den Quotienten der extremalen Eigenwerte von \(A\) und \(\norm {z}_A := \sqrt {z^t A z}\) die von \(A\) induzierte Norm.

Einschub: Konjugierte Gradienten (cg-Verfahren)

\(A\)-Skalarprodukt: Zu einer symmetrischen, positiv definiten Matrix \(A\) lässt sich durch

\begin{align*} \innerproduct {x, y}_A := x^t A y \end{align*} ein Skalarprodukt definieren.
Zwei Vektoren \(x\) und \(y\) heißen \(A\)-orthogonal, falls \(\innerproduct {x, y}_A = 0\) ist.

konjugierte Gradienten: Ausgehend von einem beliebigen Startvektor \(x_0\) und
\(g_0 := -u_1 := Ax_0 - b\) erhält man mit der Iteration

\begin{align*} x_\ell & := x_{\ell -1} + \frac {\innerproduct {g_{\ell -1}, g_{\ell -1}}}{\innerproduct {u_\ell , u_\ell }_A} u_\ell \\ g_\ell & := Ax_\ell - b \\ u_{\ell +1} & := -g_\ell + \frac {\innerproduct {g_\ell , g_\ell }}{\innerproduct {g_{\ell -1}, g_{\ell -1}}} u_\ell \end{align*} die Lösung des linearen LGS \(Ax = b\) mit symmetrischer, positiv definiter \(n \times n\)-Matrix \(A\) und \(\innerproduct {\cdot , \cdot }_A\) dem \(A\)-Skalarprodukt in maximal \(n\) Schritten.
Bei exakter Rechnung ist \(g_\ell = 0\) spätestens für \(\ell = n\). Dieses Verfahren nennt man die Methode der konjugierten Gradienten (cg-Verfahren).

Mit \(f(x) = \frac {1}{2} x^t A x - b^t x\) ist \(f(x)\) minimal genau dann, wenn \(Ax = b\) (es gilt \(f’(x) = Ax - b\)), d. h. man kann das Verfahren auch als Minimierung der quadratischen Funktion \(f\) auffassen.

Für die Gradienten \(g_\ell \) und die Suchrichtungen \(u_\ell \) gilt

\begin{align*} g_\ell = g_{\ell -1} + \alpha _\ell A u_\ell \quad \text {mit}\quad \alpha _\ell := \frac {\innerproduct {g_{\ell -1}, g_{\ell -1}}}{\innerproduct {u_\ell , u_\ell }_A}, \end{align*} d. h. es ist möglich, bei der Implementierung eines Iterationsschritts
\((x_{\ell -1}, g_{\ell -1}, u_\ell ) \rightarrow (x_\ell , g_\ell , u_{\ell +1})\) mit nur einer Matrix-Multiplikation (\(Au_\ell \)) auszukommen.

Konjugierte Gradienten von Fletcher und Reeves

Das Verfahren der konjugierten Gradienten bestimmt das Minimum einer quadratischen Funktion \(f(x) = \frac {1}{2} x^t A x - b^t x\) mit einer symmetrischen, positiv definiten \(n \times n\)-Matrix \(A\) bei exakter Rechnung in höchstens \(n\) Schritten. Fletcher und Reeves formulierten den Algorithmus um, sodass dieser auf beliebige, glatte Funktion \(f\) angewendet werden kann.

Ausgehend von Startwerten

\begin{align*} x_0, \qquad g_0 := \grad f(x_0), \qquad d_0 := -g_0 \end{align*} erzeugt man eine Folge von Näherungen \(x_\ell \) für eine Minimalstelle und Suchrichtungen \(d_\ell \) durch folgende Rekursionen:

\begin{align*} x_{\ell +1} & := x_\ell + \alpha _\ell d_\ell \\ g_{\ell +1} & := \grad f(x_{\ell +1}) \\ d_{\ell +1} & := -g_{\ell +1} + \beta _\ell d_\ell , \qquad \beta _\ell := \frac {\innerproduct {g_{\ell +1}, g_{\ell +1}}}{\innerproduct {g_\ell , g_\ell }}, \end{align*} wobei \(\alpha _\ell > 0\) bestimmt ist durch Minimierung von \(f(x_\ell + \alpha d_\ell )\) für \(\alpha > 0\).

Der einzige Unterschied zum quadratischen Fall ist, dass \(\alpha _\ell \) nicht explizit bestimmt werden kann.

Eine gute Performance kann besonders dann erzielt werden, wenn \(f\) gut durch eine quadratische konvexe Funktion approximiert wird. Ist dies nicht der Fall, so sollte in geeigneten Abständen ein Neustart des Verfahrens erfolgen. Die Konvergenzgeschwindigkeit steigt in der Nähe des Minimums rapide an (in der Nähe des Minimums ähnelt jede Funktion stark einer quadratischen Funktion).

Die eindimensionale Minimierung wird i. A. nicht exakt durchgeführt. Dann ist jedoch darauf zu achten, dass

\begin{align*} \innerproduct {g_{\ell +1}, d_{\ell +1}} < 0, \end{align*} d. h. die Suchrichtungen sind lokale Abstiegsrichtungen. Ist \(x_{\ell +1}\) ein lokales Minimum von \(f\) in Richtung \(d_\ell \), so gilt \(\innerproduct {g_{\ell +1}, d_\ell } = 0\), sodass aufgrund der Definition von \(d_{\ell +1}\) diese Bedingung automatisch erfüllt ist.

Es existieren einige Varianten bei der Parameterwahl, die ebenfalls mit dem quadratischen Fall konsistent sind. Beispielsweise definieren Polak und Ribiere

\begin{align*} \beta _\ell := \frac {\innerproduct {g_{\ell +1} - g_\ell , d_\ell }}{\innerproduct {g_\ell , d_\ell }}. \end{align*} Diese Wahl führt in der Praxis oft zu besseren Ergebnissen als die klassische Variante von Fletcher und Reeves.

Minimierung mit MATLAB

Ein lokales Minimum einer reellen Funktion auf einem Intervall \([a, b]\) kann in M ATLAB mit dem Befehl [x, fx] = fminbnd(f, a, b); bestimmt werden. Die Funktion f wird als Funktionshandle oder Inline-Funktion übergeben. Der Rückgabewert x enthält die gefundene Minimalstelle und der optionale Rückgabewert fx den entsprechenden Funktionswert. Es werden sowohl lokale Randminima als auch innere lokale Minima gefunden, jedoch nicht immer das globale Minimum.

Zur Minimierung multivariater Funktionen steht [x, fx] = fminsearch(f, x0); zur Verfügung. Damit wird ein lokales Minimum in der Nähe eines Startvektors x0 gefunden.