Benutzer:Stepri2005/Kurs:Optimierung II/Quasi-Newton-Verfahren

7.1 Einleitung

Wir beschäftigen wir uns weiterhin mit der Lösung des unrestringierten Optimierungsproblems

$(P):{\text{ Minimiere }}f(x){\mbox{ über alle }}x\in \mathbb {R} ^{n}.$

In Kapitel 6 hatten wir unter geeigneten Voraussetzungen die superlineare und quadratische Konvergenz für das lokale und das globalisierte Newton-Verfahren nachgewiesen. In Abschnitt 6.3 hatten wir aber auch Nachteile des Newton-Verfahrens zumindest für größere Probleme erwähnt, wie insbesondere die, dass in jeder Iteration des Verfahrens die Hesse-Matrix bestimmt und ein lineares Gleichungssystem gelöst werden müssen.

Diese Nachteile versucht man mit den Quasi-Newton-Verfahren zu vermeiden, ohne dabei die superlineare Konvergenzrate des Newton-Verfahrens zu verlieren. Das Newton-Verfahren konvergiert zwar unter gewissen Voraussetzungen sogar quadratisch, aber für die meisten praktischen Probleme ist eine superlineare, nicht notwendig quadratische Konvergenzrate ausreichend. Überdies ist der numerische Aufwand pro Iteration für Quasi-Newton-Verfahren erheblich geringer als der für das Newton-Verfahren.

Wir definieren hier wiederum

$g^{k}:=\nabla f(x^{k})$

und setzen in unserer Diskussion der Einfachheit halber voraus, dass $f\in C^{2}(\mathbb {R} ^{n})$ ist. (In den Lemmata und Sätzen geben wir jeweils, wie zuvor, die genauen Voraussetzungen an.) Die Quasi-Newton-Verfahren, mit denen wir uns beschäftigen werden, haben dann die folgende allgemeine Form.

Modellalgorithmus 7.1 (Quasi-Newton-Modellalgorithmus)

(0) Wähle $x^{0}\in \mathbb {R} ^{n}$, eine symmetrische, positiv definite Matrix $B_{0}\in \mathbb {R} ^{n\times n}$ und eine semieffiziente Schrittweitenregel. Setze $k:=0$.

(1) Falls $g^{k}=0$ ist, stop! ($x^{k}$ ist kritische Lösung von Problem $(P)$.)

(2) Berechne

(7.1) $p^{k}:=-B_{k}g^{k},$

eine Schrittweite $t_{k}>0$ und setze

$x^{k+1}:=x^{k}+t_{k}p^{k}.$

(3) Setze

$s^{k}:=x^{k+1}-x^{k},\quad y^{k}:=g^{k+1}-g^{k}$

und bestimme aus $B_{k},s^{k}$ und $y^{k}$ mittels einer Aufdatierungsformel eine symmetrische, positiv definite Matrix $B_{k+1}$.

(4) Setze $k:=k+1$ und gehe nach (1).

Da $B_{k}$ im Modellalgorithmus 7.1 symmetrisch und positiv definit ist und in Schritt (2) $g^{k}\neq 0$ gilt, folgt

(7.2) $(g^{k})^{T}p^{k}=-(g^{k})^{T}B_{k}g^{k}<0.$

Somit ist $p^{k}$ eine Abstiegsrichtung von $f$ in $x^{k}$ und gilt damit $f(x^{k}+t_{k}p^{k})<f(x^{k})$, da eine semieffiziente Schrittweite verwendet wird. Bei Algorithmen dieses Typs versucht man offenbar das Newton-Verfahren zu simulieren, indem man $\left[\nabla ^{2}f(x^{k})\right]^{-1}$ durch eine geeignete Matrix $B_{k}$ ersetzt (normalerweise schreibt man $B_{k}$ und nicht $B_{k}^{-1}$). Daher bezeichnet man solche Verfahren als Quasi-Newton-Verfahren.

Bemerkung 7.2

Gelegentlich werden Verfahren vom Typ des Modellalgorithmus 7.1 auch als Variable-Metrik-Verfahren bezeichnet. Dies ist dadurch begründet, dass man die Richtung $p^{k}:=-B_{k}g^{k}$ gemäß Bemerkung 2.4 als Richtung des steilsten Abstiegs für $f$ in $x^{k}$ bezüglich der Norm $\|\cdot \|_{B_{k}^{-1}}$ interpretieren kann. In dem Modellalgorithmus 7.1 wählt man demnach in jeder Iteration die Richtungsteilsten Abstiegs bezüglich einer von $k$ abhängenden Metrik. Wird $B_{0}:=I$ gewählt, so ist $p^{0}:=-g^{0}$ gerade die Richtung steilsten Abstiegs in $x^{0}$ bezüglich der Euklidischen Norm.

Die zentrale Frage ist nun, wie man aus den in der $k$-ten Iteration vorliegenden Daten $B_{k},s^{k}$ und $y^{k}$ eine geeignete Matrix $B_{k+1}$ berechnet. Man stellt diesbezüglich die folgenden drei Forderungen an $B_{k+1}$ auf:

(F1) Es gilt

$B_{k+1}:=\varrho _{k}B_{k}+T_{k}$

für eine Konstante $\varrho _{k}\neq 0$ und eine Matrix $T_{k}\in \mathbb {R} ^{n\times n}$ mit

$1\leq \operatorname {Rang} (T_{k})\leq 2.$

(F2) Es gilt die Quasi-Newton-Gleichung

(7.3) $B_{k+1}y^{k}=s^{k}.$

(F3) Mit $B_{k}$ ist auch $B_{k+1}$ eine symmetrische, positiv definite Matrix.

Eine Formel wie in (F1), gemäß der man aus Größen, die in der $k$-ten Iteration bekannt sind, eine Größe für die $(k+1)$-te Iteration eines Verfahrens bestimmt, bezeichnet man als Update-Formel. (Das englische „Update“ hat sich allgemein gegenüber der deutschen „Aufdatierung“ durchgesetzt, so dass wir es hier auch verwenden wollen.) Je nachdem, ob $T_{k}$ in (F1) normalerweise den Rang 1 oder den Rang 2 hat, spricht man von einem Rang-1- oder Rang-2-Verfahren. Die Berechnung einer Rang-1-Matrix erfordert weniger Rechenoperationen als die einer Rang-2-Matrix. Auf der anderen Seite kann man aber erwarten, dass man mittels einer Rang-2-Matrix in einem Verfahren vom Typ des Modellalgorithmus 7.1 mit wachsendem $k$ schneller eine Annäherung an das Newton-Verfahren erreichen kann.

Die Forderung (F2) lässt sich dabei auf folgende Weisen motivieren, sofern $B_{k+1}$ eine reguläre Matrix ist. Für jede quadratische Funktion

$f(x):={\frac {1}{2}}x^{T}Qx+c^{T}x+\alpha$

mit symmetrischer, positiv definiter Matrix $Q$ hat man

$\nabla f(x)=Qx+c,\quad \nabla ^{2}f(x)=Q$

und somit

$\left[\nabla ^{2}f(x^{k+1})\right]^{-1}y^{k}=Q^{-1}\left(Qx^{k+1}-Qx^{k}\right)=x^{k+1}-x^{k}=s^{k}.$

Demnach genügt $B_{k+1}:=\left[\nabla ^{2}f(x^{k+1})\right]^{-1}$ der Quasi-Newton-Gleichung. Außerdem folgt für eine beliebige Funktion $f$ mit der quadratischen Näherung

$F_{k+1}(x):=f(x^{k+1})+\nabla f(x^{k+1})^{T}\left(x-x^{k+1}\right)+{\frac {1}{2}}\left(x-x^{k+1}\right)^{T}B_{k+1}^{-1}\left(x-x^{k+1}\right)$

im Fall, dass die Quasi-Newton-Gleichung erfüllt und $B_{k+1}$ eine reguläre Matrix ist:

$F_{k+1}(x^{k+1})=f(x^{k+1}),\quad \nabla F_{k+1}(x^{k+1})=\nabla f(x^{k+1}),$

$\nabla F_{k+1}(x^{k})=\nabla f(x^{k+1})+B_{k+1}^{-1}\left(x^{k}-x^{k+1}\right)=g^{k+1}-y^{k}=\nabla f(x^{k}).$

Diese Beziehungen deuten darauf hin, dass unter den genannten Bedingungen $F_{k+1}$ eine gute quadratische Approximation für $f$ und folglich $B_{k+1}$ eine gute Näherung für $\left[\nabla ^{2}f(x^{k+1})\right]^{-1}$ ist.

Bemerkung 7.3

Im Fall $n=1$ besagt die Quasi-Newton-Gleichung:

$B_{k+1}^{-1}={\frac {y^{k}}{s^{k}}}={\frac {f'(x^{k+1})-f'(x^{k})}{x^{k+1}-x^{k}}}.$

