3.2. Unrestringierte Optimierung#
Im Rahmen dieser Vorlesung wollen wir uns zunächst auf eine bestimmte Klasse von allgemeinen Optimierungsproblemen konzentrieren, den unbeschränkten oder unrestringierten Optimierungsproblemen.
Definition 3.4 (Unrestringierte Optimierung)
Liegt ein allgemeines Optimierungsproblem der Form (3.1) ohne Nebenbedingungen vor, d.h.,
so sprechen wir von einem unrestringierten Optimierungsproblem.
3.2.1. Optimalitätsbedingungen#
Wir werden uns später intensiv mit der Bestimmung und numerischen Approximation von Optima beschäftigen. Dazu sollten wir uns zunächst einmal Gedanken darüber machen, wie wir das anstellen wollen, denn die Bedingungen aus Definition 3.3 sind keineswegs praxistauglich. Wir formulieren deshalb Optimalitätsbedingungen. Diese erlauben uns nicht nur, die Optimalität eines Punktes numerisch günstig nachzuweisen, sondern geben uns auch eine klare Vorstellung, wonach wir eigentlich suchen müssen.
Satz 3.1 (Notwendige Optimalitätsbedingungen 1. Ordnung)
Sei \({f:\R^n\to\R}\) auf der offenen Menge \(U\subset\R^n\) stetig differenzierbar. Ist \({x^\ast\in U}\) ein lokales Minimum von \(f\), dann gilt
Proof. Wir führen einen Beweis durch Widerspruch. Nehmen wir also an, dass \({x^* \in \mathbb{R}}\) ein lokaler Minimierer von \(f\) sei, aber \({\nabla f(x^*) \neq 0}\) gelte. Wir wählen den Richtungsvektor \(\vec{p} \in \R^n\) mit \(\vec{p} \coloneqq -\nabla f(x^*) \neq 0\). Es ist somit klar, dass
Da \(\nabla f\) nach Voraussetzung stetig in einer lokalen Umgebung \(U \subset \R^n\) von \(x^*\) ist, existiert ein \(T > 0\), so dass auch gilt:
Nach dem Satz von Taylor gilt aber auch für jedes \(\tilde{t} \in (0,T]\):
Somit gilt also \(f(x^* + \tilde{t}\vec{p}) < f(x^*)\) für alle \(\tilde{t} \in (0, T]\) und wir haben offenbar eine Richtung \(\vec{p} \in \mathbb{R}^n \setminus \lbrace 0\rbrace\) gefunden in der die Funktionswerte von \(f\) abnehmen. Also ist \(x^*\) kein lokaler Minimierer von \(f\).
Das ist aber ein Widerspruch zur Annahme und somit ist die Behauptung bewiesen. ◻
Die für die Optimierung interessanten Punkte \(x^* \in \Omega\), die die notwendigen Optimalitätsbedingungen aus Satz 3.1 erfüllen, nennen wir stationäre Punkte.
Definition 3.5 (Stationärer Punkt)
Sei \({f:\R^n\to\R}\) auf der offenen Menge \(U\subset\R^n\) stetig differenzierbar. Dann heißt \({x^\ast\in U}\) mit \({\nabla f(x^\ast)=0}\) stationärer Punkt von \(f\).
Mit der Definition von stationären Punkten lässt sich folgendes Korollar direkt ableiten:
Korollar 3.1
Jeder lokale Minimierer \(x^* \in \R^n\) einer stetig differenzierbaren Zielfunktion \(f\) ist ein stationärer Punkt.
Bemerkung 3.3 (Sattelpunkte)
Die Umkehrung der Aussage in Satz 3.1 gilt im Allgemeinen nicht. Man betrachte zum Beispiel die Zielfunktion \(f(x) \coloneqq x^3\). Diese besitzt einen stationären Punkt in \(x^* = 0\), d.h., es gilt \(\nabla f(0) = 0\). Dennoch handelt es sich hierbei nicht um einen lokales Minimierer, sondern lediglich um einen Sattelpunkt. Aus diesem Grund handelt es sich nur um eine notwendige, nicht um eine hinreichende Bedingungen.
Bei der Suche nach lokalen Minima lässt sich ein weiteres Kriterium anwenden, welches die zweite Ableitung der Funktion verwendet.
Satz 3.2 (Notwendige Optimalitätsbedingungen 2. Ordnung)
Sei \(f:\R^n\to\R\) auf der offenen Menge \(U\subset\R^n\) zweimal stetig differenzierbar. Ist \(x^\ast\in U\) ein lokales Minimum von \(f\), dann gilt:
\(\nabla f(x^\ast) = 0\),
Die Hessematrix \(\nabla^2 f(x^\ast)\) ist positiv semidefinit, d.h.
\[\langle \vec{p}, \nabla^2f(x^*) \vec{p} \rangle \, \geq \, 0 \qquad \forall \, \vec{p} \in \mathbb{R}^n.\]
Proof. In den Übungsaufgaben zu zeigen. ◻
Schlussendlich wollen wir auch eine hinreichende Bedingung für das Vorliegen eines lokalen Minimums angeben.
Satz 3.3 (Hinreichende Optimalitätsbedingungen 2. Ordnung)
Sei \(f:\R^n\to\R\) auf der offenen Menge \(U\subset\R^n\) zweimal stetig differenzierbar. Sei \({x^\ast\in U}\) ein Punkt, in dem gilt:
\(\nabla f(x^\ast)=0\),
Die Hessematrix \(\nabla^2 f(x^\ast)\) ist positiv definit, d.h.
\[\langle \vec{p}, \nabla^2f(x^*) \vec{p} \rangle \, > \, 0 \qquad \forall \, \vec{p} \in \mathbb{R}^n\setminus\{0\}.\]
Dann ist \(x^\ast\) ein striktes lokales Minimum von \(f\).
Proof. Da die Hessematrix \(\nabla^2 f\) nach Voraussetzung stetig und im Punkt \(x^* \in U\) positiv definit ist, können wir einen Radius \(r > 0\) finden, so dass \(\nabla^2 f(x)\) positiv definit ist für alle \(x \in B_r(x^*)\). Für jeden Vektor \(\vec{p} \in \mathbb{R}^n \setminus \lbrace 0\rbrace\) mit \(\Vert\vec{p}\Vert < r\) gilt dann nach dem Satz von Taylor:
Da \(\Vert t\vec{p}\Vert < r\) ist nach Konstruktion wissen wir, dass
gilt und somit schon \(f(x^* + \vec{p}) > f(x^*)\) gelten muss. Da \(\vec{p} \in \mathbb{R}^n \setminus \lbrace 0 \rbrace\) mit \(\Vert\vec{p}|\Vert < r\) beliebig gewählt war handelt es sich bei \(f(x^*)\) um ein striktes lokales Minimum der Zielfunktion \(f\). ◻
Bemerkung 3.4 (Definitheit der Hessematrix)
Die in Satz 3.3 genannten Bedingungen sind hinreichend, jedoch nicht notwendig für das Vorliegen eines strikten lokalen Minimums. Dies sieht man ein, wenn man beispielsweise die Zielfunktion \(f(x) \coloneqq x^4\) betrachtet. \(f\) besitzt ein striktes lokales Minimum in \(x^*=0\) und es gilt \(\nabla f(0) = 0\), jedoch verschwindet die zweite Ableitung \(\nabla^2f(0) = 0\) und ist somit nicht positiv definit.
Eine äußerst wertvolle Eigenschaft bei der Optimierung ist die Konvexität einer Zielfunktion, da jedes lokale Optimum einer konvexen Funktion bereits ein globales Optimum ist.
Definition 3.6 (Konvexität)
Sei \(\Omega \subset \mathbb{R}^n\) eine konvexe Menge und sei \({f:\Omega\to\R}\). Wir nennen \(f\) konvex, wenn für beliebige Vektoren \({x,y \in \Omega}\) die folgende Ungleichung für alle \({0 \leq \alpha \leq 1}\) gilt:
Anschaulich bedeutet Konvexität einer Funktion \(f\), dass jede Verbindungsgerade zwischen zwei Punkten \(f(x)\) und \(f(y)\) in diesem Intervall oberhalb des Graphen der Funktion \(f\) verläuft.
3.2.2. Gradientenabstiegsverfahren#
Wir wissen aus Satz 3.1, dass wir uns bei der Suche nach (lokalen) Minimierern auf die Suche nach stationären Punkten beschränken können. Leider können wir die aus der Schule bekannte Strategie Ableiten – Nullsetzen – Typ untersuchen in den meisten Fällen nicht direkt anwenden. Ein wesentlicher Grund hierfür ist, dass \({\nabla f}\) in vielen Anwendungsproblemen keine explizite Darstellung hat, sondern (z.B. durch teure numerische Simulationen und finite Differenzen) approximiert werden muss. In anderen Worten: Wir können zwar für gegebenes \({x\in\R^n}\) den Gradienten \({\nabla f(x)}\) bestimmen, den funktionalen Zusammenhang \({x\mapsto\nabla f(x)}\) haben wir aber schlicht nicht zur Verfügung.
Deshalb orientieren wir uns an den iterativen Fixpunktverfahren, die wir in der Numerik kennengelernt haben. Dazu stellen wir fest, dass
Wir können stationäre Punkte also durch das Fixpunktverfahren
bestimmen.
Beispiel 3.2
Wir betrachten die Funktion \({f(x):=x^2}\) mit (eindeutigem) stationären Punkt \({x^\ast=0}\). Für beliebiges \({x_0\neq0}\) gilt mit der Verfahrensfunktion \({G(x):=x-\nabla f(x)}\) allerdings
Das Verfahren ist also für [kein]{.underline} \({x_0\neq0}\) konvergent!
Wie wir in Beispiel 3.2 sehen, besteht für diese Form des Gradientenabstiegsverfahrens die Gefahr, einen stationären Punkt (also ein potentielles Minimum) zu überspringen. In der Tat bleibt im obigen Beispiel der Zielfunktionswert über alle Iterationen unverändert. Andererseits sehen wir mithilfe der zugehörigen Taylorapproximation
dass wir für \(\nabla f(x_k)\neq0\) immer in der Lage sind, eine Verkleinerung des Zielfunktionswerts zu erreichen, sofern wir \(\tau>0\) klein genug wählen. Das Gradientenabstiegsverfahren (im Englischen: gradient descent (GD)) ist somit zwangsläufig mit einer Schrittweite \(\tau_k>0\) verknüpft, die es passend zu bestimmen gilt. Dies lässt sich durch folgenden Algorithmus umsetzen:
Algorithmus 3.1 (Gradientenabstiegsverfahren)
function \([x^*, f(x^*)]=\,\)gradientDescent\((f,\nabla f, x_0)\) # Wahl der Schrittweite \(\tau = \tau_k\) # Update \(x_{k+1} = x_k - \tau\nabla f(x_k)\) # Ausgabe des letzten Punktes \(x^* = x_k\) \(f(x^*) = f(x_k)\)
Bevor wir auf die Wahl der Schrittweite in Algorithmus 3.1 näher eingehen wollen, stellt sich uns eine andere offene Frage: Wann beenden wir das Verfahren? Als sog. Abbruchkriterien kommen in der Praxis viele verschiedene Kandidaten in Frage. Einige geeignete Möglichkeiten wären z.B.:
Abbruch nach einer festen Anzahl an Iterationen.
Abbruch, wenn wir einen fast stationären Punkt erreicht haben, d.h. sobald \({\big\Vert\nabla f(x_k)\big\Vert<\texttt{tol}}\).
Abbruch, wenn die Zielfunktionswerte nicht weiter abnehme, d.h. \(f(x_{k+1})\ge f(x_k)\).
3.2.2.1. Konstante Schrittweiten#
So simpel die Kernidee hinter dem Gradientenabstiegsverfahren (Algorithmus 3.1) auch sein mag, die Wahl einer passenden Schrittweitenfolge \((\tau_k)_{k\in\N}\) ist keine einfache Aufgabe. Gängige moderne Gradientenverfahren nutzen dazu meinst eine sogenannte Liniensuche, bei der die Informationen über \(f(x_k)\) und \(\nabla f(x_k)\) verwendet werden, um eine gute Schätzung für die Schrittweite \(\tau_k\) zu erhalten. Diese wird dann getestet und unter Umständen angepasst. Wir überlassen solche adaptiven Verfahren anderen Vorlesungen und konzentrieren uns hier erstmal auf den einfachsten Fall: konstante Schrittweiten \(\tau_k=\tau>0\). Wie wir in Beispiel 3.2 bereits gesehen haben, können wir \(\tau\) dabei nicht beliebig groß wählen. Zunächst definieren wir also eine Klasse von angenehmen Funktionen.
Definition 3.7
Eine differenzierbare Funktion \(f:M\subseteq\R^n\to\R\) heißt \(L\)-glatt, wenn eine Konstante \(L\ge0\) existiert mit
Für \(L\)-glatte Funktionen lässt sich die folgende wichtige Eigenschaft zeigen:
Lemma 3.1 (Abstiegslemma)
Sei \(f:M\subseteq\R^n\to\R\) eine \(L\)-glatte Funktion. Dann gilt für alle \(x,y\in\operatorname{int}(M)\)
Proof. Sei \(h:=y-x\) und
Dann ist \(\phi\) differenzierbar und es gilt
wobei wir im letzten Schritt ausgenutzt haben, dass \(f\) \(L\)-glatt ist. Aus dem Hauptsatz der Integral- und Differentialrechnung folgt somit
Sortieren wir diese Ungleichung um, erhalten wir die Aussage. ◻
Mithilfe von Lemma 3.1 können wir die Konvergenz des Gradientenverfahrens mit konstanten Schrittweiten beweisen.
Satz 3.4 (Globale Konvergenz)
Sei \(f:\R^n\to\R\) von unten beschränkt und \(L\)-glatt. Außerdem sei die konstante Schrittweite \(\tau\) des Gradientenabstiegsverfahrens so gewählt, dass
Dann gilt
Proof. Mit \(x = x_k\) und \(y = x_{k+1} = x_k-\tau\nabla f(x_k)\) folgt mithilfe von Lemma 3.1
also
Wegen \(0<\tau<\tfrac{2}{L}\) ist \(\tfrac{L\tau^2}{2}-\tau < 0\). Die Folge \((f(x_k))_{k\in\N}\) ist also monoton fallend. Da \(f\) nach unserer Annahme von unten beschränkt ist, ist die Folge der Funktionswerte sogar konvergent. Demnach existiert \(f^\ast\in\R\) mit
Leider besitzen wir keine direkte Möglichkeit, Funktionswerte und Gradienten miteinander in Verbindung zu bringen. Allerdings kennen wir einen Zusammenhang zwischen Gradienten und Differenzen von Funktionswerten (3.3). Deswegen benutzen wir den üblichen Teleskopsummentrick und schreiben \(f(x_k)\) für \(k\ge 1\) wie folgt um:
Dadurch erhalten wir
Wegen \(\tfrac{L\tau^2}{2}-\tau < 0\) muss die auftretende Reihe konvergieren. Insbesondere ist \(\Vert \nabla f(x_k)\Vert\) eine Nullfolge, woraus die globale Konvergenz des Gradientenverfahrens folgt. ◻
3.2.2.2. Gradient flow#
Betrachten wir erneut die Verfahrensvorschrift für das Gradientenverfahren
so erkennen wir eine deutliche Ähnlichkeit zu den Einschrittverfahren für Anfangswertprobleme, die wir im ersten Abschnitt der Vorlesung kennengelernt haben. Genauer gesagt lässt sich (3.4) als Vorwärts-Eulerverfahren, angewendet auf das Anfangswertproblem
interpretieren. Dieses Anfangswertproblem, auch bekannt als Gradientenfluss (engl. gradient flow) können wir natürlich auch mit anderen Diskretisierungsverfahren lösen.
Von besonderem Interesse ist hierbei das Verfahren, das aus dem Rückwärts-Euler-verfahren resultiert, weil es auch Anwendung in der nichtglatten Optimierung findet, wie wir später feststellen werden. Dazu nehmen wir an, dass \(f:\R^n\to\R\) zusätzlich konvex ist. Das Rückwärts-Eulerverfahren, angewendet auf den Gradientenfluss, hat die Form
Um diese spezielle implizite Gleichung zu lösen, stellen wir fest, dass sich (3.5) als Optimalitätsbedingung für die Minimierung der (strikt konvexen!) Funktion
interpretieren lässt, denn es gilt
Dies motiviert unsere folgende Definition:
Definition 3.8 (Proximaloperator für differenzierbare Funktionen)
Für eine konvexe, differenzierbare Funktion \(f:\R^n\to\R\) definieren wir den Proximaloperator
Mithilfe des Proximaloperators können wir ein simples proximales Gradientenverfahren konstruieren:
Algorithmus 3.2 (Proximales Gradientenverfahren)
function \([x^*, f(x^*)]=\,\)ProxGradientDescent\((f, x_0)\) # Wahl der Schrittweite \(\tau = \tau_k\) # Update \(x_{k+1} = \operatorname{prox}_{\tau f}(x_k)\) # Ausgabe des letzten Punktes \(x^* = x_k\) \(f(x^*) = f(x_k)\)
Wir werden später sehen, dass sich der Proximaloperator auch auf nicht-differenzierbare Funktionen erweitern lässt.
3.2.3. Stochastisches Gradientenabstiegsverfahren#
Eine aktuell weit verbreitete Variante des Gradientenabstiegsverfahrens in (3.4) ist das stochastische Gradientenverfahren (im Englischen: stochastic gradient descent (SGD)). Dieser Name ist leider etwas irreführend, denn es handelt sich nicht um ein Abstiegsverfahren, d.h. die Zielfunktionswerte müssen während des Verfahrens nicht zwangsläufig abnehmen.
Wie der Name schon verrät handelt es sich hierbei nicht um einen deterministischen Algorithmus. Das bedeutet, dass man bei mehrmaliger Anwendung des Verfahrens bei gleichbleibenden Startbedingungen in der Regel unterschiedliche Ergebnisse in unterschiedlichen Laufzeiten erhält. Was auf den ersten Blick wie ein Nachteil wirkt, kann in manchen Fällen jedoch praktische Eigenschaften mit sich bringen. So kann die Zufallsnatur des stochastischen Gradientenverfahrens dazu führen, dass Sattelpunkte und ungewollte, lokale Minima der Funktion durch die Folge der Punkte vermieden werden. Das Verfahren findet aktuell vor allem beim Training von neuronalen Netzen bei der sogenannten Backpropagation in verschiedenen Variationen Anwendung, da man hierdurch dem bekannten Problem des Übertrainierens des neuronalen Netzes entgegen wirken kann. Außerdem lässt sich zeigen, dass es im richtigen Kontext deutlich effizienter als ein deterministisches Gradientenabstiegsverfahren ist.
Beim stochastischen Gradientenverfahren geht man davon aus, dass sich die zu minimierende Zielfunktion \(f:\R^n\to\R\) als Erwartungswert einer Größe schreiben lässt, d.h.
wobei \(\xi\) ein \(\R^d\)-wertiger Zufallsvektor mit Wahrscheinlichkeitsverteilung \(\mathbb{P}\) ist, der Werte in \(\Xi\subseteq\R^d\) annimmt. Ein einfacher (aber häufig auftretender) Spezialfall ist, dass sich \(f\) als endliche Summe
schreiben lässt.
Die Idee des stochastischen Gradientenverfahrens ist es nun, die Berechnung des vollständigen Gradienten
zu vermeiden, weil dies i.A. zu aufwendig ist. Dazu ziehen wir in jeder Iteration zufällig eine Realisierung \(\xi_k\sim\mathbb{P}\) und verwenden den zugehörigen stochastischen Gradienten \(\nabla_x F(x_k,\xi_k)\) als Richtung für den aktuell Schritt (im Falle einer endlichen Summe ziehen wir für (3.6) zufällig einen Index \(i_k\in\{1,\ldots,m\}\) und betrachten \(\nabla F_{i_k}(x_k)\)). Das stochastische Gradientenverfahren hat also die folgende Form:
Algorithmus 3.3 (Stochastisches Gradientenverfahren)
function \(x^\ast=\,\)StochasticGradient\((F, x_0)\) # Ziehe zufällige Realisierung \(\xi_k\sim\mathbb{P}\) # Wahl der Schrittweite \(\tau = \tau_k\) # Update \(x_{k+1} = x_k - \tau\nabla_x F(x_k,\xi_k)\) # Ausgabe des letzten Punktes \(x^* = x_k\)
Das stochastische Gradientenverfahren konvergiert i.A. nicht, wenn konstante Schrittweiten verwendet werden. Stattdessen wählt man hier für gewöhnlich Schrittweitenfolgen \((\tau_k)_{k\in\N}\) mit der Eigenschaft
also z.B. \(\tau_k=\tfrac{c}{k}\). Ein übliches Konvergenzresultat (auf dessen Beweis wir an dieser Stelle verzichten), lautet z.B.:
Satz 3.5
Seien die Realisierungen \(\xi_k\) unabhängig und identisch verteilt. Wir nehmen weiterhin an, dass die Zielfunktion gleichmäßig konvex ist, d.h. es existiert \(\mu>0\) mit
Sei \(M>0\) mit
Dann produziert das stochastische Gradientenverfahren (Algorithmus 3.3) mit Schrittweite \(\tau_k=\theta/k\), mit Konstante \(\theta > 1/(2\mu)\), eine Folge von Iterierten \((x_k)_{k\in\N}\) mit der Eigenschaft
wobei
Bemerkung 3.5
Das Konvergenzresultat für SG ist deutlich schwächer als unser Konvergenzresultat für das gewöhnliche Gradientenabstiegsverfahren (Satz 3.4). Allerdings müssen wir bedenken, dass die Berechnung eines vollständigen Gradienten eine Integration über \(\Xi\subseteq\R^d\) erfordert. Zwar können wir dieses Integral mit denen uns aus der Numerik bekannten Quadraturregeln numerisch approximieren, aber der Aufwand hierfür steigt exponentiell in \(\operatorname{dim}(\Xi)\) an. Im Gegensatz dazu liefert uns Satz 3.5 eine dimensionsunabhängige Konvergenzrate. Der Rechenaufwand des stochastischen Gradientenverfahrens ist ebenfalls dimensionsunabhängig, da ausschließlich stochastische Gradienten (anstelle des vollständigen Gradienten) bestimmt werden.
3.2.4. Newton Verfahren#
Bisher haben wir lediglich Verfahren untersucht, die ausschließlich Gradienten- und Zielfunktionsinformationen für das Iteriertenupdate verwenden (sog. Verfahren 1. Ordnung). Jetzt stellt sich natürlich die Frage, ob wir stärkere Konvergenzresultate erhalten, wenn wir höhere Ableitungen von \(f\) in unser Verfahren mit einfließen lassen.
Dazu wollen wir uns das bereits bekannte Newton-Verfahren in Erinnerung rufen und dieses geeignet zur Optimierung von nichtlinearen Funktionen verallgemeinern. In [TR] haben wir das Newton Verfahren zur Approximation von Nullstellen nichtlinearer Gleichungssysteme hergeleitet. Wir haben zunächst die Taylorapproximation einer nichtlinearen Nullstellengleichung \(f(x^*) = 0\) von der folgenden Form betrachtet
wobei \(f'\) die als regulär angenommene Jacobi-Matrix der differenzierbaren Funktion \(f \colon \R^n \rightarrow \R^n\) bezeichnet. Hierauf basierend haben wir die folgende Fixpunktfunktion als Approximation erster Ordnung angegeben:
Hierbei haben wir die Fixpunktgleichung als erfüllt gesehen, wenn wir ein \(x^* \in \Omega\) gefunden haben, so dass für die Fixpunktgleichung (3.7) gilt \(x^* = G(x^*)\). Unter dieser Beobachtung haben wir das Newton-Verfahren als iteratives Schema zur Bestimmung eines solchen Fixpunktes \(x^* \in \Omega\) hergeleitet:
Hierfür benötigten wir einen geeigneten Startwert \(x_0 \in \Omega\) in einer lokalen Umgebung \(U \subset \Omega\) des Fixpunktes \(x^* \in U\).
Der folgende Satz formuliert Bedingungen für die lokale Konvergenz des Newton-Verfahrens.
Satz 3.6 (Lokale Konvergenz des Newton Verfahrens)
Sei \(f: \R^n \rightarrow \R^n\) in einer Umgebung von \(\overline{x} \in \R^n\) stetig differenzierbar und \(\overline{x}\) sei eine Nullstelle von \(f\) mit \(f(\overline{x}) = 0\). Sei außerdem die Jacobi-Matrix \(f'\) lokal Lipschitz-stetig und \(f'(\overline{x})\) regulär in der Nullstelle.
Dann existiert eine lokale Umgebung \(B_R(\overline{x})\), so dass das Newton-Verfahren für jeden Startwert \(x_0 \in B_R(\overline{x})\) gegen die Nullstelle \(\overline{x}\) konvergiert, d.h. es gilt \(\lim_{x\rightarrow \infty} x_k = \overline{x}\).
Proof. Siehe [TR]. ◻
Anstatt nun eine Nullstelle der Funktion \(f\) zu suchen, wollen wir das Newton-Verfahren nutzen, um eine Nullstelle des Gradienten \(\nabla f\) (d.h einen stationären Punkt von \(f\)) zu approximieren und damit die notwendigen Optimalitätsbedingungen in Satz 3.1 zu erfüllen. Im Folgenden sei \(\Omega \subset \mathbb{R}^n\) ein offenes, zusammenhängendes Gebiet und \(f \colon \Omega \rightarrow \mathbb{R}\) eine differenzierbare, reellwertige Funktion. Wir betrachten wieder die Taylorapproximation der Funktion \(f\) in eine Richtung \(x_k + p \in \Omega\), aber berücksichtigen diesmal auch Terme von zweiter Ordnung:
Unter gewissen Bedingungen an die Hessematrix \(\nabla^2 f(x_k)\), lässt sich ein eindeutiges Minimum der Modellfunktion \(m_k(p)\) in (3.8) bestimmen, wie folgendes Theorem besagt.
Satz 3.7 (Newton-Abstiegsrichtung)
Sei \(\Omega \subset \R^n\) ein offenes, zusammenhängendes Gebiet. Sei außerdem \(f \colon \Omega \rightarrow \R\) eine in einer lokalen Umgebung eines Punktes \(x_k \in \Omega\) zweimal stetig differenzierbare Zielfunktion, deren Hessematrix \(\nabla^2 f(x_k)\) im Punkt \(x_k\) positiv definit ist.
Dann ist der Newton-Abstiegsrichtung benannte Vektor \(p_k^N \in \R^n\) mit
das eindeutige Minimum der Modellfunktion \(m_k(p)\) in (3.8).
Proof. In den Übungsaufgaben zu zeigen. ◻
Mit der Newton-Abstiegsrichtung in (3.9) lässt sich ein iteratives Abstiegsverfahren für einen initialen Punkt \(x_0 \in \Omega\), welcher geeignet in der Nähe des stationären Punktes \(x^* \in \Omega\) gewählt wird, wie folgt konstruieren:
Damit das Newton-Abstiegsverfahren in (3.10) überhaupt sinnvoll ist, müssen wir fordern, dass die Hessematrix in jedem Punkt \(x_k \in \Omega\) der Iterationsfolge regulär und somit invertierbar ist. Um sicher zu gehen, dass es sich tatsächlich um eine Abstiegsrichtung handelt müssen wir fordern, dass die Hessematrix \(\nabla^2 f(x_k)\) nicht nur invertierbar für alle \(x_k \in \Omega\) der Iterationsfolge ist, sondern auch positiv definit in jedem Punkt \(x_k\) ist. Denn dann ergibt eine Taylorapproximation zweiter Ordnung die folgende Abschätzung:
Wir sehen also, dass wir einen echten Abstieg der Funktionswerte erhalten, wenn die Hessematrix \(\nabla^2 f(x_k)\) positiv definit ist für alle \(x_k \in \Omega\) der Iterationsfolge. Sollte die Hessematrix nicht positiv definit in einem Punkt \(x_k\) der Iterationsfolge sein, so muss zumindest eine Abnahme der Funktionswerte vorliegen, d.h., es muss für die Newton-Abstiegsrichtung gelten:
Sollte dies nicht der Fall sein, so existieren Methoden um dennoch einen Abstieg zu erzwingen, siehe zum Beispiel [NW99]. Auf diese werden wir jedoch im weiteren Verlauf der Vorlesung nicht näher eingehen.
Bemerkung 3.6 (Schrittweite und Konvergenz)
Das Newton-Abstiegsverfahren in (3.10) ist ein Abstiegsverfahren, dessen Schrittweite \(\tau_k > 0\) implizit durch die lokale Krümmung und die Ableitung der Funktion \(f\) bestimmt ist. In diesem Fall können wir \(\tau_k \equiv 1\) für alle \(k \in \mathbb{N}\) setzen.
Aus der Numerik wissen wir, dass das Newton-Abstiegsverfahren in der Regel quadratisch gegen einen stationären Punkt \(x^* \in \Omega\) mit \(\nabla f(x^*)=0\) konvergiert, d.h. man erreicht sehr schnell eine hohe Genauigkeit bei der Approximation von \(x^*\).
3.2.5. Quasi-Newton Verfahren#
Im Newton Verfahren haben wir das Newton Verfahren zur iterativen Approximation eines stationären Punktes \(x^* \in \Omega\) einer Funktion \(f\) mit \(\nabla f(x^*) = 0\) hergeleitet. Hierbei haben wir im Gegensatz zu den vorherigen numerischen Verfahren auch Ableitungen höherer Ordnung hinzugezogen. Dies führt in der Regel zu einem verbesserten Konvergenzverhalten im Vergleich zu den Verfahren, die nur die lokale Ableitung \(\nabla f\) der Zielfunktion \(f\) verwenden.
Dennoch ist das Newton Verfahren aus numerischer Sicht noch nicht ideal, da es einige Probleme mit sich bringt. Zuerst mussten wir fordern, dass die Hessematrix \(\nabla^2 f(x_k)\) in jedem Punkt des Iterationsverfahrens positiv definit ist, da ansonsten kein Abstieg der Funktionswerte garantiert werden kann. Zweitens muss für die Berechnung der Newton-Richtung in (3.9) zuerst die Hessematrix bestimmt und anschließend invertiert werden. Dies ist aus Effizienzgründen unerwünscht, da die Inversion einer \(n \times n\)-Matrix bekanntlich in \(\mathcal{O}(n^3)\) Rechenoperationen liegt (siehe [TR]. In vielen Fällen ist sogar allein die Berechnung von \(\nabla^2 f\) schon zu zeitaufwendig. Da die Bestimmung und die Inversion der Hessematrix in jedem Iterationsschritt passieren müssen, ist das Newton Verfahren nur eingeschränkt empfehlenswert für die numerische Optimierung.
Eine naheliegende Idee ist es nun die echte Hessematrix in jedem Iterationsschritt durch eine geeignete Matrix zu approximieren, so dass der numerische Aufwand geringer wird, d.h., wir suchen nach einer Matrix \(B_k \in \R^{n \times n}\)
Damit können wir die Modellfunktion \(m_k(p)\) in (3.8) schreiben als:
das heißt, wir approximieren die Zielfunktion \(f\) im \(k\)-ten Iterationsschritt entlang der Richtung \(p \in \mathbb{R}^n\) lokal durch eine quadratische Funktion. Für sehr kleine Schrittweiten können wir davon ausgehen, dass der Fehler dieser Approximation gering ist, da wir davon ausgehen, dass \(f\) stetig differenzierbar in einer lokalen Umgebung \(U \subset \Omega\) des stationären Punktes \(x^* \in \Omega\) ist und für \(p = 0\) die Approximation exakt ist, da gilt
Wenn wir fordern, dass \(B_k\) in (3.11) eine positiv definite Matrix ist, so lässt sich ein Abstiegsschritt analog zur Herleitung des Newton Abstiegsverfahrens in Newton Verfahren angeben als:
Die sogenannten Quasi-Newton Verfahren verfolgen diesen Ansatz.
Bemerkung 3.7 (Konvergenzgeschwindigkeit Quasi-Newton Verfahren)
Durch die Approximation der echten Hessematrix verlieren Quasi-Newton Verfahren an Genauigkeit, wodurch ihre Konvergenzgeschwindigkeit superlinear anstatt quadratisch ist. Dafür gewinnen sie zusätzliche Geschwindigkeit durch die Vermeidung der Bestimmung und Inversion von \(\nabla^2 f(x_k)\). Der Vorteil der Quasi-Newton Methoden ist es, dass man nur den Gradienten \(\nabla f\) für einen Schritt des numerischen Optimierungsverfahrens benötigt und keine expliziten Informationen über die zweiten Ableitungen. Dadurch werden sie in bestimmten Problemen sogar effizienter bei der Approximation eines stationären Punktes als das Newton Abstiegsverfahren in Newton Verfahren.
3.2.5.1. Sekantengleichung und Krümmungsbedingung#
Die entscheidende Frage bei der Konstruktion eines Quasi-Newton Abstiegsverfahrens der Form (3.12) ist es, wie die positiv definite Matrix \(B_k\) in jedem Schritt möglichst effizient bestimmt werden kann. Anstatt die Näherung \(B_k\) der Hessematrix \(\nabla^2 f(x_k)\) in jedem Schritt von Grund auf neu zu berechnen, wäre es wünschenswert ein initiales \(B_0\) zu bestimmen, das in jedem Schritt des Iterationsverfahrens nur aktualisiert werden muss. Hierbei ist es möglich die durch den Iterationsschritt erhaltenen Informationen über den Gradienten \(\nabla f\) zu Hilfe zu nehmen.
Wir nehmen an, wir haben bereits einen Abstiegsschritt durchgeführt und so einen neuen Punkt \(x_{k+1} = x_k + \tau_kp\) erhalten. Unsere quadratische Approximation in diesem neuen Punkt für eine neue Richtung \(p \in \mathbb{R}^n\) sieht dementsprechend wie folgt aus:
Es ist leicht einzusehen, dass die Modellfunktion \(m_{k+1}\) im Punkt \(x_{k+1} \in \R^n\) zentriert ist und für \(p=0\) mit dem Funktionswert der Zielfunktion im Punkt \(x_{k+1}\) übereinstimmt, d.h., es gilt \(m_{k+1}(0) = f(x_{k+1})\).
Da wir uns bei der Wahl der positiv definiten Matrix \(B_{k+1}\) noch nicht festgelegt haben, wird durch (3.13) eine Funktionenschar beschrieben. Eine Forderung, die man nun die Modellfunktion \(m_{k+1}\) stellen kann, um eine sinnvolle Matrix \(B_{k+1}\) zu bestimmen, ist, dass ihre Ableitung \(\nabla m_{k+1}\) mit der Ableitung der Zielfunktion \(f\) in den letzten beiden Punkten \(x_k\) und \(x_{k+1}\) übereinstimmt. Dies bedeutet, dass man die Matrix \(B_{k+1}\) versucht so zu bestimmen, dass die Modellfunktion \(m_{k+1}\) die Krümmung der Zielfunktion \(f\) gut approximiert. Da bereits gilt
ist eine der beiden Forderungen automatisch erfüllt. Für die zweite Forderung können wir nutzen, dass \(x_k = x_{k+1} - \tau_k p_k\) gilt und wir erhalten somit:
Durch Umstellen von (3.14) erhalten wir die Bedingung
Eine vernünftige Wahl der Matrix \(B_{k+1}\) in (3.11) sollte versuchen, diese Eigenschaft, auch bekannt als Sekantengleichung, zu imitieren. Im eindimensionalen Fall mit \(f \colon \Omega \subset \mathbb{R} \rightarrow \mathbb{R}\) bedeutet die Sekantengleichung nichts anderes, als dass der Faktor \(B_{k+1}\) eine Approximation des zweiten Ableitung von \(f\) im Sinne eines Differenzenquotienten ist, d.h., im Fall \(n=1\) soll gelten:
Für unser allgemeines Quasi-Newton Verfahren in (3.12) suchen wir also einen Weg die bereits bekannte Approximation der Hessematrix \(B_k \approx \nabla^2f(x_k)\) zu einer Matrix \(B_{k+1}\) zu aktualisieren, so dass der folgende Zusammenhang für den nächsten Punkt \(x_{k+1} \in \Omega\) erfüllt wird:
wobei
Es wird klar, dass diese Forderung alleine nicht für die Konstruktion eines Abstiegsverfahrens genügt, da die Sekantengleichung in (3.15) für \(n > 1\) unterbestimmt ist, d.h., dass es mehr unbekannte Einträge der Matrix \(B_{k+1} \in \R^{n \times n}\) gibt als durch die \(n\) Gleichungen festgelegt werden. Daher versuchen wir im Folgenden weitere Forderungen an die Matrix \(B_{k+1}\) zu stellen.
Um die positive Definitheit der Matrix \(B_{k+1}\) in Schrittrichtung \(x_{k+1} - x_k = \tau_k p_k \in \R^n\) zu gewährleisten müssen wir fordern, dass die Vektoren \(y_k\) und \(s_k\) die sogenannte Krümmungsbedingung erfüllen:
Dies ist eine hinreichende Bedingung für die positive Definitheit von \(B_{k+1}\) bezüglich der Richtung \(\tau_k p_k\), da wir einfach die Sekantengleichung (3.15) von links mit dem Vektor \(s_k^\top\) multiplizieren können und so erhalten wir mit der Forderung (3.16) schon:
Bemerkung 3.8 (Krümmungsbedingung und Konvexität)
Falls die Zielunktion \(f\) strikt konvex ist, so ist die Krümmungsbedingung (3.16) für alle Punktepaare \(x_k, x_{k+1} \in \Omega\) erfüllt und die Matrix \(B_{k+1}\) wird damit positiv definit. Für nichtkonvexe Funktionen hingegen muss man die Krümmungsbedingung explizit forcieren, um ein Abstiegsverfahren zu erhalten.
Falls die Krümmungsbedingung (3.16) erfüllt ist, so existiert mindestens eine Lösung \(B_{k+1}\) der Sekantengleichung (3.15). Man sieht ein, dass es in der Tat sogar unendlich viele Lösungen \(B_{k+1}\) gibt, da eine symmetrische \(n \times n\) Matrix \(n(n+1)/2\) Freiheitsgrade besitzt und die Sekantengleichung (3.15) nur \(n\) Bedingungen an \(B_{k+1}\) stellt. Zusätzlich erhält man \(n\) Bedingungen an \(B_{k+1}\) durch die Forderung von positiver Definitheit, da alle \(n\) Hauptminoren von \(B_{k+1}\) positiv sein müssen. Dies reicht jedoch nicht für die eindeutige Bestimmung der Matrix \(B_{k+1}\). Hierfür müssen wir zusätzlich fordern, dass die Matrix \(B_{k+1}\) diejenige Matrix unter allen möglichen Lösungen ist, die der vorherigen Matrix \(B_k\) am nächsten bezüglich eines geeigneten Maßes ist. Das heißt wir suchen eine Lösung des folgenden Optimierungsproblems:
wobei \(s_k\) und \(y_k\) definiert sind wie in der Sekantengleichung (3.15). Man beachte, dass man eine unterschiedliche Lösung \(B_{k+1}\) des Optimierungsproblems (3.17) in Abhängigkeit der gewählten Matrixnorm erhält und somit auch ein unterschiedliches Quasi-Newton Verfahren herleiten kann.
3.2.5.2. Das Davidon-Fletcher-Powell Verfahren#
Im ursprünglich im Jahr \(1959\) von Davidon vorgeschlagenen Verfahren [Dav59] (das im Übrigen bei der Erstbegutachtung abgelehnt wurde) wählt man für die Norm im Optimierungsproblem (3.17) eine gewichtete Frobeniusnorm der Form
Die Gewichtungsmatrix \(W\) dient dazu, dass das implizierte Quasi-Newton Verfahren zur Approximation eines stationären Punktes \(x^* \in \Omega\) skalierungs-invariant wird. Hierzu wählt man eine beliebige Matrix für die die Relation \(W y_k = s_k\) gilt, d.h., eine Matrix \(W\), die sich wie die Inverse der Matrix \(B\) in (3.17) verhält. Ein konkretes Beispiel für solch eine Gewichtungsmatrix wäre \(W \coloneqq G_k^{-1}\), wobei \(G_k\) die durchschnittliche Hessematrix von \(f\) entlang des letzten Abstiegsschritts von \(x_k \rightarrow x_{k+1}\) ist mit
Mit der konkreten Wahl dieser Gewichtungsmatrix \(W = G_k^{-1}\) wird die gewichtete Frobeniusnorm dimensionslos und man erhält eine eindeutige Lösung des Optimierungsproblems (3.17) wie folgt:
Die Gleichung (3.18) wird auch DFP-Schritt genannt, da sie zuerst von Davidon vorgeschlagen und später von Fletcher und Powell untersucht und verbreitet wurde.
Obwohl wir die explizite Berechnung der Hessematrix \(\nabla^2 f(x_k)\) vermieden haben und die Aktualisierung der Matrix \(B_k\) zu \(B_{k+1}\) lediglich auf den Gradienteninformationen von \(f\) basiert ist der numerische Aufwand bei direkter Verwendung von \(B_{k+1}\) in (3.18) noch zu hoch. Das liegt daran, dass wir für einen Schritt des Quasi-Newton Verfahrens in (3.12) die Inverse der Matrix \(B_k\) benötigen und die Inversion einen numerischen Aufwand von \(\mathcal{O}(n^3)\) besitzt. Glücklicherweise gibt es einen Trick, wie wir die Inverse von \(B_k\) in jedem Schritt des Iterationsverfahrens numerisch günstig erhalten können. Sei \(H_k \coloneqq B_k^{-1}\), dann können wir die sogenannte Sherman-Morrison-Woodbury Formel (siehe [Wikb]) auf Gleichung (3.18) anwenden um die neue Inverse \(H_{k+1}\) durch eine Aktualisierung der Matrix \(H_k\) zu berechnen:
Wie man einssieht liegt der numerische Rechenaufwand für das Update von \(H_{k+1}\) in (3.19) in \(\mathcal{O}(n^2)\). Es fällt außerdem auf, dass \(H_k\) nur durch die Addition zweier Matrizen mit Rang \(1\) verändert wird, also insgesamt eine Änderung von höchstens Rang \(2\) erfährt. Das passt gut zu der Forderung, dass wir erwarten, dass sich die Approximation der Hessematrix \(\nabla^2 f\) in einer lokalen Umgebung nur wenig ändert.
3.2.5.3. Das Broyden–Fletcher–Goldfarb–Shanno Verfahren#
Das Davidon–Fletcher–Powell Verfahren wurde trotz seiner Effektivität bald schon durch ein Verfahren abgelöst, das noch besser war und bis heute zu den effizientesten Quasi-Newton Verfahren gehört: das Broyden–Fletcher–Goldfard–Shanno (BFGS) Verfahren in [Bro70]. Die Idee des BFGS Verfahrens leitet sich unmittelbar aus der Idee des DFP Verfahren ab. Anstatt das Optimierungsproblem (3.17) mit bestimmten Bedingungen an die Approximation \(B_{k+1}\) der Hessematrix \(\nabla^2 f(x_k)\) zu stellen, versucht man direkt die Inverse der Hessematrix \((\nabla^2 f(x_k))^{-1}\) geeignet zu approximieren. Hierfür nehmen wir an, dass wir eine Matrix \(H_{k+1}\) als geringfügige Aktualisierung einer bereits vorher bestimmten Matrix \(H_k\) suchen, die gleichzeitig symmetrisch und positiv definit ist und zusätzlich die Sekantenbedingung in umgeschriebener Form erfüllt:
Hierzu formuliert man ein analoges Optimierungsproblem zu (3.17) von der Form:
Unter der Verwendung der gewichteten Frobeniusnorm und einer beliebigen Gewichtsfunktion, die die Sekantengleichung \(Ws_k = y_k\) erfüllt, erhält man wiederum die eindeutige Lösung des Minimierungsproblems (3.20) als:
Das Update der Matrix \(H_k\) in (3.21) kann numerisch in \(\mathcal{O}(n^2)\) durchgeführt werden, was man schnell einsieht, wenn man das Produkt ausschreibt:
In dieser Schreibweise sieht man gut, dass man lediglich Skalarprodukte in \(\mathcal{O}(n)\), Matrix-Vektor Multiplikationen in \(\mathcal{O}(n^2)\) und dyadische Produkte in \(\mathcal{O}(n^2)\) berechnen muss. Im Gegensatz hierzu würde eine naive Implementierung des BFGS-Updates in (3.21) zu einem numerischen Aufwand von \(\mathcal{O}(n^3)\) führen.
Abschließend bleibt die Frage was eine gute Initialisierung der Matrix \(H_0\) ist. Idealerweise hat man bereits Informationen über die Inverse der Hessematrix \((\nabla^2 f(x_0))^{-1}\) im Initialisierungspunkt \(x_0 \in \Omega\), zum Beispiel durch eine numerische Approximation mittels finiter Differenzen (später in der Vorlesung!). Andererseits erwarten wir, dass die Aktualisierung von \(H_k\) im \(k\)-ten Schritt des Iterationsverfahrens (3.21) zu \(H_{k+1}\) die aktuellen Informationen über den Verlauf der Gradienten \(\nabla f(x_k)\) und \(\nabla f(x_{k+1})\) berücksichtigt. Darum ist eine häufige Wahl von \(H_0\) die Initialisierung als Einheitsmatrix \(I_n\) oder ein Vielfaches der Einheitsmatrix, wobei die Vorfaktoren der Diagonaleinträge entsprechend der Skalierung der Variablen gewählt werden.