In diesem Fall entspricht also $B_{k+1}^{-1}\in \mathbb {R}$ gerade der Steigung der Sekante für $f'$ durch die Punkte $(x_{k+1},f'(x_{k+1}))$ und $(x_{k},f'(x_{k}))$ und ist somit $B_{k+1}^{-1}\approx f''(x_{k+1})$. Wählt man im Modellalgorithmus 7.1 für jedes $k$ die Schrittweite $t_{k}:=1$, so stimmt dieser also für $n=1$ mit dem aus der Numerischen Mathematik bekannten Sekanten-Verfahren zur Bestimmung einer Nullstelle der Gleichung $f'(x)=0$ überein. Quasi-Newton-Verfahren verallgemeinern demnach in gewisser Weise das Sekantenverfahren für Funktionen in einer Veränderlichen auf solche in n Veränderlichen, weswegen sie gelegentlich auch als Sekanten-Verfahren bezeichnet werden. Entsprechend spricht man bei (F2) auch manchmal von der Sekanten-Gleichung.

In Abschnitt 7.2 werden wir allgemein Rang-1-Update-Formeln für Quasi-Newton-Verfahren aus den Forderungen (F1) - (F3) ableiten. Im Fall von Rang-2-Update-Formeln werden wir nicht so axiomatisch vorgehen, sondern in Abschnitt 7.3 gleich die wichtige Klasse der Broyden-Update-Formeln angeben und genauer betrachten. Die Formel aus dieser Klasse, die sich im Laufe der Zeit als die vom numerischen Verhalten her in den meisten Situationen beste Formel herausgestellt hat, ist die sog. BFGS-Update-Formel. Die Konvergenz des dadurch definierten BFGS-Verfahrens sowie Aussagen über dessen Konvergenzgeschwindigkeit werden wir in Abschnitt 7.4 beweisen. Mit einigen Bemerkungen zur Numerik werden wir in Abschnitt 7.5 dieses Kapitel abschließen.

7.2 Rang-1-Update-Formeln

Als erstes wollen wir nun untersuchen, inwieweit die drei Forderungen (F1) - (F3) im Fall $\operatorname {Rang} (T_{k})=1$ eine Update-Formel für $B_{k+1}$ festlegen. Dazu leiten wir zunächst Gestalt und Eigenschaften von Rang-1-Matrizen her.

Lemma 7.4

Sei $A\in \mathbb {R} ^{n\times n}$ eine Matrix mit $\operatorname {Rang} (A)=1$.

(i) Dann gilt $A=uv^{T}$ für gewisse $u,v\in \mathbb {R} ^{n}\setminus \{0\}$.

(ii) Ist $A=A^{T}$, dann gilt $A=\varepsilon ww^{T}$ für ein $w\in \mathbb {R} ^{n}\setminus \{0\}$ und $\varepsilon \in \{-1,1\}$.

Beweis.

Nach Voraussetzung hat der Bildraum von $A$ die Dimension 1. Demnach gibt es einen Vektor $u\neq 0$, so dass für jedes $x\in \mathbb {R} ^{n}$ mit einem $v_{x}\in \mathbb {R}$ gilt: $Ax=v_{x}u$. Ist $e^{j}$ der $j$-te Standardeinheitsvektor, so folgt insbesondere für $x:=e^{j}$ mit gewissen $v_{j}\in \mathbb {R}$

(7.4) $a^{j}=v_{j}u\quad (j=1,\ldots ,n),$

wobei $a^{j}$ die $j$-te Spalte von $A$ ist. Also hat man mit $v:=(v_{1},\ldots ,v_{n})^{T}$

$A=(u_{i}v_{j})_{i,j=1,\ldots ,n}=uv^{T}.$

Offenbar ist $v\neq 0$, da anderenfalls nach (7.4) $A=0$ und nicht $\operatorname {Rang} (A)=1$ wäre.

Im Fall $A=A^{T}$ folgt somit $uv^{T}=vu^{T}$ und mit $\alpha :=(u^{T}v/v^{T}v)\neq 0$ daher

$uv^{T}v=vu^{T}v\Leftrightarrow u={\frac {u^{T}v}{v^{T}v}}v=\alpha v.$

Also ergibt sich für $w:={\sqrt {|\alpha |}}v$

$A=uv^{T}=vv^{T}=\operatorname {sgn}(\alpha )\left({\sqrt {|\alpha |}}v\right)\left({\sqrt {|\alpha |}}v\right)^{T}=\operatorname {sgn}(\alpha )ww^{T}.$

q.e.d.

Man beachte, dass eine symmetrische Matrix maximal $(n^{2}+n)/2$ unterschiedliche Elemente enthält, so dass das Aufstellen einer Rang-1-Matrix

$uv^{T}=(u_{i}v_{j})_{i,j=1,\ldots ,n}$

die Berechnung von höchstens $(n^{2}+n)/2$ reellen Produkten erfordert. Wir benötigen weiter folgende Aussagen für Rang-1-Matrizen, wobei die Spur einer Matrix die Summe ihrer Diagonalelemente ist:

Lemma 7.5

Es seien $u,v\in \mathbb {R} ^{n}\setminus \{0\}$. Dann gilt:

(i) Die Matrix $uv^{T}$ hat im Fall $u^{T}v=0$ den $n$-fachen Eigenwert 0 und im Fall $u^{T}v\neq 0$ den $(n-1)$-fachen Eigenwert 0 und den 1-fachen Eigenwert $u^{T}v$.

(ii) $\operatorname {Spur} (uv^{T})=u^{T}v.$

(iii) $\det(I+uv^{T})=1+u^{T}v.$

Beweis.

Nach Voraussetzung hat man $u\neq 0$ und $v\neq 0$. Sei nun $\lambda$ ein Eigenwert von $uv^{T}$ mit zugehörigem Eigenvektor $x$, d. h., es sei

(7.5) $(uv^{T})x=\lambda x\Leftrightarrow (v^{T}x)u=\lambda x.$

Dann kann also nur $x$ orthogonal zu $v$ und damit $\lambda =0$ sein oder (nicht im ausschließenden Sinne) es sind $x$ und $u$ parallel zueinander, wobei dann $\lambda =0$ oder $\lambda \neq 0$ möglich ist.

Ist $x$ orthogonal zu $v$, so ist $\lambda =0$ und die Menge

$\left\{x\in \mathbb {R} ^{n}{\big |}\left(uv^{T}\right)x=0\right\}$

der zugehörige Eigenraum. Wegen $\operatorname {Rang} (uv^{T})=1$ hat dieser die Dimension $n-1$. Da die geometrische Vielfachheit eines Eigenwertes kleiner oder gleich seiner algebraischen Vielfachheit ist, besitzt der Eigenwert $\lambda =0$ von $uv^{T}$ somit mindestens die Vielfachheit $n-1$. Sind andererseits $x$ und $u$ parallel zueinander, d. h. gilt $x=\mu u$ für ein $\mu \neq 0$, so folgt aus (7.5) die Beziehung

$\left(uv^{T}\right)u=\lambda u.$

Wegen $u\neq 0$ schließt man daraus, dass $\lambda =v^{T}u$ ist (wobei $v^{T}u=0$ möglich ist).

Aussage (ii) folgt wegen

$\operatorname {Spur} (uv^{T})=\sum _{i=1}^{n}u_{i}v_{i}=u^{T}v.$

Schließlich hat die Matrix $I+uv^{T}$ im Fall $u^{T}v=0$ den $n$-fachen Eigenwert 1 und im Fall $u^{T}v\neq 0$ den $(n-1)$-fachen Eigenwert 1 und den 1-fachen Eigenwert $1+v^{T}u$. Da die Determinante einer reellen $(n\times n)$-Matrix gleich dem Produkt ihrer $n$ Eigenwerte ist, folgt damit das gewünschte Ergebnis in (iii).

q.e.d.

Ist die Bedingung (F1) für eine Matrix $T_{k}$ mit $\operatorname {Rang} (T_{k})=1$ erfüllt und soll $B_{k+1}$, wie es in (F3) gefordert wird, symmetrisch sein, so muss auch $T_{k}$ symmetrisch sein und muss gemäß Lemma 7.4 mit $\varrho _{k}\neq 0,\varepsilon _{k}\in \{-1,1\}$ und einem Vektor $v^{k}\neq 0$ gelten:

(7.6) $B_{k+1}:=\varrho _{k}B_{k}+\varepsilon _{k}v^{k}(v^{k})^{T}.$

Die Quasi-Newton-Gleichung (F2) hat demnach in diesem Fall die Gestalt

$\varrho _{k}B_{k}y^{k}+\varepsilon _{k}v^{k}(v^{k})^{T}y^{k}=s^{k}.$

Unter der Voraussetzung

(7.7) $(v^{k})^{T}y^{k}\neq 0\Leftrightarrow s^{k}-\varrho _{k}B_{k}y^{k}\neq 0$

impliziert die letztere Beziehung

(7.8) $v^{k}=\lambda _{k}\left(s^{k}-\varrho _{k}B_{k}y^{k}\right),\quad \lambda _{k}:={\frac {1}{\varepsilon _{k}(v^{k})^{T}y^{k}}},$

womit man erhält:

(7.9) $1=\lambda _{k}\varepsilon _{k}(v^{k})^{T}y^{k}=\varepsilon _{k}\lambda _{k}^{2}\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}y^{k}.$

Mit (7.8) und (7.9) bekommt man schließlich

$\varepsilon _{k}v^{k}(v^{k})^{T}=\varepsilon _{k}\lambda _{k}^{2}\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}={\frac {\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}}{\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}y^{k}}}.$

Also ergibt sich aus (7.6) die Formel

(7.10) $B_{k+1}:=\varrho _{k}B_{k}+{\frac {\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}}{\left(s^{k}-\varrho _{k}B_{k}y^{k}\right)^{T}y^{k}}}.$

Für $\varrho _{k}:=1$ erhält man aus (7.10) gerade die Broyden-Rang-1- bzw. SR1-Formel („symmetric-rank-1 formula“)

(7.10) $B_{k+1}:=B_{k}+{\frac {\left(s^{k}-B_{k}y^{k}\right)\left(s^{k}-B_{k}y^{k}\right)^{T}}{\left(s^{k}-B_{k}y^{k}\right)^{T}y^{k}}}.$

Diese Formel hat den Vorteil, eine recht einfache Update-Formel zu sein. Für Algorithmen mit Schrittweitenstrategien, wie wir sie hier untersuchen, hat sie jedoch den schwerwiegenden Nachteil, dass nicht auszuschließen ist, dass der Nenner darin verschwindet und dass die Forderung (F3) der positiven Definitheit von B_{k+1} nicht erfüllt ist. (Insbesondere ist offenbar der Nenner in (7.11) identisch Null und dann diese Update-Formel nicht definiert, wenn die Voraussetzung (7.7) für den hier gewählten Weg ihrer Herleitung nicht erfüllt ist.)

Trotz einiger attraktiver Eigenschaften, wie z. B. der endlichen Konvergenz bei quadratischer Zielfunktion auch bei inexakten Schrittweiten, wenn das Verfahren durchführbar ist, wollen wir daher das Quasi-Newton-Verfahren mit der SR1-Formel nicht weiter untersuchen. Wir verweisen aber darauf, dass sich herausgestellt hat, dass die (modifizierte) SR1-Formel in manchen Zusammenhängen, insbesondere im Zusammenhang mit Trust-Region-Verfahren, die wir in Kapitel 8 diskutieren werden, hervorragende und oft bessere Ergebnisse als die unten untersuchte BFGS-Formel liefert (siehe [NoWri06] und die Literaturverweise in [GeiKa99]).

Wir wollen nun als nächstes der Frage nachgehen, ob nicht für $\varrho _{k}\neq 1$ eine Rang-1-Update-Formel vom Typ (7.11) existiert, für die auch (F3) erfüllt ist. Dazu führen wir zunächst für eine symmetrische, positiv definite Matrix $A\in \mathbb {R} ^{n\times n}$ auf folgende Weise eine ebenfalls symmetrische, positiv definite Matrix $A^{1/2}\in \mathbb {R} ^{n\times n}$ ein: es sei $\langle x,y\rangle _{A}$ das durch $A$ definierte Skalarprodukt

(7.12) $\langle x,y\rangle _{A}:=x^{T}Ay,\quad x,y\in \mathbb {R} ^{n}$

und

$\|x\|_{A}:=sqrt{\langle x,x\rangle _{A}}:={\sqrt {x^{T}Ax}},\quad x\in \mathbb {R} ^{n}$

die durch dieses Skalarprodukt induzierte elliptische Norm (vgl. Bemerkung 2.4). Weiter seien $D:=\operatorname {diag} (d_{1},\ldots ,d_{n})$ mit $d_{i}>0$ und $P\in \mathbb {R} ^{n\times n}$ die zu $A$ existierende Diagonalmatrix und orthogonale Matrix, so dass $P^{-1}AP=D$ gilt. Dann definieren wir

$D^{1/2}:=\operatorname {diag} ({\sqrt {d_{1}}},\ldots ,{\sqrt {d_{n}}})$

und

$A^{1/2}:=P^{-1}D^{1/2}P.$

Die Matrix $A^{1/2}$ ist wegen $P^{-1}=P^{T}$ ebenfalls symmetrisch und, da sie die Eigenwerte ${\sqrt {d_{i}}}>0$ hat, ebenfalls positiv definit. Ferner folgt für sie

$A^{1/2}A^{1/2}=\left(P^{-1}D^{1/2}P\right)\left(P^{-1}D^{1/2}P\right)=P^{-1}DP=A$

sowie

$A^{-1/2}:=\left(A^{1/2}\right)^{-1}=P^{-1}D^{-1/2}P.$

Für die durch den Modellalgorithmus 7.1 erzeugten Größen gilt nun

$s^{k}=x^{k+1}-x^{k}=t_{k}p^{k}$

und mit $\eta _{k}:=t_{k}/\varrho _{k}$

(7.13) $z^{k}:=y^{k}+\eta _{k}g^{k}=y^{k}+{\frac {t_{k}}{\varrho _{k}}}g^{k}.$

Wegen $B_{k}g^{k}=-p^{k}$ hat man folglich

$B_{k}z^{k}=B_{k}y^{k}-{\frac {t_{k}}{\varrho _{k}}}p^{k}=B_{k}y^{k}{\frac {1}{\varrho _{k}}}s^{k}.$

Also impliziert (7.10) unter Verwendung der Beziehung $p^{k}=-B_{k}g^{k}$

(7.14) $B_{k+1}=\varrho _{k}B_{k}-{\frac {\varrho _{k}^{2}\left(B_{k}y^{k}-{\frac {1}{\varrho _{k}}}s^{k}\right)\left(B_{k}y^{k}-{\frac {1}{\varrho _{k}}}s^{k}\right)^{T}}{\varrho _{k}\left(B_{k}y^{k}-{\frac {1}{\varrho _{k}}}s^{k}\right)^{T}y^{k}}}={\frac {t_{k}}{\eta _{k}}}\left\{B_{k}-{\frac {\left(B_{k}z^{k}\right)\left(B_{k}z^{k}\right)^{T}}{\left(B_{k}z^{k}\right)^{T}\left(z^{k}-\eta _{k}g^{k}\right)}}\right\}$

$={\frac {t_{k}}{\eta _{k}}}B_{k}^{1/2}\left\{I-{\frac {\left(B_{k}^{1/2}z^{k}\right)\left(B_{k}^{1/2}z^{k}\right)^{T}}{(z^{k})^{T}B_{k}z^{k}+\eta _{k}(z^{k})^{T}p^{k}}}\right\}B_{k}^{1/2}={\frac {t_{k}}{\eta _{k}}}B_{k}^{1/2}K_{k}B_{k}^{1/2}$

mit

(7.15) $K_{k}:=I-{\frac {\left(B_{k}^{1/2}z^{k}\right)\left(B_{k}^{1/2}z^{k}\right)^{T}}{(z^{k})^{T}B_{k}z^{k}+\eta _{k}(z^{k})^{T}p^{k}}}$

Im Fall $\varrho _{k}=t_{k}/\eta _{k}>0$ ist $B_{k+1}$ offenbar genau dann positiv definit, wenn $K_{k}$ dies ist.

Also stellt sich die Frage, wann $K_{k}$ positiv definit ist. Die Matrix $K_{k}$ hat gemäß Lemma 7.5 $(n-1)$-mal den Eigenwert 1 sowie den Eigenwert

(7.16) $\nu _{k}:=1-{\frac {(z^{k})^{T}B_{k}z^{k}}{(z^{k})^{T}B_{k}z^{k}+\eta _{k}(z^{k})^{T}p^{k}}}$

Also ist $K_{k}$ genau dann positiv definit, wenn $\nu _{k}>0$ gilt. Insbesondere ist dies unter den Voraussetzungen des folgenden Lemmas der Fall.

Lemma 7.6

Es seien $f\in C^{1}(\mathbb {R} ^{n})$ sowie $\eta _{k}\in (0,1)$ und es sei $t_{k}>0$ eine exakte Schrittweite. Dann folgt $s^{k}-\varrho _{k}B_{k}y^{k}\neq 0$ (vgl. (7.7)) sowie $\nu _{k}>0$.

Beweis.

Gemäß Lemma 5.9 hat man

$\nabla f(x^{k}+t_{k}p^{k})^{T}p^{k}=(g^{k+1})^{T}p^{k}=0.$

Unter Verwendung von (7.2) und (7.13) ergibt sicht somit für $\eta _{k}\in (0,1)$

$(z^{k})^{T}p^{k}=(g^{k+1}-g^{k}+\eta _{k}g^{k})^{T}p^{k}=-(1-\eta _{k})(g^{k})^{T}p^{k}>0.$

Da $B_{k}$ positiv definit ist, gilt $(z^{k})^{T}B_{k}z^{k}\geq 0$ und folgt demnach in diesem Fall $\nu _{k}>0$. Da $\varrho _{k}\neq 0$ ist und der Nenner in (7.16), wie die Herleitung von (7.15) zeigt, gerade dem Produkt

$\left(B^{k}y^{k}-{\frac {1}{\varrho _{k}}}s^{k}\right)^{T}y^{k}$

entspricht, ist demnach in diesem Fall auch (7.7) erfüllt.

q.e.d.

Die durch (7.14) gegebene Rang-1-Update-Formel

$B_{k+1}:={\frac {t_{k}}{\eta _{k}}}\left\{B_{k}-{\frac {\left(B_{k}z^{k}\right)\left(B_{k}z^{k}\right)^{T}}{\left(B_{k}z^{k}\right)^{T}y^{k}}}\right\}$ mit $z^{k}:=y^{k}-\eta _{k}g^{k}$

und mit einem frei wählbaren $\eta _{k}\in (0,1)$ wurde 1981 von Kleinmichel angegeben, wobei Kleinmichel aufgrund seiner numerischen Experimente die Wahl $\eta _{k}:=0.5$ für alle $k$ vorschlägt. Das dadurch bestimmte Kleinmichel-Verfahren funktioniert im Allgemeinen gut (s. [GeiKa99]). Für eine Restart-Version konnte Kleinmichel die $n$-Schritt quadratische Konvergenz beweisen.

7.3 Rang-2-Update-Formeln

Nach ca. 1970 wurden zahlreiche Rang-2-Update-Formeln für Quasi-Newton-Verfahren vorgeschlagen. Die wenigsten dieser Formeln haben sich in der Praxis durchgesetzt. Wir wollen hier nur die ein-parametrische Broyden-Klasse von Update-Formeln untersuchen. Sie ist die für die Praxis wichtigste Unterklasse der zwei-parametrischen Oren-Luenberger-Klasse von Update-Formeln und enthält die sog. BFGS-Formel, die sich unter den Rang-2-Formeln als die im Allgemeinen numerisch effizienteste herausgestellt hat.

Für die Rang-2-Update-Formeln wollen wir nicht so axiomatisch vorgehen wie im Fall der Rang-1-Formeln, obwohl das möglich wäre, sondern wir geben hier direkt die durch einen Parameter $\theta _{k}\geq 0$ bestimmte Broyden-Klasse von Update-Formeln an, welche durch

(7.17) $B_{k+1}:=B_{k}+\left(1+\theta _{k}{\frac {\tau _{k}}{\sigma _{k}}}\right){\frac {s^{k}(s^{k})^{T}}{\sigma _{k}}}-(1-\theta _{k}){\frac {B_{k}y^{k}(B_{k}y^{k})^{T}}{\tau _{k}}}-{\frac {\theta _{k}}{\sigma _{k}}}\left[s^{k}\left(B_{k}y^{k}\right)^{T}+\left(B_{k}y^{k}\right)(s^{k})^{T}\right]$

definiert ist, wobei

$s^{k}:=x^{k+1}-x^{k},\quad y^{k}:=g^{k+1}-g^{k}$

und

(7.18) $\sigma _{k}:=(s^{k})^{T}y^{k},\quad \tau _{k}:=(y^{k})^{T}B_{k}y^{k}$

seien. Eine spezielle Update-Formel der Broyden-Klasse erhält man also durch Festlegung von $\theta _{k}$, wobei $\theta _{k}$ typischerweise für alle $k$ konstant gewählt wird. (Auch andere Schreibweisen der Formel in (7.17) sind gebräuchlich.)

Ein Quasi-Newton-Verfahren vom Typ des Modellalgorithmus 7.1 mit einer Update-Formel (7.17) für $\theta _{k}\geq 0$ bezeichnen wir als Quasi-Newton-Verfahren der Broyden-Klasse. Damit ein solches Verfahren überhaupt durchführbar ist, ist sicherzustellen, dass $\sigma _{k}\neq 0$ für alle $k$ gilt. Da mit $\sigma _{k}\neq 0$ auch $y^{k}\neq 0$ ist, folgt mit der positiven Definitheit von $B_{k}$ dann auch $\tau _{k}>0$ für alle $k$.

Wir wollen nun zunächst diskutieren, unter welchen Voraussetzungen die Forderungen (F1) - (F3) für die Formel in (7.17) erfüllt sind. Dazu setzen wir insbesondere voraus:

(7.19) $B_{i}$ $(i=0,\ldots ,k)$ ist eine symmetrische, positiv definite Matrix,

wobei dies für $B_{0}$ bereits im Modellalgorithmus 7.1 vorausgesetzt wurde. Somit ist $p^{i}$ $(i=0,\ldots ,k)$ eine Abstiegsrichtung für $f$ in $x^{i}$ und damit insbesondere $p^{i}\neq 0$.

Sind die Vektoren $s^{k}$ und $B_{k}y^{k}$ linear unabhängig und gilt $\sigma _{k}\neq 0$, so handelt es sich bei (7.17) offenbar um eine Rang-2-Update-Formel mit $\varrho _{k}:=1$ und ist somit (F1) erfüllt; im anderen Fall reduziert sich die Formel auf eine Rang-1-Update-Formel. (Übung!) Ferner gilt für $B_{k+1}$ aus (7.17) die Quasi-Newton-Gleichung, wie das folgende Lemma besagt.

Lemma 7.7

Es sei $f\in C^{1}(\mathbb {R} ^{n})$ sowie $\sigma _{k}\neq 0$ und es sei (7.19) erfüllt. Dann genügt $B_{k+1}$ aus (7.17) der Quasi-Newton-Gleichung (7.3).

Beweis.

Da nach Voraussetzung $\sigma _{k}\neq 0$ und $B_{k}$ positiv definit ist, ist auch $\tau _{k}\neq 0$. Somit ist die Update-Formel in (7.17) wohldefiniert. Weiter ist $B_{k+1}y^{k}$ identisch mit

$B_{k}y^{k}+\left(1+\theta _{k}{\frac {\tau _{k}}{\sigma _{k}}}\right){\frac {s^{k}(s^{k})^{T}y^{k}}{\sigma _{k}}}-(1-\theta _{k}){\frac {B_{k}y^{k}(B_{k}y^{k})^{T}y^{k}}{\tau _{k}}}-{\frac {\theta _{k}}{\sigma _{k}}}\left[s^{k}\left(B_{k}y^{k}\right)^{T}y^{k}+\left(B_{k}y^{k}\right)(s^{k})^{T}y^{k}\right].$

Verwendung der Definitionen von $\sigma _{k}$ und $\tau _{k}$ aus (7.18) liefert

$B_{k+1}y^{k}=B_{k}y^{k}+\left(1+\theta _{k}{\frac {\tau _{k}}{\sigma _{k}}}\right)s^{k}-(1-\theta _{k})B_{k}y^{k}-\theta _{k}{\frac {\tau _{k}}{\sigma _{k}}}s^{k}-\theta _{k}B_{k}y^{k}=s^{k}.$

q.e.d.

Bleibt noch neben $\sigma _{k}\neq 0$ zu garantieren, dass unter der Voraussetzung in (7.19) für $B_{k+1}$ in (7.17) auch die Bedingung (F3) gesichert ist. Weil Matrizen vom Typ $uu^{T}$ und $[uv^{T}+vu^{T}]$ mit $u,v\in \mathbb {R} ^{n}$ symmetrisch sind, ist mit $B_{k}$ offenbar auch $B_{k+1}$ symmetrisch. Im Fall, dass $B_{k+1}$ positiv definit ist, ist dies auch die Matrix $B_{k+1}^{-1}$. Da weiter $B_{k+1}$ nach Lemma 7.7 der Quasi-Newton-Gleichung (7.3) genügt und da $p^{k}\neq 0$ und somit $s^{k}\neq 0$ gilt, folgt dann

$\sigma _{k}=(s^{k})^{T}y^{k}=(s^{k})^{T}B_{k+1}^{-1}s^{k}>0.$

Das folgende Lemma sagt nun aus, dass die Bedingung $\sigma _{k}>0$ auch hinreichend dafür ist, dass $B_{k+1}$ in (7.17) positiv definit ist, wobei wir den Beweis dafür hier nur skizzieren wollen. (Für das genauer diskutierte BFGS-Verfahren werden wir die positive Definitheit von $B_{k+1}$ weiter unten direkt nachweisen.)

Lemma 7.8

Es sei $f\in C^{1}(\mathbb {R} ^{n})$ und es gelte (7.19). Ist $\sigma _{k}>0$, so ist $B_{k+1}$ in (7.17) symmetrisch und positiv definit.

Beweis.

Dass $B_{k+1}$ symmetrisch ist, hatten wir schon erwähnt. Für den Nachweis der positiven Definitheit setze man

$u^{k}:=B_{k}^{-1/2}s^{k},\quad v^{k}:=B_{k}^{1/2}y^{k}.$

Dann gilt $B_{k+1}=B_{k}^{1/2}{\tilde {K}}_{k}B_{k}^{1/2}$ für

${\tilde {K}}_{k}:=I+(1+\theta _{k}{\frac {\tau _{k}}{\sigma _{k}}}){\frac {u^{k}(u^{k})^{T}}{\sigma _{k}}}-(1-\theta _{k}){\frac {v^{k}(v^{k})^{T}}{\tau _{k}}}-{\frac {\theta _{k}}{\sigma _{k}}}(u^{k}(v^{k})^{T}+v^{k}(u^{k})^{T}).$

Die Matrix $B_{k+1}$ ist demnach genau dann positiv definit, wenn ${\tilde {K}}_{k}$ dies ist. Letzteres ist der Fall, wenn alle Eigenwerte von ${\tilde {K}}_{k}$ positiv sind.

Sind $u^{k}$ und $v^{k}$ linear unabhängig, so besitzt ${\tilde {K}}_{k}$, wie man zeigen kann, den $(n-2)$-fachen Eigenwert 1 und zwei Eigenwerte mit Eigenvektoren aus dem von $u^{k}$ und $v^{k}$ erzeugten Teilraum des $\mathbb {R} ^{n}$. Für diese Eigenvektoren mache man den Ansatz $x=\alpha u^{k}+\beta v^{k}$ und rechne man die beiden unbekannten Eigenwerte aus. Die gemachten Voraussetzungen garantieren dann, dass diese Eigenwerte positiv sind. Abschließend diskutiere man den einfacheren Fall, dass $u^{k}$ und $v^{k}$ linear abhängig sind.

q.e.d.

Insofern sind Voraussetzungen von Interesse, unter denen die Bedingung $\sigma _{k}>0$ erfüllt ist.

Lemma 7.9

Es seien (V1) - (V4) sowie die Voraussetzung in (7.19) erfüllt und es sei $t_{k}>0$. Dann folgt $\sigma _{k}>0$.

Beweis.

Wie wir oben festgestellt haben, ist unter der Voraussetzung in (7.19) $p^{i}$ $(i=0,\ldots ,k)$ eine Abstiegsrichtung für $f$ in $x^{i}$, so dass gemäß Lemma 2.13 $x^{i}\in N_{0}$ $(i=0,\ldots ,k+1)$ für die Menge $N_{0}$ aus (2.9) folgt. Wegen $p^{k}\neq 0$ und $t_{k}>0$ kann man daher insbesondere mit Lemma 2.9 (ii) schließen:

$\sigma _{k}=(g^{k+1}-g^{k})^{T}(x^{k+1}-x^{k})\geq \beta \left\|x^{k+1}-x^{k}\right\|^{2}=\beta t_{k}^{2}\left\|p^{k}\right\|^{2}>0.$

q.e.d.

Für gleichmäßig konvexe Funktionen ist also für alle $k$ die Bedingung $\sigma _{k}>0$ und damit die Durchführbarkeit eines Quasi-Newton-Verfahrens der Broyden-Klasse gesichert. Ist jedoch $f$ keine gleichmäßig konvexe Funktion, so kann die Bedingung $\sigma _{k}>0$ für alle Verfahren der Broyden-Klasse gemeinsam nur im Fall der Verwendung bestimmter Schrittweitenregeln garantiert werden. Und zwar gilt (Übung!):

Lemma 7.10

Es sei $f\in C^{1}(\mathbb {R} ^{n})$ und es gelte (7.19). Ist $t_{k}>0$ eine exakte Schrittweite, eine Wolfe-Powell- oder eine strenge Wolfe-Powell-Schrittweite, so folgt $\sigma _{k}>0$.

Im Fall, dass die Matrix $B_{k+1}$ in (7.17) für das gewählte $\theta _{k}$ positiv definit ist, kann ihre Inverse mit

$\varepsilon _{k}:=(s^{k})^{T}B_{k}^{-1}s^{k}$

explizit angegeben werden in der Form

$B_{k+1}^{-1}=B_{k}^{-1}+\left(1+{\frac {\varepsilon _{k}\sigma _{k}(1-\theta _{k})}{(\sigma _{k}^{2}+\theta _{k}(\varepsilon _{k}\tau _{k}-\sigma _{k}^{2}))}}\right){\frac {y^{k}(y^{k})^{T}}{\sigma _{k}}}-{\frac {\tau _{k}\theta _{k}}{(\sigma _{k}^{2}+\theta _{k}(\varepsilon _{k}\tau _{k}-\sigma _{k}^{2}))}}B_{k}^{-1}s^{k}(B_{k}^{-1}s^{k})^{T}$

(7.20) $-{\frac {(1-\theta _{k})\sigma _{k}}{(\sigma _{k}^{2}+\theta _{k}(\varepsilon _{k}\tau _{k}-\sigma _{k}^{2}))}}\left[y^{k}(B_{k}^{-1}s^{k})^{T}+B_{k}^{-1}s^{k}(y^{k})^{T}\right].$

Man prüft leicht, aber mit einiger Schreibarbeit nach, dass tatsächlich $B_{k+1}B_{k+1}^{-1}=I$ gilt. Diese Formel ist offenbar eine Update-Formel für $B_{k+1}^{-1}$. Ist $B_{k+1}$ positiv definit, so ist auch $B_{k+1}^{-1}$ positiv definit.

Bemerkung 7.11

In unserer Formulierung von Quasi-Newton-Verfahren (siehe den Modellalgorithmus 7.1) wird $p^{k}$ als das Matrix-Vektor-Produkt „$-B_{k}g^{k}$“ bestimmt. Die Berechnung dieses Produktes erfordert ${\mathcal {O}}(n^{2})$ arithmetische Rechenoperationen. (Der Ausdruck ${\mathcal {O}}(n^{p})$ steht für $Cn^{p}$, wobei $C>0$ eine von $n$ unabhängige Konstante ist.) Im Hinblick auf das Newton-Verfahren kann man $p^{k}$ mit $H_{k}:=B_{k}^{-1}$ auch durch $H_{k}p^{k}=-g^{k}$ definieren und $H_{k+1}$ aus $H_{k}$ gemäß (7.20) berechnen. Ähnlich wie beim Newton-Verfahren wird $p^{k}$ bei einer solchen Vorgehensweise dann auch durch Lösung des linearen Gleichungssystems $H_{k}p^{k}=-g^{k}$ gewonnen.

Ein derartiges Vorgehen wäre nicht sinnvoll, wenn man dann in jeder Iteration eine vollständige Cholesky-Zerlegung der Matrix $H_{k+1}$ mit dem aus der Numerischen Mathematik bekannten Verfahren vornehmen würde, da eine solche Zerlegung ${\mathcal {O}}(n^{3})$ arithmetische Rechenoperationen erfordert. Es ist jedoch möglich, aus einer gegebenen Cholesky-Zerlegung $H_{k}=L_{k}L_{k}^{T}$ (die Startmatrix $H_{0}$ wird oft als Diagonalmatrix gewählt) durch ein Rang-1-Update von $L_{k}$ in ${\mathcal {O}}(n^{2})$ Rechenoperationen zu einer Cholesky-Zerlegung für $H_{k+1}$ zu gelangen (z. B. [Wer92], [GeiKa99]).

Der Gewinn dabei ist, dass man dann anders als bei der zuerst beschriebenen Vorgehensweise quasi umsonst feststellen kann, ob $H_{k+1}$ tatsächlich positiv definit und $p^{k}$ damit eine Abstiegsrichtung für $f$ in $x^{k}$ ist (vgl. (7.2)). Denn in der Praxis kann man ja z. B. exakte Schrittweiten nicht genau berechnen oder kann man im Allgemeinen nicht verifizieren, ob man sich tatsächlich in einem Bereich von $f$ befindet, in dem $f$ gleichmäßig konvex ist (s. die Voraussetzungen von Lemma 7.9 und Satz 7.14).

In der nachfolgenden Bemerkung fassen wir außerdem ohne Beweis einige verblüffende Ergebnisse zu Quasi-Newton-Verfahren der Broyden-Klasse zusammen. Beweise dafür sind z.B. in [Fle91] und zum Teil auch in [Wer92] und [JaSt04] zu finden.

Bemerkung 7.12

Der Modellalgorithmus 7.1 sei mit einer exakten Schrittweitenregel (!) und einer Update-Formel der Broyden-Klasse in (7.17) versehen. (Mit den Lemmata 7.10 und 7.8 schließt man induktiv, dass dann $B_{k}$ für alle $k$ symmetrisch und positiv definit ist.)

(i) Ist $f$ eine quadratische Funktion, d. h. ist

$f(x):={\frac {1}{2}}x^{T}Qx+c^{T}x+\alpha ,\quad x\in \mathbb {R} ^{n}$

und ist $Q$ positiv definit, so bricht jedes solche (durch die Wahl der $\theta _{k}$ bestimmte) Quasi-Newton-Verfahren nach $m\leq n$ Schritten ab und sind die von ihm erzeugten Richtungen $p^{0},\ldots ,p^{m-1}$ $Q$-konjugiert. Wird insbesondere $B_{0}:=I$ gewählt, so erzeugt jedes derartige Verfahren unabhängig von der Wahl der $\theta _{k}\geq 0$ bei gleichem Startwert $x^{0}$ dieselben Iterierten wie das Fletcher-Reeves-Verfahren bzw. wie das im quadratischen Fall ja damit identische Polak-Ribière-Verfahren.

(ii) Für jede, also auch jede nichtquadratische Funktion $f\in C^{1}(\mathbb {R} ^{n})$ erzeugen Quasi-Newton-Verfahren des oben spezifizierten Typs bei gleichen Startwerten für alle $\theta _{k}\geq 0$ dieselben Iterierten $x^{k}$.

Aus der allgemeinen Update-Formel der Broyden-Klasse in (7.17) erhält man insbesondere auch die Broyden-Rang-1- bzw. SR1-Formel (7.11), wenn man

$\theta _{k}:={\frac {\sigma _{k}}{\sigma _{k}-\tau _{k}}}$

setzt, wobei $\sigma _{k}-\tau _{k}\neq 0$ nicht gesichert ist. Weiter leitet man für $\theta _{k}:=0$ aus (7.17) die sog. DFP-Formel

(7.21) $B_{k+1}:=B_{k}+{\frac {s^{k}(s^{k})^{T}}{\sigma _{k}}}-{\frac {B_{k}y^{k}(B_{k}y^{k})^{T}}{\tau _{k}}}$

ab, welche nach Davidon (1959) und Fletcher und Powell (1963) benannt wird, die diese Formel unabhängig voneinander angegeben hatten.

Schließlich gewinnt man aus (7.17) für $\theta _{k}:=1$ die wichtigste Update-Formel für Quasi-Newton-Verfahren, die BFGS-Formel

(7.22) $B_{k+1}:=B_{k}+\left(1+{\frac {\tau _{k}}{\sigma _{k}}}\right){\frac {s^{k}(s^{k})^{T}}{\sigma _{k}}}-{\frac {1}{\sigma _{k}}}\left[s^{k}\left(B_{k}y^{k}\right)^{T}+\left(B_{k}y^{k}\right)(s^{k})^{T}\right]$

Diese Formel bezieht ihren Namen aus den Anfangsbuchstaben der Namen von Broyden, Fletcher, Goldfarb und Shanno, die sie im Jahre 1970 unabhängig voneinander und nahezu zeitgleich mit unterschiedlichen Begründungen vorgeschlagen hatten. Die Quasi-Newton-Verfahren vom Typ des Modellalgorithmus 7.1 mit der DFP- und der BFGS-Formel heißen DFP- bzw. BFGS-Verfahren. Man kann leicht nachprüfen, dass die gesamte Broyden-Klasse von Update-Formeln in (7.17) in der Form

$B_{k+1}=(1-\theta _{k})B_{k+1}^{DFP}+\theta _{k}B_{k+1}^{BFGS}\quad (\theta _{k}\geq 0)$

beschrieben werden kann, wobei $B_{k+1}^{DFP}$ und $B_{k+1}^{BFGS}$ die Matrizen $B_{k+1}$ in (7.21) und (7.22) bezeichnen mögen.

In der Praxis hat sich das BFGS-Verfahren als das verlässlichste und effizienteste unter den Rang-2-Quasi-Newton-Verfahren durchgesetzt. Im nächsten Abschnitt wollen wir daher die Konvergenz des BFGS-Verfahrens für gleichmäßig konvexe Funktionen beweisen, was sogar für beliebige semieffiziente Schrittweitenregeln gelingt. Man beachte, dass damit nach Bemerkung 7.12 bei Verwendung einer exakten Schrittweitenregel für gleichmäßig konvexe Funktionen auch die Konvergenz jedes anderen Verfahrens der Broyden-Klasse, also insbesondere die des DFP-Verfahrens bewiesen ist.

Es hat sich jedoch herausgestellt und numerische Beispiele in der angegebenen Literatur veranschaulichen dies, dass das DFP-Verfahren sehr viel empfindlicher auf die Wahl der Schrittweiten als das BFGS-Verfahren reagiert (z.B. [Fle91], [GeiKa99]). So konnte für das DFP-Verfahren im Hinblick auf die Schrittweitenregel auch kein so allgemeines Ergebnis wie das aus dem nächsten Abschnitt für das BFGS-Verfahren bewiesen werden. Hinweise darauf, warum das BFGS-Verfahren insgesamt bessere numerische Ergebnisse als z. B. das DFP-Verfahren liefert, findet man in der neueren Literatur (z. B. [NoWri06], [JaSt04, S. 180]).

7.4 Das BFGS-Verfahren

Das BFGS-Verfahren, versehen mit einer semieffizienten Schrittweitenregel, lautet:

Algorithmus 7.13 (BFGS-Verfahren)

(0) Wähle $x^{0}\in \mathbb {R} ^{n}$, eine symmetrische, positiv definite Matrix $B_{0}\in \mathbb {R} ^{n\times n}$ und eine semieffiziente Schrittweitenregel. Setze $k:=0$.

(1) Falls $g^{k}=0$ ist, stop! ($x^{k}$ ist kritische Lösung von Problem $(P)$.)

(2) Berechne

$p^{k}:=-B_{k}g^{k},$

eine Schrittweite $t_{k}>0$ und setze

$x^{k+1}:=x^{k}+t_{k}p^{k}.$

(3) Berechne

${\begin{array}{ll}s^{k}:=x^{k+1}-x^{k},&y^{k}:=g^{k+1}-g^{k},\\\sigma _{k}:=(s^{k})^{T}y^{k},&\tau _{k}:=(y^{k})^{T}B_{k}y^{k}\end{array}}$

und

(7.23) $B_{k+1}:=B_{k}+\left(1+{\frac {\tau _{k}}{\sigma _{k}}}\right){\frac {s^{k}(s^{k})^{T}}{\sigma _{k}}}-{\frac {1}{\sigma _{k}}}\left[s^{k}(B_{k}y^{k})^{T}+B_{k}y^{k}(s^{k})^{T}\right].$

(4) Setze $k:=k+1$ und gehe nach (1).

Wir wollen zunächst beweisen, dass die Matrix $B_{k+1}$ im Fall der BFGS-Formel für eine gleichmäßig konvexe Zielfunktion $f$ für jedes $k$ positiv definit ist. (Den Beweis der positiven Definitheit von $B_{k+1}$ für die ganze Broyden-Klasse hatten wir ja nicht vollständig ausgeführt.) Unter den entsprechenden Voraussetzungen fällt Algorithmus 7.13 also in das Schema des Modellalgorithmus 2.5.

Satz 7.14

Es seien (V1) - (V4) erfüllt, und es sei $B_{i}$ $(i=0,1,\ldots ,k)$ symmetrisch und positiv definit. In der $k$-ten Iteration von Algorithmus 7.13 hat man dann $\sigma _{k}>0$ und $\tau _{k}>0$ und ist auch $B_{k+1}$ in (7.23) positiv definit. Mit $\varepsilon _{k}:=(s^{k})^{T}B_{k}^{-1}s^{k}$ folgt ferner

(7.24) $B_{k+1}^{-1}=B_{k}^{-1}+{\frac {y^{k}(y^{k})^{T}}{\sigma _{k}}}-{\frac {B_{k}^{-1}s^{k}(B_{k}^{-1}s^{k})^{T}}{\varepsilon _{k}}}.$

Beweis.

Lemma 7.9 impliziert unter den gegebenen Voraussetzungen $\sigma _{k}>0$. Damit hat man $y^{k}\neq 0$ und demnach $\tau _{k}>0$.

Die Matrix $B_{k+1}hatnunmitZahlen<math>\alpha ,\beta \in \mathbb {R}$ und Vektoren $a,b\in \mathbb {R} ^{n}$ die Gestalt

$B_{k+1}=B_{k}+\alpha aa^{T}-\beta \left[ab^{T}+ba^{T}\right].$

Somit folgt die Symmetrie von $B_{k+1}$ wegen

$B_{k+1}^{T}=B_{k}^{T}+\alpha aa^{T}-\beta \left[ba^{T}+ab^{T}\right]=B_{k+1}.$

Bleibt zu zeigen, dass $z^{T}B_{k+1}z>0$ für alle $z\neq 0$ gilt.

Mit

$u:=B_{k}^{1/2}z,\quad v:=B_{k}^{1/2}y^{k}$

erhalten wir wegen $\tau _{k}=\|v\|^{2}$ und $v^{T}u=(B_{k}y^{k})^{T}z$

$z^{T}B_{k+1}z=z^{T}B_{k}z+\left(1+{\frac {\tau _{k}}{\sigma _{k}}}\right){\frac {(z^{T}s^{k})^{2}}{\sigma _{k}}}-{\frac {1}{\sigma _{k}}}z^{T}\left[s^{k}(B_{k}y^{k})^{T}+B_{k}y^{k}(s^{k})^{T}\right]z$

(7.25) $=\|u\|^{2}+{\frac {1}{\sigma _{k}}}(z^{T}s^{k})^{2}+{\frac {\|v\|^{2}}{\sigma _{k}^{2}}}(z^{T}s^{k})^{2}-{\frac {2}{\sigma _{k}}}(v^{T}u)(z^{T}s^{k}).$

Ist $(v^{T}u)(z^{T}s^{k})\leq 0$, so folgt weiter $z^{T}B_{k+1}z\geq \|u\|^{2}>0$ für $u\neq 0$ bzw. $z\neq 0$, so dass in diesem Fall alles gezeigt ist. Gilt andererseits $(v^{T}u)(z^{T}s^{k})>0$ und damit insbesondere $z^{T}s^{k}\neq 0$, so kann man folgendermaßen abschätzen:

$z^{T}B_{k+1}z=\|u\|^{2}+{\frac {1}{\sigma _{k}}}(z^{T}s^{k})^{2}+{\frac {\|v\|^{2}}{\sigma _{k}^{2}}}(z^{T}s^{k})^{2}-{\frac {2}{\sigma _{k}}}\left|v^{T}u\right|\left|z^{T}s^{k}\right|$

$=\left({\frac {\|v\|}{\sigma _{k}}}\left|z^{T}s^{k}\right|-\|u\|\right)^{2}+{\frac {1}{\sigma _{k}}}(z^{T}s^{k})^{2}+{\frac {2}{\sigma _{k}}}\left(\|u\|\|v\|-\left|v^{T}u\right|\right)\left|z^{T}s^{k}\right|\geq {\frac {1}{\sigma _{k}}}(z^{T}s^{k})^{2}>0.$

Also ist $B_{k+1}$ positiv definit. Dass die in (7.24) angegebene Matrix die Inverse von $B_{k+1}$ ist, zeigt man durch Verifizierung der Gleichung $B_{k+1}B_{k+1}^{-1}=I$, was etwas aufwändig, aber nicht schwierig ist.

q.e.d.

Der erste Nachweis der Konvergenz des BFGS-Verfahrens wurde von Powell erbracht und stellt eine der großen Leistungen in der Optimierung dar. Da das BFGS-Verfahren eines der wichtigsten Optimierungsverfahren ist, wollen wir hier auch einen Konvergenzbeweis, der auf Werner zurückgeht, für gleichmäßig konvexe Funktionen und eine beliebige semieffiziente Schrittweitenregel führen. Dieser Beweis ist aber recht lang und technisch, so dass der eher numerisch orientierte Leser ihn überspringen mag. Wir benötigen dazu:

Lemma 7.15

(i) Es seien $A,B\in \mathbb {R} ^{n\times n}$. Dann gilt

$\operatorname {Spur} (AB)=\operatorname {Spur} (BA).$

(ii) Es sei $A\in \mathbb {R} ^{n\times n}$ eine symmetrische Matrix und $\lambda _{i}(A)$ $(i=1,\ldots ,n)$ seien ihre Eigenwerte. Dann gilt

$\operatorname {Spur} (A)=\sum _{i=1}^{n}\lambda _{i}(A).$

Damit können wir nun die Konvergenz des BFGS-Verfahrens für gleichmäßig konvexe Funktionen beweisen.

Satz 7.16

Es seien (V1) - (V4) erfüllt. Algorithmus 7.13 bricht entweder nach endlich vielen Schritten mit der Lösung $x^{*}$ von $(P)$ ab oder er liefert eine Folge $\left\{x^{k}\right\}$, die gegen $x^{*}$ konvergiert. Genauer gilt sogar

$\sum _{k=0}^{\infty }\left\|x^{k}-x^{*}\right\|<\infty .$

Beweis.

Wir nehmen an, dass Algorithmus 7.13 nicht nach endlich vielen Schritten abbricht. Aus Satz 7.14 schließt man dann induktiv für alle $k$, dass $\sigma _{k}>0,\tau _{k}>0,s^{k}\neq 0$ und $y^{k}\neq 0$ gilt und dass $B_{k}$ positiv definit ist. Für alle $k$ ist damit ferner $p^{k}$ eine Abstiegsrichtung (vgl. (7.2)) und gemäß Definition 2.12 einer semieffizienten Schrittweitenregel mit $x^{k-1}$ auch $x^{k}\in N_{0}$ für die Menge $N_{0}$ aus (2.9). Sei nun $k\in \mathbb {N}$ beliebig.

Es wird eine semieffiziente Schrittweitenregel verwendet, so dass gemäß Definition 2.12 eine Konstante $\vartheta >0$ existiert mit

(7.26) $f(x^{k})-f(x^{k+1})\geq \vartheta \min \left\{\left({\frac {(g^{k})^{T}p^{k}}{\|p^{k}\|}}\right)^{2},\left|(g^{k})^{T}p^{k}\right|\right\}=\vartheta \delta _{k}\left\|g^{k}\right\|^{2}$

für

$\delta _{k}:=\min \left\{\left({\frac {(g^{k})^{T}p^{k}}{\|g^{k}\|\|p^{k}\|}}\right)^{2},{\frac {\left|(g^{k})^{T}p^{k}\right|}{\|g^{k}\|^{2}}}\right\}$

Nach Aussage (v) von Lemma 2.9 hat man

$\left\|g^{k}\right\|^{2}\geq 2\beta (f(x^{k})-f(x^{*})).$

Folglich gilt mit (7.26)

$f(x^{k})-f(x^{k+1})\geq 2\beta \vartheta \delta _{k}\left(f(x^{k})-f(x^{*})\right).$

Unter Verwendung der Beziehung $1+x\leq e^{x}$ für $x\in \mathbb {R}$ gewinnt man daraus

$0\leq f(x^{k+1})-f(x^{*})\leq (1-2\beta \vartheta \delta _{k})\left(f(x^{k})-f(x^{*})\right)$

$\leq \prod _{j=0}^{k}(1-2\beta \vartheta \delta _{j})\left(f(x^{0})-f(x^{*})\right)\leq \exp(-2\beta \vartheta \sum _{j=0}^{k}\delta _{j})\left(f(x^{0})-f(x^{*})\right).$

Wenn wir nun zeigen können, dass ein $\delta >0$ existiert mit

(7.27) $\delta \leq {\frac {1}{k+1}}\sum _{j=0}^{k}\delta _{j},$

so folgte aus Letzterem

$0\leq f(x^{k+1})-f(x^{*})\leq \exp(-2\beta \vartheta \delta (k+1))\left(f(x^{0})-f(x^{*})\right)$

Da nach Lemma 2.9 (iii)

${\frac {\beta }{2}}\left\|x^{k}-x^{*}\right\|^{2}\leq f(x^{k})-f(x^{*})$

gilt, hätte man dann weiter

$\left\|x^{k}-x^{*}\right\|\leq \exp(-\beta \vartheta \delta k)\left({\frac {2(f(x^{0})-f(x^{*}))}{\beta }}\right)^{1/2}.$

Summation über alle k lieferte schließlich

(7.28) $\sum _{k=0}^{\infty }\left\|x^{k}-x^{*}\right\|=\left({\frac {2(f(x^{0})-f(x^{*}))}{\beta }}\right)^{1/2}\sum _{k=0}^{\infty }\exp(-\beta \vartheta \delta k)<\infty ,$

wobei die Konvergenz der rechten Reihe wegen $\beta >0,\vartheta >0$ und $\delta >0$ leicht mit dem Quotientenkriterium erschlossen werden kann. Da die Konvergenz einer Reihe nach sich zieht, dass ihre Glieder gegen 0 streben, folgte aus (7.28) auch die gewünschte Konvergenz $x^{k}\to x^{*}$ $(k\to \infty )$.

Wir wollen also (7.27) für ein $\delta >0$ nachweisen. Dazu zeigen wir, dass für Konstanten $c_{1},c_{2}>0$ die Ungleichungen

(7.29) $\sum _{j=0}^{k}{\frac {\|g^{j}\|^{2}}{|(g^{j})^{T}p^{j}|}}\leq c_{1}(k+1)$

und

(7.30) $prod_{j=0}^{k}\left\{{\frac {\|g^{j}\|\|p^{j}\|}{|(g^{j})^{T}p^{j}|}}\right\}^{2}\leq c_{2}^{k+1}$

gelten. Denn hat man dies, so kann man folgenden, auf Powell zurückgehenden „Trick“ anwenden ($|I_{k}|$ bezeichnet dabei die Anzahl der Elemente von $I_{k}$).

Es gilt: Sind $\gamma _{j}\geq 0$ und $c>0$ Zahlen mit $\sum _{j=0}^{k}\gamma _{j}\leq (k+1)c$, dann hat man für

$I_{k}:=\{j\in \{0,\ldots ,k\}{\big |}\gamma _{j}\leq 3c\}$

die Abschätzung $|I_{k}|\geq {\frac {2}{3}}(k+1)$.

Setzt man ${\bar {I}}_{k}:=\{0,\ldots ,k\}\setminus I_{k}$, so folgt dies für ${\bar {I}}_{k}=\emptyset$ sofort und für ${\bar {I}}_{k}\neq \emptyset$ aus

$(k+1)c\geq \sum _{j=0}^{k}\gamma _{j}\geq \sum _{j\in {\bar {I}}_{k}}\gamma _{j}>3c\left|{\bar {I}}_{k}\right|=3c(k+1-|I_{k}|).$

Ferner gilt: Sind $\delta _{j}\geq 1$ und $d>0$ Zahlen mit $\prod _{j=0}^{k}\delta _{j}\leq d^{k+1}$, dann hat man für

$J_{k}:=\left\{j\in \{0,\ldots ,k\}{\big |}\delta _{j}\leq d^{3}\right\}$

die Abschätzung $|J_{k}|\geq {\frac {2}{3}}(k+1)$.

Letztere Aussage folgt offenbar aus der ersteren durch Logarithmisierung mit $\gamma _{j}:=\ln(\delta _{j})$ und $c:=\ln(d)$.

Gibt es nun von einer $(k+1)$-elementigen Menge eine Teilmenge $I_{k_{1}}$ mit $|I_{k_{1}}|\geq {\frac {2}{3}}(k+1)$ und einer „Eigenschaft 1“ sowie eine Teilmenge $I_{k_{2}}$ mit $|I_{k_{2}}|\geq {\frac {2}{3}}(k+1)$ und einer „Eigenschaft 2“, so existiert offenbar eine Teilmenge $I_{k}$ der Menge mit $|I_{k}|\geq {\frac {1}{3}}(k+1)$, welche beide Eigenschaften besitzt. Aus (7.29) und (7.30) folgert man daher bei festem $k$ die Existenz einer Teilmenge $I_{k}\subseteq \{0,\ldots ,k\}$ mit $|I_{k}|\geq {\frac {1}{3}}(k+1)$, d. h. insbesondere $|I_{k}|\neq 0$ und

${\frac {\|g^{j}\|^{2}}{|(g^{j})^{T}p^{j}|}}\leq 3c_{1},\quad \left\{{\frac {\|g^{j}\|\|p^{j}\|}{|(g^{j})^{T}p^{j}|}}\right\}^{2}\leq c_{2}^{3},\quad j\in I_{k}.$

Damit bekommt man weiter

$\sum _{j=0}^{k}\delta _{j}=\sum _{j=0}^{k}\min \left(\left\{{\frac {(g^{j})^{T}p^{j}}{\|g^{j}\|\|p^{j}\|}}\right\}^{2},{\frac {\left|(g^{j})^{T}p^{j}\right|}{\|g^{j}\|^{2}}}\right)\geq \sum _{j\in I_{k}}\min \left(\left\{{\frac {(g^{j})^{T}p^{j}}{\|g^{j}\|\|p^{j}\|}}\right\}^{2},{\frac {\left|(g^{j})^{T}p^{j}\right|}{\|g^{j}\|^{2}}}\right)$

$\geq {\frac {1}{3}}(k+1)\min \left(c_{2}^{-3},{\frac {1}{3c_{1}}}\right),$

so dass man für $\delta :=\min(c_{2}^{-3}/3,1/(9c_{1}))$ die Ungleichung (7.27) erschließt. Folglich müssen wir noch (7.29) und (7.30) beweisen.

Zunächst wollen wir die Gültigkeit von (7.29) nachweisen. Mit $\varepsilon _{k}:=(s^{k})^{T}B_{k}^{-1}s^{k}$ ist

$B_{k+1}^{-1}=B_{k}^{-1}+{\frac {y^{k}(y^{k})^{T}}{\sigma _{k}}}-{\frac {B_{k}^{-1}s^{k}(B_{k}^{-1}s^{k})^{T}}{\varepsilon _{k}}}$

und daher unter Verwendung der Lemmata (7.5) und (7.15)

(7.31) $0<\operatorname {Spur} (B_{k+1}^{-1})=\operatorname {Spur} (B_{k}^{-1})+{\frac {\left\|y^{k}\right\|^{2}}{\sigma _{k}}}-{\frac {\left\|B_{k}^{-1}s^{k}\right\|^{2}}{\varepsilon _{k}}}=\operatorname {Spur} \left(B_{0}^{-1}\right)+\sum _{j=0}^{k}{\frac {\left\|y^{j}\right\|^{2}}{\sigma _{j}}}-\sum _{j=0}^{k}{\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}}{\varepsilon _{j}}}.$

Nach Aussage (ii) von Lemma 2.9 hat man

(7.32) $\sigma _{j}=(y^{j})^{T}s^{j}\geq \beta \left\|s^{j}\right\|^{2}$

so dass mit (V3) folgt:

(7.33) ${\frac {\|y^{j}\|^{2}}{\sigma _{j}}}\leq {\frac {\gamma ^{2}\|s^{j}\|^{2}}{\beta \|s^{j}\|^{2}}}={\frac {\gamma ^{2}}{\beta }}$

Ferner gilt

(7.34) $B_{j}^{-1}s^{j}=t_{j}B_{j}^{-1}p^{j}=-t_{j}g^{j},$

(7.35) $\varepsilon _{j}=(s^{j})^{T}B_{j}^{-1}s^{j}=(t_{j}p^{j})^{T}(-t_{j}g^{j})=t_{j}^{2}\left|(g^{j})^{T}p^{j}\right|.$

Somit kann man mit

$c_{1}:=2\max \left(\operatorname {Spur} (B_{0}^{-1}),{\frac {\gamma ^{2}}{\beta }}\right)$

aus (7.31) schließen:

$\sum _{j=0}^{k}{\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}}{\varepsilon _{j}}}=\sum _{j=0}^{k}{\frac {\|g^{j}\|^{2}}{|(g^{j})^{T}p^{j}|}}\leq \operatorname {Spur} (B_{0}^{-1})+\sum _{j=0}^{k}{\frac {\|y^{j}\|^{2}}{\sigma _{j}}}$

(7.36) $\leq {\frac {1}{2}}c_{1}+{\frac {1}{2}}(k+1)c_{1}\leq c_{1}(k+1).$

Also gilt (7.29) und haben wir nur noch (7.30) zu zeigen. In diesem Zusammenhang halten wir aber noch die folgende Abschätzung fest, die sich aus (7.31) mit (7.33), (7.36) und der Definition von $c_{1}$ ergibt:

(7.37) $\operatorname {Spur} (B_{k+1}^{-1})\leq {\frac {c_{1}}{2}}+{\frac {k+1}{2}}c_{1}+c_{1}(k+1).$

Offenbar ist $B_{k+1}^{-1}=B_{k}^{-1/2}K_{k}B_{k}^{-1/2}$ mit

$K_{k}=I+{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}-{\frac {B_{k}^{-1/2}s^{k}(B_{k}^{-1/2}s^{k})^{T}}{\varepsilon _{k}}}.$

Für $v\neq 0$ gilt nun, wie man leicht überprüft (siehe auch Aufgabenblatt 5):

$\left(I+vv^{T}\right)^{-1}=I-{\frac {1}{1+\|v\|^{2}}}vv^{T},$

so dass speziell für $v:=(1/{\sqrt {\sigma _{k}}})B_{k}^{1/2}y^{k}$ folgt: $\|v\|^{2}=\tau _{k}/\sigma _{k}$ und somit

$\left(I+{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}\right)^{-1}=I-{\frac {\sigma _{k}}{\sigma _{k}+\tau _{k}}}{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}.$

Demnach kann man schreiben

$K_{k}=\left(I+{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}\right)\cdot \left\{I-\left(I-{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}\right){\frac {B_{k}^{-1/2}s^{k}(B_{k}^{-1/2}s^{k})^{T}}{\varepsilon _{k}}}\right\}$

$=\left(I+{\frac {B_{k}^{1/2}y^{k}(B_{k}^{1/2}y^{k})^{T}}{\sigma _{k}}}\right)\cdot \left\{I+\left({\frac {\sigma _{k}}{\varepsilon _{k}(\sigma _{k}+\tau _{k})}}B_{k}^{1/2}y^{k}-{\frac {1}{\varepsilon _{k}}}B_{k}^{-1/2}s^{k}\right)\left(B_{k}^{-1/2}s^{k}\right)^{T}\right\}.$

Berücksichtigt man nun Lemma 7.5, so bekommt man

$\det(K_{k})=\left(1+{\frac {\tau _{k}}{\sigma _{k}}}\right)\left(1+{\frac {\sigma _{k}^{2}}{\varepsilon _{k}(\sigma _{k}+\tau _{k})}}-{\frac {\varepsilon _{k}}{\varepsilon _{k}}}\right)={\frac {\sigma _{k}+\tau _{k}}{\sigma _{k}}}{\frac {\sigma _{k}^{2}}{\varepsilon _{k}(\sigma _{k}+\tau _{k})}}={\frac {\sigma _{k}}{\varepsilon _{k}}}.$

Damit erreicht man:

$\det(B_{k+1}^{-1})=\det(B_{k}^{-1/2})\det(K_{k})\det(B_{k}^{-1/2})=\det(K_{k})\det(B_{k}^{-1})$

(7.38) $={\frac {\sigma _{k}}{\varepsilon _{k}}}\det(B_{k}^{-1})=\det(B_{0}^{-1})\prod _{j=0}^{k}{\frac {\sigma _{j}}{\varepsilon _{j}}}$

Sind $\lambda _{i}(B_{k+1}^{-1})dieEigenwertevon<math>B_{k+1}^{-1}$, so schließt man außerdem mit der Ungleichung vom geometrisch-arithmetischen Mittel, mit Lemma 7.15 sowie mit der Abschätzung in (7.37)

$\det(B_{k+1}^{-1})=\prod _{i=1}^{n}\lambda _{i}(B_{k+1}^{-1})\leq \left\{{\frac {1}{n}}\sum _{i=1}^{n}\lambda _{i}(B_{k+1}^{-1})\right\}^{n}=\left\{{\frac {1}{n}}\operatorname {Spur} (B_{k+1}^{-1})\right\}^{n}\leq \left\{{\frac {1}{n}}\left[{\frac {1}{2}}c_{1}+{\frac {k+1}{2}}c_{1}+c_{1}(k+1)\right]\right\}^{n}$

(7.39) $\leq \left\{{\frac {2}{n}}c_{1}(k+1)\right\}^{n}.$

Die Ungleichung vom geometrisch-arithmetischen Mittel besagt, dass für Zahlen $a_{j}\geq 0$ und $n\in \mathbb {N}$ gilt:

${\sqrt[{n}]{a_{1}\cdot a_{2}\cdot \ldots \cdot a_{n}}}\leq {\frac {1}{n}}(a_{1}+a_{2}+\ldots +a_{n}).$

Zusammen implizieren (7.38) und (7.39)

(7.40) $\prod _{j=0}^{k}{\frac {\sigma _{j}}{\varepsilon _{j}}}\leq {\frac {1}{\det(B_{0}^{-1})}}\left\{{\frac {2}{n}}c_{1}(k+1)\right\}^{n}.$

Aus einer weiteren Anwendung der Ungleichung vom geometrisch-arithmetischen Mittel erhält man ferner mit (7.36)

(7.41) $\prod _{j=0}^{k}{\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}}{\varepsilon _{j}}}\leq \left\{{\frac {1}{k+1}}\sum _{j=0}^{k}{\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}}{\varepsilon _{j}}}\right\}^{k+1}\leq c_{1}^{k+1}.$

Unter Ausnutzung von (7.34), (7.35), (7.32), (7.41) und (7.40) gelangt man schließlich zu

$\prod _{j=0}^{k}{\frac {\|g^{j}\|^{2}\|p^{j}\|^{2}}{|(g^{j})^{T}p^{j}|^{2}}}=\prod _{j=0}^{k}{\frac {{\frac {1}{t_{j}^{2}}}\left\|B_{j}^{-1}s^{j}\right\|^{2}{\frac {1}{t_{j}^{2}}}\|s^{j}\|^{2}}{{\frac {1}{t_{j}^{4}}}\varepsilon _{j}^{2}}}=\prod _{j=0}^{k}{\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}\|s^{j}\|^{2}}{\varepsilon _{j}^{2}}}$

$\leq {\frac {1}{\beta ^{k+1}}}\prod _{j=0}^{k}\left({\frac {\left\|B_{j}^{-1}s^{j}\right\|^{2}}{\varepsilon _{j}}}{\frac {\sigma _{j}}{\varepsilon _{j}}}\right)\leq {\frac {1}{\det(B_{0}^{-1})}}{\frac {1}{\beta ^{k+1}}}c_{1}^{k+1}\left\{{\frac {2}{n}}c_{1}(k+1)\right\}^{n}\leq c_{2}^{k+1},$

wobei $c_{2}>0$ eine gewisse Konstante ist. Damit ist auch die Gültigkeit von (7.30) bewiesen.

q.e.d.

Damit haben wir die globale Konvergenz des BFGS-Verfahrens für gleichmäßig konvexe Funktionen bewiesen. Obwohl das BFGS-Verfahren in der Praxis auch für andere nichtlineare Funktionen sehr robust ist, ist bisher nicht bekannt, ob es für beliebige Funktionen $f\in C^{1}(\mathbb {R} ^{n})$ und zumindest für einige Schrittweitenregeln in dem Sinne global konvergent ist, dass jeder Häufungspunkt einer durch das Verfahren erzeugten Folge ein stationärer Punkt von $f$ ist. Letzteres kann man z. B. für das globalisierte Newton-Verfahren mit der Armijo-Schrittweitenregel zeigen (s. [GeiKa99]). Es sei aber bemerkt, dass jede Funktion in der Umgebung eines lokalen Minimalpunktes, in dem die hinreichenden Optimalitätsbedingungen zweiter Ordnung aus Satz 1.14 erfüllt sind, gleichmäßig konvex ist. Denn es gilt (Beweis: Übung!):

Lemma 7.17

Es sei $f\in C^{2}(\mathbb {R} ^{n})$. Man zeige, dass die Bedingung

(7.42) ${\frac {1}{M}}\|u\|^{2}\leq u^{T}\nabla ^{2}f(x)u\leq {\frac {1}{m}}\|u\|^{2},\quad u\in \mathbb {R} ^{n},\quad x\in {\mathcal {U}}_{\varepsilon }(x^{*})$

genau dann für ein $x^{*}$ und ein $\varepsilon >0$ erfüllt ist, wenn $\nabla ^{2}f(x^{*})$ positiv definit ist.

Der Nachweis der superlinearen Konvergenz des BFGS-Verfahrens ist noch aufwändiger als der seiner Konvergenz. Es werden dazu auch stärkere Glattheitsvoraussetzungen benötigt.

Satz 7.18

Es seien (V1) - (V4) erfüllt und für die dann existierende eindeutige Lösung $x^{*}$ von $(P)$ gebe es Konstanten $L>0$ und $\varepsilon >0$, für die $f\in C^{2}({\mathcal {U}}_{\varepsilon }(x^{*}))$ gilt sowie

(7.43) $\left\|\nabla ^{2}f(x)-\nabla ^{2}f(x^{*})\right\|\leq L\|x-x^{*}\|,\quad x\in {\mathcal {U}}_{\varepsilon }(x^{*})$

Weiter gelte

(7.44) $\lim _{k\to \infty }{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}=0$

und

(7.45) $\lim _{k\to \infty }t_{k}=1.$

Dann konvergiert die durch Algorithmus 7.13 erzeugte Folge $\left\{x^{k}\right\}$ superlinear gegen die eindeutige Lösung $x^{*}$ von $(P)$.

Beweis.

Wir wollen zeigen, dass

(7.46) $\lim _{k\to \infty }{\frac {\left\|g^{k+1}\right\|}{\left\|x^{k+1}-x^{k}\right\|}}=0$

gilt. Nach Lemma 2.9 (iv) hat man nämlich

$\left\|g^{k+1}\right\|\geq \beta \left\|x^{k+1}-x^{*}\right\|,$

so dass wegen

$\left\|x^{k+1}-x^{k}\right\|\leq \left\|x^{k+1}-x^{*}\right\|+\left\|x^{k}-x^{*}\right\|$

aus (7.46) folgt:

$\lim _{k\to \infty }{\frac {\left\|x^{k+1}-x^{*}\right\|}{\left\|x^{k+1}-x^{*}\right\|+\left\|x^{k}-x^{*}\right\|}}=\lim _{k\to \infty }{\frac {\left\|{\frac {x^{k+1}-x^{*}}{x^{k}-x^{*}}}\right\|}{\left\|{\frac {x^{k+1}-x^{*}}{x^{k}-x^{*}}}\right\|+1}}=0.$

Daraus erschließt man die superlineare Konvergenz von $\left\{x^{k}\right\}$, da

$\lim _{k\to \infty }{\frac {z_{k}}{z_{k}+1}}=0$

genau dann für eine Folge von Zahlen $z_{k}\geq 0$ gilt, wenn $z_{k}\to 0$ $(k\to \infty )$ ist.

Wegen

$g^{k}=-B_{k}^{-1}p^{k}=-{\frac {1}{t_{k}}}B_{k}^{-1}(x^{k+1}-x^{k})$

erhalten wir nun im Hinblick auf den Zähler in (7.46)

(7.47) $\left\|g^{k+1}\right|=\left\|g^{k+1}-g^{k}-\nabla ^{2}f(x^{*})(x^{k+1}-x^{k})-\left({\frac {1}{t_{k}}}B_{k}^{-1}-\nabla ^{2}f(x^{*})\right)(x^{k+1}-x^{k})\right\|.$

Für alle hinreichend großen $k$ gilt weiter gemäß Satz 7.16 $x^{k},x^{k+1}\in {\mathcal {U}}_{\varepsilon }(x^{*})$ und liefern daher eine Taylor-Entwicklung von $g^{k+1}$ mit einem $\vartheta _{k}\in (0,1)$ und die Anwendung von (7.43)

$\left\|g^{k+1}-g^{k}-\nabla ^{2}f(x^{*})(x^{k+1}-x^{k})\right\|=\left\|\left[\nabla ^{2}f(x^{k}+\vartheta _{k}(x^{k+1}-x^{k}))-\nabla ^{2}f(x^{*})\right](x^{k+1}-x^{k})\right\|$

(7.48) Fehler beim Parsen (Syntaxfehler): {\displaystyle \le L \left\| x^k + \vartheta_k(x^{k+1} - x^k) - x^* \right\| \left\| x^{k+1} - x^k \right\| \right\| \le L \left\{ \left\|x^k - x^*\right\| + \left\| x^{k+1} - x^k \right\| \right\} \left\| x^{k+1} - x^k \right\|.}

Außerdem hat man

$\left\|{\frac {\left({\frac {1}{t_{k}}}B_{k}^{-1}-\nabla ^{2}f(x^{*})\right)\left(x^{k+1}-x^{k}\right)}{x^{k+1}-x^{k}}}\right\|=\left\|{\frac {\left({\frac {1}{t_{k}}}B_{k}^{-1}-B_{k}^{-1}+B_{k}^{-1}-\nabla ^{2}f(x^{*})\right)p^{k}}{p^{k}}}\right\|$

(7.49) $\leq \left|{\frac {1}{t_{k}}}-1\right|{\frac {\left\|B_{k}^{-1}p^{k}\right\|}{\left\|p^{k}\right\|}}+{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}.$

Der zweite Summand in (7.49) konvergiert gemäß der Voraussetzung (7.44) gegen 0 für $k\to \infty$ und der erste tut dies aufgrund der Voraussetzung (7.45) und aufgrund der Beschränktheit der Folge $\left\{\left\|B_{k}^{-1}p^{k}\right\|/\left\|p^{k}\right\|\right\}$, welche sich wegen

${\frac {\left\|\nabla ^{2}f(x^{*})p^{k}\right\|}{\left\|p^{k}\right\|}}\leq \left\|\nabla ^{2}f(x^{*})\right\|$

aus (7.44) ergibt. Den Grenzwert in (7.46) gewinnt man nun mit den Abschätzungen (7.47), (7.48) und (7.49) und wegen der aus Satz 7.16 hervorgehenden Konvergenz $x^{k}\to x^{*}$ $(k\to \infty )$.

q.e.d.

Der folgende Satz zeigt, dass es für den Nachweis der superlinearen Konvergenz des BFGS-Verfahrens genügt, die Bedingung (7.44) nachzuweisen, wenn das Verfahren mit einer der Schrittweitenregeln aus Kapitel 3 verknüpft wird. Wir verwenden in diesem Zusammenhang wieder die Voraussetzung (V5), die gemäß Bemerkung 6.8 das Erfülltsein der Bedingungen (V1) - (V4) nach sich zieht (siehe auch Lemma 7.17):

(V5) Es ist $f\in C^{2}(\mathbb {R} ^{n})$, die Menge $N_{0}$ aus (2.9) ist konvex und es existieren Konstanten $0<m\leq M$ mit

(7.50) ${\frac {1}{M}}\|u\|^{2}\leq u^{T}\nabla ^{2}f(x)u\leq {\frac {1}{m}}\|u\|^{2},\quad u\in \mathbb {R} ^{n},\quad x\in N_{0}.$

Satz 7.19

Es sei (V5) erfüllt und es sei

(7.51) $\lim _{k\to \infty }{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}=0.$

Dann gilt bei Verwendung der Minimum-, Curry-, Armijo-, Wolfe-Powell- oder strengen Wolfe-Powell-Schrittweitenregel im Algorithmus 7.13:

(i) $\lim _{k\to \infty }t_{M}(x^{k},p^{k})=\lim _{k\to \infty }t_{C}(x^{k},p^{k})=1.$

(ii) $t_{A}(x^{k},p^{k})=1$ für alle hinreichend großen $k\in \mathbb {N}$.

(iii) Für alle hinreichend großen $k$ ist die Wahl

$t_{WP}(x^{k},p^{k})=t_{SWP}(x^{k},p^{k})=1$

möglich.

Beweis.

Sei für jedes $k$ zunächst $t_{k}:=t_{M}(x^{k},p^{k})$ oder $t_{k}:=t_{C}(x^{k},p^{k})$. Bekanntlich gilt dann die Beziehung

$\nabla f(x^{k}+t_{k}p^{k})^{T}p^{k}=\left(g^{k+1}\right)^{T}p^{k}=0.$

Für alle $k$ liegen nun offenbar $x^{k}$ und $x^{k+1}$ und damit auch alle Punkte auf der Verbindungsstrecke zwischen diesen Vektoren in der gemäß (V5) konvexen Menge $N_{0}$. Für ein ${\tilde {x}}^{k}\in N_{0}$ zwischen $x^{k}$ und $x^{k+1}$ gewinnt man daher mittels einer Taylor-Entwicklung

$-(g^{k})^{T}p^{k}=(g^{k+1}-g^{k})^{T}p^{k}=t_{k}(p^{k})^{T}\nabla ^{2}f({\tilde {x}}^{k})p^{k}.$

Also bekommt man

$t_{k}={\frac {-(g^{k})^{T}p^{k}}{(p^{k})^{T}\nabla ^{2}f({\tilde {x}}^{k})p^{k}}}=1-{\frac {\left[\nabla ^{2}f({\tilde {x}}^{k})p^{k}+g^{k}\right]^{T}p^{k}}{(p^{k})^{T}\nabla ^{2}f({\tilde {x}}^{k})p^{k}}}$

und damit unter Verwendung von (V5)

$|t_{k}-1|\leq {\frac {M\left\|\nabla ^{2}f({\tilde {x}}^{k})p^{k}+g^{k}\right\|\left\|p^{k}\right\|}{\left\|p^{k}\right\|^{2}}}={\frac {M\left\|\nabla ^{2}f({\tilde {x}}^{k})p^{k}+g^{k}\right\|}{\left\|p^{k}\right\|}}$

$\leq M\left\{{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}+\left\|\nabla ^{2}f({\tilde {x}}^{k})-\nabla ^{2}f(x^{*})\right\|\right\}.$

Aus (7.51) und aus der mit Satz 7.16 gesicherten Konvergenz $x^{k}\to x^{*}$ $(k\to \infty )$ erschließt man somit den Grenzwert in Aussage (i).

Sei nun $t_{k}:=t_{A}(x^{k},p^{k})$ für alle $k$. Im Hinblick auf Aussage (ii) beachte man zunächst, dass wegen $-g^{k}=B_{k}^{-1}p^{k}$, Voraussetzung (V5) und (7.51) für alle $k\geq k_{0}$ mit einem $k_{0}\in \mathbb {N}$ gilt:

${\frac {-(g^{k})^{T}p^{k}}{\|p^{k}\|^{2}}}={\frac {\left[-g^{k}-\nabla ^{2}f(x^{*})p^{k}\right]^{T}p^{k}}{\left\|p^{k}\right\|^{2}}}+{\frac {(p^{k})^{T}\nabla ^{2}f(x^{*})p^{k}}{\left\|p^{k}\right\|^{2}}}$

(7.52) $\geq -{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}+{\frac {\left\|p^{k}\right\|^{2}}{M\left\|p^{k}\right\|^{2}}}\geq {\frac {1}{2M}}.$

Demzufolge hat man

(7.53) ${\frac {1}{2M}}\left\|p^{k}\right\|^{2}\leq (g^{k})^{T}p^{k}\leq \left\|g^{k}\right\|\left\|p^{k}\right\|,\quad k\geq k_{0}$

bzw.

${\frac {1}{2M}}\left\|p^{k}\right\|\leq \left\|g^{k}\right\|,\quad k\geq k_{0}$

Da die Konvergenz $x^{k}\to x^{*}$ $(k\to \infty )$ (vgl. Satz 7.16) impliziert, dass

$\lim _{k\to \infty }g^{k}=\nabla g(x^{*})=0$

gilt, folgt daraus $\lim _{k\to \infty }p^{k}=0$.

Mittels einer Taylor-Entwicklung schließen wir weiter für ein $\eta _{k}\in (0,1)$ und für alle $k\geq k_{0}$

${\frac {f(x^{k})-f(x^{k}+p^{k})}{-(g^{k})^{T}p^{k}}}={\frac {-\nabla f(x^{k})^{T}p^{k}-{\frac {1}{2}}(p^{k})^{T}\nabla ^{2}f(x^{k}+\eta _{k}p^{k})p^{k}}{-(g^{k})^{T}p^{k}}}=1-{\frac {1}{2}}{\frac {(p^{k})^{T}\nabla ^{2}f(x^{k}+\eta _{k}p^{k})p^{k}}{-(g^{k})^{T}p^{k}}}$

(7.54) $={\frac {1}{2}}-{\frac {1}{2}}\delta _{k}$

mit

$\delta _{k}:={\frac {\left[\nabla ^{2}f(x^{k}+\eta _{k}p^{k})p^{k}+g^{k}\right]^{T}p^{k}]}{-(g^{k})^{T}p^{k}}}.$

Dieses $\delta _{k}$ können wir mit (7.53) nach oben abschätzen durch

$|\delta _{k}|\leq {\frac {\left|\left[\nabla ^{2}f(x^{k}+\eta _{k}p^{k})p^{k}+g^{k}\right]^{T}p^{k}\right|]}{{\frac {1}{2M}}\left\|p^{k}\right\|^{2}}}$

$={\frac {2M}{\left\|p^{k}\right\|}}\left\|\nabla ^{2}f(x^{k}+\eta _{k}p^{k})p^{k}-\nabla ^{2}f(x^{*})p^{k}+\nabla ^{2}f(x^{*})p^{k}-B_{k}^{-1}p^{k}\right\|$

$\leq 2M\left\{\left\|\nabla ^{2}f(x^{k}+\eta _{k}p^{k})-\nabla ^{2}f(x^{*})\right\|+{\frac {\left\|\left(B_{k}^{-1}-\nabla ^{2}f(x^{*})\right)p^{k}\right\|}{\left\|p^{k}\right\|}}\right\}.$

Da $\left\{x^{k}\right\}$ gegen $x^{*}$ und $\left\{p^{k}\right\}$ gegen 0 konvergieren, strebt der erste Ausdruck auf der rechten Seite der letzten Ungleichung gegen 0, so dass in Verbindung mit (7.51) $\lim _{k\to \infty }\delta _{k}=0$ folgt. Aus (7.54) ergibt sich somit für jedes $\zeta \in (0,1/2)$ und alle hinreichend großen $k$ die Abschätzung

${\frac {f(x^{k})-f(x^{k}+p^{k})}{-(g^{k})^{T}p^{k}}}>\zeta .$

Gemäß der Definition 3.10 der Armijo-Schrittweitenregel ist für diese $k$ daher $t_{k}=1$.

Die erste Ungleichung der Wolfe-Powell- bzw. strengen Wolfe-Powell-Schrittweitenregel ist gemäß Aussage (ii) für alle hinreichend großen $k$ für $t_{k}:=1$ erfüllt. Weiter sind auch hier die Abschätzungen (7.52) und (7.53) gültig, so dass man ebenfalls $\lim _{k\to \infty }p^{k}=0$ hat. Unter Verwendung von (7.53) und der Identität $p^{k}=-B_{k}g^{k}$ schließt man nun mit einem $\theta _{k}\in (0,1)$

$\left|{\frac {-\nabla f(x^{k}+p^{k})^{T}p^{k}}{-\nabla f(x^{k})^{T}p^{k}}}\right|=\left|{\frac {-(g^{k})^{T}p^{k}-(p^{k})^{T}\nabla ^{2}f(x^{k}+\theta _{k}p^{k})p^{k}}{-(g^{k})^{T}p^{k}}}\right|$

$\leq {\frac {2M}{\left\|p^{k}\right\|}}\left\|-\nabla ^{2}f(x^{k}+\theta _{k}p^{k})p^{k}+\nabla ^{2}f(x^{*})p^{k}-\nabla ^{2}f(x^{*})p^{k}+B_{k}^{-1}p^{k}\right\|$

$\leq 2M\left\{\left\|-\nabla ^{2}f(x^{k}+\theta _{k}p^{k})+\nabla ^{2}f(x^{*})\right\|+{\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}\right\}.$

Grenzwertbildung für $k\to \infty$ liefert wegen der Konvergenz $x^{k}\to x^{*}$ $(k\to \infty )$ und wegen (7.51) das gewünschte Ergebnis (siehe den letzten Teil des Beweises von Satz 6.10).

q.e.d.

Der Beweis dafür, dass die Bedingung (7.51) für das BFGS-Verfahren erfüllt ist, ist der schwierigste und aufwändigste Teil des Nachweises der superlinearen Konvergenz des BFGS-Verfahrens (siehe [Wer92] für einen solchen Beweis). Die Bedingung (7.51) geht auf Dennis und Moré zurück und ist zentral für die superlineare Konvergenz von Quasi-Newton-Verfahren. Zusammen mit der Forderung an die Schrittweiten, dass $t_{k}:=1$ für alle hinreichend großen $k$ gewählt wird, ist sie nicht nur hinreichend, sondern auch notwendig für die superlineare Konvergenz von Quasi-Newton-Verfahren (z. B. [GeiKa99, S. 55ff.], [SuYu06, S. 241ff.]).

Da

${\frac {\left\|\left[B_{k}^{-1}-\nabla ^{2}f(x^{*})\right]p^{k}\right\|}{\left\|p^{k}\right\|}}\leq \left\|B_{k}^{-1}-\nabla ^{2}f(x^{*})\right\|$

gilt, bedeutet die Bedingung (7.51) nicht, dass die Matrizen $B_{k}^{-1}$ für $k\to \infty$ gegen $\nabla ^{2}f(x^{*})$ konvergieren müssen, sondern nur, dass sie die Matrix $\nabla ^{2}f(x^{*})$ entlang der Richtungen $p^{k}$ immer genauer anzunähern haben. Allerdings kann man für einige Funktionenklassen tatsächlich auch zeigen, dass $B_{k}^{-1}\to \nabla ^{2}f(x^{*})$ für $k\to \infty$ gilt (siehe die Hinweise in [GeiKa99]).

7.5 Bemerkungen zur Numerik

Als Startmatrix $B_{0}$ für das BFGS-Verfahren wählt man häufig ein skalares Vielfaches der Einheitsmatrix, wobei der Skalar dazu dient, die Variablen geeignet zu skalieren, oder man wählt eine Näherung der Hesse-Matrix von $f$ in $x^{0}$, indem man z. B. die partiellen Ableitungen mittels sog. finiter Differenzen approximiert. Als Schrittweitenregel für das BFGS-Verfahren wird in der Praxis häufig die Wolfe-Powell- oder die strenge Wolfe-Powell-Regel verwendet. Denn z. B. bei Verwendung der Armijo-Schrittweitenregel kann für Funktionen, die nicht gleichmäßig konvex sind, nicht garantiert werden, dass $\sigma _{k}>0$ ist (vgl. Lemma 7.8). Für Details und weitere Hinweise zur Implementation des Verfahrens verweisen wir z. B. auf [NoWri06, S. 142ff.].

Das BFGS-Verfahren zeigt in der Praxis meist sehr viel schnellere Konvergenz als das Gradientenverfahren. Für nichtquadratische Funktionen ist es auch schneller als die CG-Verfahren, wobei es aber im Vergleich mit diesen erheblich mehr Speicherplatz benötigt. Letzteres ist ein Nachteil bei sehr großen Problemen.

Überdies ist das BFGS-Verfahren häufig auch numerisch effizienter als das Newton-Verfahren, obwohl dieses unter geeigneten Voraussetzungen eine höhere Konvergenzrate aufweist. Dies liegt daran, dass für die Berechnung der Matrix $B_{k}$ und der Richtung $p^{k}$ nur jeweils (und damit auch zusammen nur) ${\mathcal {O}}(n^{2})$ arithmetische Rechenoperationen erforderlich sind, während beim Newton-Verfahren in jeder Iteration ein lineares Gleichungssystem gelöst werden muss, wofür ${\mathcal {O}}(n^{3})$ Rechenoperationen benötigt werden. Überdies muss beim Newton-Verfahren in jeder Iteration die Hesse-Matrix von $f$ in der aktuellen Näherung exakt oder auf numerischem Wege ermittelt werden. Für Hinweise zum numerischen Verhalten anderer Quasi-Newton-Verfahren wie dem des Kleinmichel- und des DFP-Verfahrens verweisen wir auf die angegebene Literatur (z.B. [GeiKa99], [NoWri06]).

Ein Nachteil von Quasi-Newton-Verfahren insbesondere im Hinblick auf die Lösung großer Optimierungsprobleme ist es, dass die Update-Matrizen häufig nicht die Struktur der Hesse-Matrizen widerspiegeln, also (fast) vollbesetzt sein können, auch wenn die Hesse-Matrizen selbst nur dünn besetzt sind. Daher sollte man Quasi-Newton-Verfahren nur für Optimierungsprobleme kleiner und mittlerer Größe mit maximal wenigen Hundert Variablen verwenden.

Im Hinblick auf die Lösung großer und sehr großer unrestringierter Optimierungsprobleme sind Quasi-Newton-Methoden auf verschiedene Weisen modifiziert worden (z. B. [SuYu06], [NoWri06]). Eine bekannte Klasse von Verfahren ist die der Limited-Memory-Quasi-Newton-Verfahren. Wir verweisen für solche Modifikationen von Quasi-Newton-Verfahren wieder auf die Literatur (z. B. [GeiKa99], [NoWri06]).