3.3. Verfahren der konjugierten Gradienten

Im Folgenden wollen wir uns mit einem besonders eleganten Verfahren der Optimierung beschäftigen: dem Verfahren der konjugierten Gradienten. Ursprünglich wurde das Verfahren von Hestenes und Stiefel in [HS52] im Jahr \(1952\) vorgeschlagen. Obwohl das Verfahren im Allgemeinen für die nichtlineare Optimierung eingesetzt werden kann, wird es insbesondere zur Lösung von großen linearen Gleichungssystemen \(Ax = b\) mit symmetrischer, dünn besetzter, positiv definiter Matrix \(A \in \R^{n\times n}\) eingesetzt. Solche Gleichungssysteme treten zum Beispiel bei der numerischen Modellierung und Lösung partieller Differentialgleichungen auf. Das Verfahren lässt sich in diesem Fall besonders anschaulich motivieren und herleiten. Darum wollen wir uns im Folgenden zunächst auf das Lösen von großen linearen Gleichungssystemen \(Ax=b\) konzentrieren. Wir folgen bei der Herleitung des Verfahren der konjugierten Gradienten der didaktisch sehr gelungenen Arbeit von Jonathan Shewchuk in [She94]. Für eine ansprechende, interaktive Visualisierung des Verfahren der konjugierten Gradienten empfehlen wir den Mathematik Blog von Philipp Wacker [Wac].

3.3.1. Problemstellung

Sei im Folgenden also \(A \in \mathbb{R}^{n \times n}\) eine sehr große Matrix und \(b \in \mathbb{R}^n\) ein reeller Vektor. Wir suchen einen unbekannten Vektor \(x \in \mathbb{R}^n\), der das lineare Gleichungssystem

(3.22)\[Ax \ = \ b\]

löst. Wir suchen also nach denjenigen Koeffizienten, mit denen sich der Vektor \(b\) als Linearkombination aus Spaltenvektoren der Matrix \(A\) darstellen lässt. Diese Koeffizienten entsprechen den Einträgen des unbekannten Vektors \(x\).

Aus der Vorlesung Einführung in die Numerik  in [TR] ist bekannt, dass es genau dann eine eindeutige Lösung \(x \in \mathbb{R}^n\) für die Gleichung (3.22) gibt, falls die Determinante \(\operatorname{det}(A) \neq 0\) ist. Eine hinreichende Bedingung für die Eindeutigkeit des Lösungsvektors \(x \in \mathbb{R}^n\) ist es also zu fordern, dass die Matrix \(A\) symmetrisch und positiv definit ist.

Wir gehen aus diesem Grund im Folgenden immer davon aus, dass \(A \in \mathbb{R}^{n\times n}\) eine symmetrische und positiv definite Matrix ist. In diesem Fall ist das Bestimmen einer Lösung von (3.22) ein gut-gestelltes Problem und die Lösung lässt sich direkt angeben als:

\[x \ = \ A^{-1}b.\]

Wie wir jedoch ebenfalls aus [TR] wissen ist die Inversion einer Matrix \(A \in \mathbb{R}^{n\times n}\) numerisch sehr aufwendig und selbst unter Ausnutzung der Symmetrie lässt sich höchstens ein Verfahren mit Rechenaufwand \(\mathcal{O}(\frac{1}{6}n^3)\) angeben. Sollte die Dimension des Problems jedoch sehr groß sein (d.h. wir nehmen \(n \gg 1\) an), so ist eine direkte Lösung von (3.22) mittels Inversion nicht durchführbar. Glücklicherweise liefert uns das Verfahren der konjugierten Gradienten (neben anderen iterativen Lösungsalgorithmen) eine Möglichkeit, das lineare Gleichungssystem numerisch zu lösen. Man beachte, dass wir explizit darauf verzichten zu fordern, dass die Matrix \(A\) dünnbesetzt ist. In diesem Fall könnten wir nämlich ebenfalls die numerischen Iterationsverfahren aus [TR] anwenden.

Wir betrachten zunächst das folgende konvexe, quadratische Optimierungsproblem der Form

(3.23)\[\min_{x\in \mathbb{R}^n} \, \frac{1}{2} \langle x, Ax \rangle - \langle b, x \rangle + c,\]

wobei \(A\) und \(b\) wie im Fall des linearen Gleichungssystems in (3.22) gewählt sind und \(c \in \mathbb{R}\) eine beliebige, reelle Konstante ist. Der folgende Satz liefert uns eine hilfreiche Aussage zur Lösung des ursprünglichen Problems.

Satz 3.8 (Äquivalenzaussage für lineares Gleichungssystem)

Das konvexe, quadratische Minimierungsproblem in (3.23) ist äquivalent zum ursprünglichen linearen Gleichungssystem in (3.22), d.h., jede Lösung von (3.23) ist schon Lösung von (3.22) und anders herum.

Proof. Für die erste Richtung des Beweises nehmen wir an, dass \(x^* \in \mathbb{R}^n\) eine Lösung des linearen Gleichungssystems \(Ax = b\) sei. Wir betrachten die hinreichenden Optimalitätsbedingungen zweiter Ordnung aus Satz 3.3 für das Minimierungsproblem (3.23). Hierzu suchen wir zunächst die stationären Punkte der Funktion \(f \colon \mathbb{R}^n \rightarrow \mathbb{R}\) mit

\[f(x) \ \coloneqq \ \frac{1}{2} \langle x, Ax \rangle - \langle b, x \rangle + c.\]

Der Gradient von \(f\) lässt sich wegen der Symmetrie von \(A\) bestimmen als

\[\nabla f(x) \ = \ \frac{1}{2}(A + A^\top)x - b \ = \ Ax - b \ \overset{!}{=} \ 0.\]

Alle stationären Punkte \(x \in \mathbb{R}^n\) von \(f\) mit \(\nabla f(x) = 0\) sind also gerade die Lösungen des linearen Gleichungssystems \(Ax = b\). Damit ist \(x^*\) nach Voraussetzung also einziger stationärer Punkt von \(f\). Um zu zeigen, dass \(x^* \in \mathbb{R}^n\) auch schon ein lokales Minimum von \(f\) ist müssen wir noch die Hessematrix von \(f\) betrachten, welche gegeben ist durch:

\[\nabla^2 f(x) \ = \ A.\]

Da \(A\) nach Voraussetzung positiv definit ist, ist auch die Hessematrix \(\nabla^2 f(x) = A\) positiv definit und somit sind die hinreichenden Kriterien für das Vorliegen eines lokalen Minimums von \(f\) im Punkt \(x^* \in \mathbb{R}^n\) erfüllt.

Für die Rückrichtung des Beweises nehmen wir, dass \(x^* \in \mathbb{R}^n\) ein lokales Minimum der Zielfunktion \(f\) ist. Damit folgt direkt, dass \(x^*\) ein stationärer Punkt von \(f\) ist und somit muss gelten:

\[\nabla f(x^*) \ = \ A x^* - b \ = \ 0.\]

Das bedeutet aber schon, dass \(x^* \in \mathbb{R}^n\) Lösung des linearen Gleichungssystems \(Ax = b\) ist. ◻

Satz 3.8 erlaubt es uns also, ein quadratisches Optimierungsproblem der Form (3.23) numerisch zu lösen, anstatt einen unbekannten Lösungsvektor für ursprüngliche lineare Gleichungssystem (3.22) zu finden.

Wir interessieren uns nun also für ein iteratives Verfahren, welches eine Folge von Punkten \(x_0,x_1,\ldots \in \mathbb{R}^n\) konstruiert, die gegen ein Minimum von (3.23) und somit gegen die eindeutige Lösung des linearen Gleichungssystems (3.22) konvergiert. Hierfür benötigen wir noch zusätzliche Notation, um das angestrebte Verfahren vernünftig zu beschreiben.

Definition 3.9 (Fehler und Residuum)

Sei \(x_{k+1} = G(x_k)\) ein Iterationsverfahren, dass gegen ein lokales Minimum \(x^* \in \mathbb{R}^n\) der quadratischen Funktion \(f\) in (3.23) konvergiert, d.h., \(x_k \rightarrow x^*\) für \(k \rightarrow \infty\). Dann können wir die beiden folgenden Begriffe definieren:

1. Wir bezeichnen den Vektor \(e_k \in \mathbb{R}^n\) mit

   \[e_k \ \coloneqq \ x_k - x^*\]

   als den aktuellen Fehler, den man durch den aktuellen Punkt \(x_k \in \mathbb{R}^n\) macht.

2. Wir bezeichnen den Vektor \(r_k \in \mathbb{R}^n\) mit

   \[r_k \ \coloneqq \ b - Ax_k\]

   als das aktuelle Residuum, das man durch den aktuellen Punkt \(x_k \in \mathbb{R}^n\) erhält.

Bemerkung 3.9 (Fehler und Residuum)

In Bezug auf Definition 3.9 lassen sich folgende Aussagen festhalten:

1. Der Fehler \(e_k \in \mathbb{R}^n\) ist eher abstrakter Natur und dient zur besseren Analyse des Verfahrens der konjugierten Gradienten. Explizit werden wir diesen Vektor innerhalb des Iterationsverfahrens jedoch nie bestimmen können, da wir dann schon fertig wären mit einem einfachen Update der Form \(x^* = x_k - e_k\).

2. Wie wir bereits im Beweis von Satz 3.8 gesehen haben, lässt sich das Residuum \(r_k \in \mathbb{R}^n\) außerdem wie folgt umschreiben:

   \[r_k \ = \ \underbrace{b - Ax_k}_{= -\nabla f(x_k)} \ = \ Ax^* - Ax_k \ = \ A (x^* - x_k) \ = \ -Ae_k.\]

   Daher lässt sich das Residuum \(r_k\) auch als die Richtung des stärksten Abstiegs interpretieren und es ist klar, dass \(r_k\) immer orthogonal zu den Niveaulinien der Funktion \(f\) steht.

Abb. 3.2 illustriert anschaulich die geometrische Bedeutung der beiden in Definition 3.9 eingeführten Vektoren. Wie man unschwer erkennt zeigen Fehler und Residuum im Allgemeinen nicht in die selbe Richtung. Das erklärt auch warum das Gradientenabstiegsverfahren in Gradientenabstiegsverfahren selbst bei optimaler Schrittweite \(\tau_k > 0\) nicht in einem Schritt die gesuchte Lösung \(x^* \in \mathbb{R}^n\) erreicht.

Abb. 3.2 Visualisierung des Fehlers \(e_0 \in \mathbb{R}^n\) und des Residuums \(r_0 \in \mathbb{R}^n\) für einen Startpunkt \(x_0 \in \mathbb{R}^n\).

3.3.2. Motivation

Um das Vorgehen beim Verfahren der konjugierten Gradienten zu motivieren rufen wir uns noch einmal das Gradientenabstiegsverfahren aus Gradientenabstiegsverfahren in Erinnerung. Nehmen wir an wir befinden uns im \(k\)-Schritt des Gradientenabstiegsverfahrens in Algorithmus 3.1 in einem Punkt \(x_k \in \mathbb{R}^n\) und es sei eine Schrittweite \(\tau_k > 0\) gegeben. Dann erhalten wir den nächsten Punkt \(x_{k+1} \in \mathbb{R}^n\) der Iterationsfolge durch folgendes Update:

\[x_{k+1} \ = \ x_k - \tau_k \nabla f(x_k) \ = \ x_k + \tau_k r_k,\]

wobei \(r_k \in \mathbb{R}^n\) das aktuelle Residuum des Punktes \(x_k\) bezeichnet. Wir machen also in Richtung des steilsten Gradientenabstiegs einen Schritt skaliert mit der Schrittweite \(\tau_k > 0\). Da die Abstiegsrichtung in jedem Schritt \(x_k \rightarrow x_{k+1}\) orthogonal zu den Niveaulinien von \(f\) steht, erhält man durch das Gradientenabstiegsverfahren typischerweise einen Zickzack-Pfad.

Um dieses typische Verhalten besser zu verstehen können wir eine Vorüberlegung zur Schrittweitenwahl für das quadratische Optimierungsproblem in (3.23) machen. Hierzu halten wir die Schrittrichtung \(p_k \coloneqq -\nabla f(x_k)\) fest und wollen bezüglich der unbekannten Schrittweite optimieren.

Wir gehen davon aus, dass wir das lokale Minimum von \(f\) noch nicht erreicht haben, denn dann wäre \(\tau_k = 0\). Wir suchen also eine Schrittweite \(\tau > 0\), so dass der Funktionswert \(f(x_{k+1})\) entlang der Linie \(x_k - \tau \nabla f(x_k)\) minimal wird. Da \(f\) eine quadratische Funktion ist, wissen wir, dass ein eindeutiges Minimum \(\tau_k\) entlang dieser Linie existieren muss. Wir nutzen also die notwendigen Optimalitätsbedingungen aus Satz 3.1 für das totale Differential, um folgenden Zusammenhang herzustellen:

\[\begin{split}\begin{split} \frac{\mathrm{d}}{\mathrm{d}\tau} f(x_{k+1}) \ &= \ \Bigl\langle \nabla f(x_{k+1}), \frac{\mathrm{d}x_{k+1}}{\mathrm{d}\tau} \Bigr\rangle \ = \ \Bigl\langle \nabla f(x_{k+1}), \frac{\mathrm{d} (x_k - \tau \nabla f(x_k))}{\mathrm{d}\tau} \Bigr\rangle \\ \ &= \ \langle \nabla f(x_{k+1}), -\nabla f(x_k) \rangle \ = \ \langle \nabla f(x_{k+1}), p_k \rangle \ \overset{!}{=} \ 0. \end{split}\end{split}\]

Das Ergebnis ist durchaus interessant. Die optimale Schrittweite \(\tau > 0\) muss so gewählt werden, dass der nächste Punkt \(x_{k+1} \in \mathbb{R}^n\) der Iterationsfolge an der Stelle liegt an der unsere Abstiegsrichtung orthogonal auf den Gradienten der Funktion \(\nabla f(x_{k+1})\) trifft. Das bedeutet, dass die optimale Abfolge der Abstiegsrichtungen im quadratischen Fall eine Menge von \(90\)-Grad Zickzack-Linien ergibt. Da jedoch der Punkt \(x_{k+1} \in \mathbb{R}^n\) bislang noch unbekannt ist, können wir das optimale \(\tau_k\) nicht in dieser Form angeben. Das folgende Lemma bestimmt die optimale Schrittweite im Fall der quadratischen Optimierung in (3.23).

Lemma 3.2 (Optimale Schrittweite)

Sei \(f \colon \mathbb{R}^n \rightarrow \mathbb{R}\) die quadratische Funktion aus (3.23). Wir betrachten das Gradientenabstiegsverfahren im \(k\)-ten Iterationsschritt mit einer unbekannten Schrittweite \(\tau_k > 0\), die jedoch so gewählt werden muss, dass

\[\langle \nabla f(x_{k+1}), \nabla f(x_k) \rangle \ \overset{!}{=} \ 0.\]

Sei außerdem \(r_k = b - Ax_k\) das Residuum im aktuellen Iterationsschritt. Dann lässt sich die optimale Schrittweite \(\tau_k\) berechnen als:

(3.24)\[\tau_k \ = \ \frac{\langle r_k, r_k \rangle}{\langle r_k, Ar_k \rangle}.\]

Proof. Wir erinnern uns daran, dass \(r_k = -\nabla f(x_k) = b - Ax_{k}\) ist und somit können wir folgern:

\[\begin{split}\begin{split} 0 \ &\overset{!}{=} \ \langle \nabla f(x_{k+1}), \nabla f(x_k) \rangle \ = \ \langle r_{k+1}, r_k \rangle \ = \ \langle b - Ax_{k+1}, r_k \rangle \\ \ &= \ \langle b - A(x_k + \tau_kr_k), r_k \rangle \ = \ \langle b - Ax_k, r_k \rangle - \tau_k \langle Ar_k, r_k \rangle \\ \ &= \langle r_k, r_k \rangle - \tau_k \langle r_k, Ar_k \rangle \end{split}\end{split}\]

Da wir \(A\) als positiv definit angenommen haben, können wir die Gleichung umstellen und erhalten so die behauptete Berechnungsformel für \(\tau_k\) in (3.24). ◻

Obwohl wir die optimale Schrittweite \(\tau_k\) für das quadratische Optimierungsproblem (3.23) bestimmen konnten, ist das Gradientenabstiegsverfahren weit davon entfernt, optimal zu sein. Trotz optimaler Schrittweiten und optimaler Abstiegsrichtungen erhalten wir eine Folge von Richtungsvektoren, die immer wieder in die gleiche Richtung zeigen. Das ist numerisch gesehen äußerst ineffizient. Man könnte sich also fragen, warum man nicht einfach nur zwei orthogonale Schritte macht und die Schrittweiten als Summe der optimalen Schrittweiten der geraden bzw. ungeraden Iterationsschritte \(k \in \mathbb{N}\) wählt. In der Tat würde man für \(N \in \mathbb{N}\) Schritte des Gradientenabstiegsverfahren im selben Punkt \(x_N \in \mathbb{R}^n\) mit nur zwei Iterationen landen, wie in fig:zig-zag illustriert ist.

Vergleich des Gradientenabstiegsverfahrens mit optimaler Schrittweite τ_k > 0 aus (links) mit einem idealen Abstiegsverfahren (rechts), bei dem alle orthogonalen Teilschritte zusammengefasst sind.

Leider können wir nicht alle Schrittweiten aufaddieren, da wir zur Berechnung der optimalen Schrittlänge \(\tau_k > 0\) bereits alle vorangegangen Punkte \(x_k k=0,\ldots,k-1\) kennen müssten. Außerdem würde ein großer, zusammengefasster Schritt in die erste der Richtungen eventuell dazu führen, dass man keinen Abstieg der Funktionswerte von \(f\) entlang der Abstiegsrichtung bis zum Minimum mehr realisiert, sondern im Allgemeinen dieses Minimum überschreitet. Diese Beobachtung ist in Abb. 3.3 illustriert.

Abb. 3.3 Illustration eines idealen Abstiegsverfahrens mit zwei orthogonalen Richtungen. Man beachte, dass die Schrittweite \(\tau_0 > 0\) so gewählt werden muss, dass man im ersten Schritt nicht in einem Punkt \(x_1 \in \mathbb{R}^2\) mit minimalen Funktionswert \(f(x_1)\) entlang der Richtung \(x_0 - \tau_0 \nabla f(x_0)\) endet.

Die Ideallösung wäre natürlich von einem Startpunkt \(x_0 \in \mathbb{R}^n\) in nur einem Schritt zum lokalen Optimum \(x^* \in \mathbb{R}^n\) zu gelangen. Da wir aber den Punkt \(x^*\) a-priori nicht kennen ist das eine unrealistische Forderung. Dennoch lässt sich zeigen, dass das Gradientenabstiegsverfahren mit der optimalen Schrittweite \(\tau_k\) in (3.24) im Fall des quadratischen Optimierungsproblems (3.23) in genau einem Schritt zum lokalen Minimum \(x^* \in \mathbb{R}^n\) führt, wenn der Fehler \(e_0 = x_0 - x^*\) ein Eigenvektor von \(A\) ist. Wir wissen nämlich aus Bemerkung 3.9, dass \(r_k = -Ae_k\) gilt und somit erhalten wir für das Gradientenabstiegsverfahren im ersten Schritt:

\[x_1 \ = \ x_0 - \tau_0 \nabla f(x_0) \ = \ x_0 + \tau_0 r_0 \ = \ x_0 - \tau_0 A e_0 \ = \ x_0 - \tau_0 \lambda e_0.\]

Man müsste bei der Wahl des Startpunktes \(x_0\in \R^n\) jedoch viel Glück haben, um diese Forderung zu erfüllen. Darum wollen wir uns mit alternativen Ideen beschäftigen.

3.3.3. Orthogonale Abstiegsrichtungen

Wir wünschen uns einen Algorithmus, der ähnlich dem Gradientenabstiegsverfahren nur orthogonale Richtungen \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) mit \(d_i \in \mathbb{R}^n\) für \(0 \leq i \leq n-1\) verwendet, jedoch mit der Einschränkung, dass diese nur ein einziges Mal genutzt werden können. Ziel dieses Verfahrens soll es außerdem sein durch \(n\) Schritte in die jeweils \(n\) orthogonalen Richtungen \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) das lokale Minimum der Funktion zu erreichen. Damit hätten wir ein Itrerationsverfahren der Form

(3.25)\[x_{k+1} \ = \ x_k + \tau_k d_k, \quad \tau_k > 0, \ k = 0,\ldots,n-1\]

gewonnen. Wir könnten dies erzwingen indem wir im \(k\)-ten Schritt des Iterationsverfahrens fordern, dass ein Schritt in Richtung \(d_k \in \mathbb{R}^n\) dazu führt, dass der Fehler \(e_{k+1} \in \mathbb{R}^n\) keinerlei Komponenten dieser Richtung mehr enthält, d.h. wir fordern

(3.26)\[\langle e_{k+1}, d_k \rangle \ \overset{!}{=} \ 0.\]

Da wir den zu erwartenden Fehler \(e_{k+1}\) in Bezug auf den aktuellen Punkt \(x_k \in \mathbb{R}^n\) folgendermaßen umschreiben können:

\[e_{k+1} \ = \ x_{k+1} - x^* \ = \ x_k + \tau_k d_k - x^* \ = \ e_k + \tau_k d_k,\]

können wir die Forderung (3.26) umformulieren zu:

\[\langle e_k + \tau_k d_k, d_k \rangle \ \overset{!}{=} \ 0.\]

Hieraus können wir die optimale Schrittweite \(\tau_k > 0\) in Richtung \(d_k \in \mathbb{R}^n\) ableiten als

(3.27)\[\tau_k \ = \ - \frac{\langle e_k, d_k \rangle}{\langle d_k, d_k \rangle}.\]

Obwohl wir in (3.27) eine optimale Schrittweite \(\tau_k\) für das Verfahren mit orthogonalen Abstiegsrichtungen in (3.25) bestimmen konnten, hilft und diese nicht in der praktischen Anwendung des Verfahrens, da sie von dem unbekannten Fehlervektor \(e_k \in \mathbb{R}^n\). Dieser hängt natürlich von der unbekannten Lösung \(x^* \in \mathbb{R}^n\) ab und wenn wir diese kennen würden, so müssten wir kein iteratives Verfahren konstruieren. Selbst wenn man den Fehler \(e_k\) weiter rekursiv umschreibst, so würde man schlussendlich doch bei einer Abhängigkeit des initialen Fehlers \(e_0\) landen. Wir müssen uns also vorerst von dieser Idee verabschieden und nach einer alternativen Möglichkeit suchen.

3.3.4. Konjugierte Abstiegsrichtungen

Obwohl unsere Idee von orthogonalen Abstiegsrichtungen in Orthogonale Abstiegsrichtungen nicht zum Ziel geführt hat, so war die Idee gar nicht schlecht. Das Hauptproblem an der ursprünglichen Idee liegt in der Forderung (3.26), nämlich dass der Fehlervektor \(e_{k+1}\) orthogonal zur aktuellen Richtung \(d_k\) stehen soll. Diese Forderung führt nämlich dazu, dass man orthogonale Vektoren erhält, die nicht an die Geometrie des quadratischen Minimierungsproblems (3.23) angepasst sind.

Wenn man sich die Niveaulinien der Funktion \(f\) genauer anschaut (siehe zum Beispiel Abb. 3.3), so erkennt man, dass es Richtungen gibt entlang derer die Abstiegsrichtung zum lokalen Minimum \(x^* \in \mathbb{R}^n\) steiler verläuft als entlang der anderen Richtungen. Die geometrischen Eigenschaften des Graphen von \(f\) sind maßgeblich durch die Gestalt der Matrix \(A\), genauer gesagt durch deren Eigenvektoren bestimmt. Daher wollen wir diese Eigenschaften bei der Konstruktion eines iterativen Abstiegsverfahren berücksichtigen. Hierzu führen wir folgendes hilfreiche Konzept ein.

Definition 3.10 (Konjugierte Vektoren)

Sei \(A \in \mathbb{R}^{n\times n}\) eine symmetrische, positiv definite Matrix und \(u,v \in \mathbb{R}^n / \lbrace 0\rbrace\) zwei Vektoren. Wir nennen \(v\) und \(w\) konjugiert bezüglich \(\mathbf{A}\) oder auch \(\mathbf{A}\)-orthogonal falls gilt

\[\langle v, Aw \rangle \ = \ \langle w, Av \rangle \ = \ 0.\]

Anstatt nun also die Orthogonalität unserer Richtungsvektoren \(\lbrace d_0, \ldots, d_{n-1}\rbrace\) zu erzwingen wie in Orthogonale Abstiegsrichtungen, fordern wir nun, dass diese Vektoren konjugiert bezüglich der Matrix \(A\) und damit besser an das Problem angepasst sind.

Bemerkung 3.10

Es ist leicht einzusehen, dass \(A\)-orthogonal und orthogonal die selbe Eigenschaft beschreiben, falls die Matrix \(A\) ein Vielfaches der Einheitsmatrix \(I_n \in \mathbb{R}^{n \times n}\) ist. In diesem Fall ist das quadratische Optimierungsproblem (3.23) symmetrisch in alle Richtungen.

Anschaulich lässt sich die Forderung nach \(A\)-Orthogonalität auch so deuten, dass wir ein Paar von Vektoren \(v,w \in \mathbb{R}^n \setminus \lbrace 0\rbrace\) suchen, welche in einem Winkel so zueinander stehen, dass wenn man die Niveaulinien der Funktion \(f\) symmetrisch reskaliert, diese Vektoren anschließend orthogonal zueinander stehen. Diese Idee ist in Abb. 3.4 dargestellt.

Abb. 3.4 Illustration der Geometrie von konjugierten Vektoren im Referenzsystem \(\mathbb{R}^2\) (links) und der selben Vektoren in einem symmetrisierten System bezüglich der Matrix \(A\) (rechts).

Anstatt also ein Abstiegsverfahren der Form (3.25) mit orthogonalen Vektoren zu verwenden, wollen wir ein Abstiegsverfahren mit \(A\)-orthogonalen Vektoren \(\lbrace d_0, \ldots, d_{n-1}\rbrace\) konstruieren, d.h., wir verwenden das Iterationsschema

(3.28)\[x_{k+1} \ = \ x_k + \tau_k d_k, \quad \tau_k > 0, \ k=0,\ldots,n-1\]

wobei für die Abstiegsrichtungen \(d_k \in \mathbb{R}^n\) gelten soll:

\[\langle d_i, Ad_j \rangle = 0 \qquad \text{ für } \ 0 \leq i,j \leq n-1 \ \text{ mit } i\neq j.\]

Wir nehmen für den Moment an, dass wir einen numerischen Algorithmus kennen mit dem wir eine Menge von \(A\)-orthogonalen Vektoren \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) konstruieren können. Wie man diese Menge konkret erhält werden wir uns im Anschluss erschließen. Sei also nun im Folgenden \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) eine gegebene Menge von \(A\)-orthogonalen Vektoren. Dann stellen wir uns die Frage, wie die optimalen Schrittweiten \(\tau_k > 0\) in (3.28) gewählt werden müssen, um in \(n\) Schritten das lokale Minimum \(x^* \in \mathbb{R}^n\) der Funktion \(f\) zu erhalten. Man beachte hierbei, dass wir nicht nur daran interessiert sind den Punkt \(x^*\) genügend gut zu approximieren, sondern wir fordern die eindeutige Lösung des linearen Gleichungssystems \(Ax = b\) in \(n\) Schritten zu finden, d.h., wir nehmen explizit \(x_n = x^*\) an.

Um das lokale Minimum wirklich in \(n\) Schritten zu erreichen müssen wir fordern, dass wir in jede Richtung \(d_k\) nur einmal gehen und der entstehende Fehler \(e_{k+1}\) \(A\)-orthogonal hierzu ist. Das entspricht der Forderung, dass man im entzerrten Problem auf der rechten Seite von Abb. 3.4 nur orthogonale Richtungen verwendet. Wir wollen also folgende Eigenschaft erzwingen:

(3.29)\[\langle Ae_{k+1}, d_k \rangle \ = \ 0.\]

Analog zur Idee der orthogonalen Richtungen in Orthogonale Abstiegsrichtungen können wir den Fehler \(e_{k+1}\) in (3.29) wieder entwickeln, um die optimale Schrittweitenlänge \(\tau_k > 0\) zu bestimmen

\[\begin{split}\begin{split} 0 \ &\overset{!}{=} \ \langle d_k, A e_{k+1} \rangle \ = \ \langle d_k, A(x_{k+1} - x^*) \rangle \ = \ \langle d_k, A(x_k + \tau_k d_k - x^*) \rangle \\ \ &= \ \langle d_k, A(e_k + \tau_k d_k) \rangle \ = \ \langle d_k, -r_k + \tau_k Ad_k \rangle \ = \ \tau_k \langle d_k, Ad_k \rangle - \langle d_k, r_k\rangle. \end{split}\end{split}\]

Da wir \(A\) als positiv definit vorausgesetzt haben, können wir die folgende Gleichung umstellen zu

(3.30)\[\tau_k \ = \ \frac{\langle d_k, r_k \rangle}{\langle d_k, A d_k \rangle}.\]

Im Gegensatz zur Idee der orthogonalen Richtungen in (3.27) lässt sich der Ausdruck in (3.30) explizit berechnen und hängt nicht von dem unbekannten lokalen Minimum \(x^* \in \mathbb{R}^n\) ab. Zusammenfassend heißt das, dass wir aus der Bedingung, dass die Abstiegsrichtung \(d_k \in \mathbb{R}^n\) \(A\)-orthogonal zum Fehlervektor \(e_{k+1} \in \mathbb{R}^n\) sein soll, eine Schrittlänge \(\tau_k > 0\) finden konnten, welche diese Bedingung erfüllt.

Andersherum könnte man fragen, welche Bedingung man aus der Optimalität einer unbekannten Schrittlänge \(\tau > 0\) folgern könnte. Dazu betrachen wir wieder die notwendigen Optimalitätsbedingungen im totalen Differential

\[\begin{split}\begin{split} 0 \ &\overset{!}{=} \ \frac{\mathrm{d}}{\mathrm{d}\tau}f(x_{k+1}) \ = \ \langle \nabla f(x_{k+1}), \frac{\mathrm{d}}{\mathrm{d}\tau}x_{k+1} \rangle \ = \ \langle -r_{k+1}, \frac{\mathrm{d}}{\mathrm{d}\tau} (x_k + \tau d_k) \rangle \\ \ &= \ \langle -r_{k+1}, d_k \rangle \ = \ \langle Ae_{k+1}, d_k \rangle. \end{split}\end{split}\]

Wir erhalten also für die Optimalität der unbekannten Schrittlänge \(\tau > 0\), dass die Abstiegsrichtung \(d_k \in \mathbb{R}^n\) und der Fehlervektor \(e_{k+1}\) konjugiert bezüglich der Matrix \(A\) sein müssen. Das ist aber genau die Eigenschaft, die wir bereits in (3.29) gefordert hatten. Wegen der positiven Definitheit der Matrix \(A\) können wir ebenfalls folgern, dass die Forderung, dass \(e_{k+1}\) und \(d_k\) konjugiert bezüglich \(A\) sind, bei optimaler Schrittweite \(\tau_k\) aus (3.30) in jedem Schritt zu einem Abstieg in Richtung \(d_k\) führt. Wir erhalten also für eine gegebene Menge von \(A\)-orthogonalen Vektoren \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) ein Iterationsverfahren mit optimalen Schrittlangen \(\tau_k > 0\), die wir in (3.30) angeben können und die uns einen Abstieg garantieren.

Zunächst benötigen wir die Einsicht aus folgendem Lemma, die es uns ermöglicht eine Basis aus \(A\)-orthogonalen Vektoren des \(\R^n\) zu betrachten.

Lemma 3.3 (Basis von A-orthogonalen Vektoren)

Sei \(A \in \mathbb{R}^{n\times n}\) eine symmetrische, positiv definite Matrix und sei \(\lbrace d_0, \ldots, d_{n-1}\rbrace\) eine Menge von \(A\)-orthogonalen Vektoren mit \(d_k \in \R^n \setminus \lbrace 0\rbrace\), d.h., es gilt \(\langle d_i, A d_j \rangle = \langle A d_i, d_j \rangle = 0\) für alle \(i \neq j\).

Dann bilden die Vektoren \(d_0, \ldots, d_{n-1} \in \R^n \setminus \lbrace 0\rbrace\) eine Basis des \(\R^n\).

Proof. In den Übungsaufgaben zu zeigen. ◻

Folgender Satz zeigt uns, dass das Verfahren für eine gegebene Menge von \(A\)-konjugierten Vektoren in der Tat in \(n\) Schritten das lokalen Minimum \(x^* \in \mathbb{R}^n\) von \(f\) erreicht.

Satz 3.9 (Konvergenz des CG-Verfahrens)

Sei eine Menge von \(A\)-konjugierten Vektoren \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) mit \(d_k \in \mathbb{R}^n / \lbrace 0 \rbrace\) gegeben. Dann konvergiert das Abstiegsverfahren in konjugierten Richtungen

(3.31)\[x_{k+1} \ = \ x_k + \tau_k d_k, \qquad \tau_k \ = \ \frac{\langle r_k, d_k \rangle}{\langle d_k, Ad_k\rangle}\]

in genau \(n\) Schritten gegen die Lösung \(x^* \in \mathbb{R}^n\) des quadratischen Optimierungsproblems (3.23).

Proof. Für den Beweis der Konvergenz des Iterationsverfahrens (3.31) betrachten wir zunächst den initialen Fehler \(e_0\in \R^n\) durch die Wahl eines beliebigen Startpunktes \(x_0 \in \mathbb{R}^n\). Da wir \(A\) als positiv definit angenommen haben folgt mit Lemma 3.3, dass die Menge \(\lbrace d_k \rbrace_{k=0,\ldots,n-1}\) eine Basis des \(\mathbb{R}^n\) bildet. Daher können wir den initialen Fehler \(e_0 \in \mathbb{R}^n\) als Linearkombination in dieser Basis darstellen als:

(3.32)\[e_0 \ = \ \sum_{k=0}^{n-1} \delta_k d_k.\]

Um die unbekannten Koeffizienten \(\delta_k \in \mathbb{R}^n\) zu bestimmen können wir obige Gleichung nun jeweils von links mit einem Vektor \(d_i^\top A, i=0,\ldots,n-1\) multiplizieren und erhalten so für jeden Index eine Gleichung

\[\langle d_i^\top A, e_0 \rangle \ = \ \langle d_i^\top A, \sum_{k=0}^{n-1} \delta_k d_k \rangle \ = \ \sum_{k=0}^{n-1} \delta_k \langle d_i^\top A, d_k \rangle \ = \ \delta_i \langle d_i^\top A, d_i \rangle.\]

Hierbei haben wir die Bilinearität des Skalarproduktes in \(\mathbb{R}^n\) ausgenutzt und verwendet, dass die Vektoren \(\lbrace d_k \rbrace_{k=0,\ldots,n-1}\) konjugiert bezüglich der Matrix \(A\) sind. Damit können wir nach den unbekannten Koeffizienten \(\delta_i \in \R\) in jeder Gleichung auflösen und erhalten so einen Ausdruck für die unbekannten Koeffizienten:

\[\delta_i \ = \ \frac{\langle d_i^\top A, e_0 \rangle}{\langle d_i^\top A, d_i \rangle}.\]

Man beachte, dass dieser Ausdruck wohldefiniert ist, da wir angenommen haben, dass die Matrix \(A\) positiv definit ist. Wir addieren eine Null hinzu, indem wir Terme hinzufügen, die \(A\)-konjugiert zur Richtung \(d_i \in \R^n\) sind:

(3.33)\[\delta_i \ = \ \frac{\langle d_i^\top A, e_0 \rangle}{\langle d_i^\top A, d_i \rangle} + \underbrace{\frac{\langle d_i^\top A, \sum_{k=0}^{i-1} \tau_k d_k \rangle}{\langle d_i^\top A, d_i \rangle}}_{=\: 0} \ = \ \frac{\langle d_i^\top A, e_0 +\sum_{k=0}^{i-1} \tau_k d_k \rangle}{\langle d_i^\top A, d_i \rangle}.\]

Wir verwenden wieder den Trick, dass sich der Fehlervektor \(e_{i+1} \in \mathbb{R}^n\) entwickeln lässt zu \(e_{i+1} = e_i + \tau_i d_i\) und somit können wir rekursiv herleiten, dass

(3.34)\[e_i \ = \ e_0 + \sum_{k=0}^{i-1} \tau_k d_k.\]

Nun können wir die Gleichung (3.34) in die Darstellung der Koeffizienten \(\delta_i\) in (3.33) einsetzen und erhalten:

\[\delta_i \ = \ \frac{\langle d_i^\top A, e_0 +\sum_{k=0}^{i-1} \tau_k d_k \rangle}{\langle d_i^\top A, d_i \rangle} \ = \ \frac{\langle d_i^\top A, e_i \rangle}{\langle d_i^\top A, d_i \rangle} \ = \ \frac{\langle d_i, Ae_i \rangle}{\langle d_i^\top A, d_i \rangle} \ = \ - \frac{\langle d_i, r_i \rangle}{\langle d_i^\top A, d_i \rangle} \ = \ -\tau_i.\]

Das bedeutet, dass die Koeffizienten \(\delta_i\) in (3.33) gerade den negativen optimalen Schrittweiten \(\tau_i\) in (3.30) entsprechen, d.h., \(\delta_i\) = \(- \tau_i\). Aus der Basisdarstellung des initialen Fehlers \(e_0 = x_0 - x^*\) in (3.32) können wir somit die Behauptung des Satzes folgern:

\[x^* \ = \ x_0 - e_0 \ = \ x_0 - \sum_{k=0}^{n-1} \delta_k d_k \ = \ x_0 + \sum_{k=0}^{n-1} \tau_k d_k \ = \ x_n.\]

◻

Bemerkung 3.11 (Veränderung des Fehlers im Iterationsverfahren)

Anstatt im Beweis von Satz 3.9 zu zeigen, dass sich das eindeutige Minimum \(x^* \in \mathbb{R}^n\) durch das Iterationsverfahren zerlegen lässt, hätte man auch zeigen können, dass der Fehlervektor \(e_i \in \mathbb{R}^n\) in jedem Schritt des Iterationsverfahren kleiner wird. Es gilt nämlich nach (3.34):

\[e_i \ = \ e_0 + \sum_{k=0}^{i-1} \tau_k d_k \ = \ \sum_{k=0}^{n-1}\delta_k d_k + \sum_{k=0}^{i-1}-\delta_k d_k \ = \ \sum_{k=i}^{n-1} \delta_k d_k.\]

Man sieht also, dass für eine wachsende Anzahl an Iterationen \(i=0,\ldots,n-1\) der Fehlerterm \(e_i \in \mathbb{R}^n\) immer weniger Terme hat, bis er schlussendlich ganz verschwindet.

Außerdem sagt es uns, dass der Abstieg mit konjugierten Richtungen in dem Sinne optimal ist, als dass der Fehlerterm \(e_i = \sum_{k=i}^{n-1} \delta_k d_k\) keine Anteile der Richtungen \(\lbrace d_j \rbrace_{j=0,\ldots,i-1}\) mehr besitzt. Wir müssen also nicht mehr entlang dieser Richtungen gehen, um zum lokalen Minimum \(x^* \in \mathbb{R}^n\) von \(f\) zu gelangen. Aus Sicht der Numerik ist das eine sehr schöne Eigenschaft, da wir nicht gezwungenermaßen \(n\) Iterationen des Abstiegsverfahrens (3.31) durchführen müssen, sondern bereits nach \(k < n\) abbrechen können, um eventuell eine gute Approximation des lokalen Minimums \(x_k \approx x^* \in \mathbb{R}^n\) zu erhalten. Dies spielt insbesondere bei sehr großen Dimensionen \(n >\!\!> 1\) eine wichtige Rolle.

Beispiel 3.3 (Konjugierte Abstiegsrichtungen)

Wir wollen im Folgenden ein Beispiel zur Durchführung eines Abstiegsverfahrens mit gegebenen konjugierten Richtungen angeben. Seien folgende Werte für das lineare Gleichungssystem \(Ax = b\) gegeben:

\[\begin{split}A \ = \ \begin{pmatrix} 3 & 2\\ 2 & 6 \end{pmatrix}, \quad b \ = \ \begin{pmatrix} 2\\ -8 \end{pmatrix}.\end{split}\]

Als Startwert für unser Iterationsverfahren wählen wir \(x_0 = (-2, 2)^\top\). Wir nehmen eine Menge von zwei \(A\)-orthogonalen Vektoren \(d_0, d_1 \in \mathbb{R}^2 \setminus \lbrace 0 \rbrace\) als gegeben an mit:

\[\begin{split}d_0 \ = \ \begin{pmatrix} 0\\ -1 \end{pmatrix}, \quad d_1 \ = \ \begin{pmatrix} 3\\ -1 \end{pmatrix}.\end{split}\]

Wir sehen ein, dass die Vektoren \(d_0\) und \(d_1\) konjugiert bezüglich der Matrix \(A\) sind, denn es gilt:

\[\begin{split}\langle d_0, Ad_1 \rangle \ = \ (0,-1)\cdot \begin{pmatrix} 3 & 2\\ 2 & 6 \end{pmatrix} \begin{pmatrix} 3\\ -1 \end{pmatrix} \ = \ (0,-1) \cdot \begin{pmatrix} 7\\ 0 \end{pmatrix} \ = \ 0.\end{split}\]

Für den ersten Schritt des Iterationsverfahren berechnen wir zuerst das aktuelle Residuum

\[\begin{split}r_0 \ = \ b - Ax_0 \ = \ \begin{pmatrix} 2 \\ -8 \end{pmatrix} - \begin{pmatrix} 3 & 2\\ 2 & 6 \end{pmatrix} \begin{pmatrix} -2\\ 2 \end{pmatrix} \ = \ \begin{pmatrix} 2 \\ -8 \end{pmatrix} - \begin{pmatrix} -2 \\ 8 \end{pmatrix} \ = \ \begin{pmatrix} 4 \\ -16 \end{pmatrix}.\end{split}\]

Nun können wir die optimale Schrittweite \(\tau_0 > 0\) für den ersten Schritt durch den Ausdruck (3.30) bestimmen mit:

\[\tau_0 \ = \ \frac{\langle d_0, r_0\rangle}{\langle d_0, A d_0 \rangle} \ = \ \frac{8}{3}.\]

Hiermit können wir den ersten Abstieg durchführen und erhalten so den nächsten Iterationspunkt

\[\begin{split}x_1 \ = \ x_0 + \tau_0 d_0 \ = \ \begin{pmatrix} -2 \\ - 2/3 \end{pmatrix}.\end{split}\]

Wir wollen nur den zweiten Schritt des Verfahrens angehen und benötigen wiederum das aktuelle Residuum

\[\begin{split}r_1 \ = \ b - Ax_1 \ = \ \begin{pmatrix} 28/3 \\ 0 \end{pmatrix}.\end{split}\]

Wir berechnen wieder die neue optimale Schrittweite mittels (3.30):

\[\tau_1 \ = \ \frac{\langle d_1, r_1\rangle}{\langle d_1, A d_1 \rangle} \ = \ \frac{28}{21}.\]

Mit dieser können wir den letzten Abstiegsschritt für \(n=2\) berechnen und erhalten somit:

\[\begin{split}x_2 \ = \ x_1 + \tau_1 d_1 \ = \ \begin{pmatrix} 2\\ -2 \end{pmatrix} \ = \ x^*.\end{split}\]

Der folgende Satz hilft uns zu verstehen, warum ein Abstiegsverfahren mit konjugierten Richtungen besser funktioniert als das Gradientenabstiegsverfahren in Gradientenabstiegsverfahren.

Satz 3.10 (Orthogonalität des Residuums)

Das Residuum \(r_{i+1} = b - Ax_{i+1}\) des Abstiegsverfahren mit konjugierten Richtungen in (3.31) ist orthogonal zu allen bisherigen Abstiegsrichtungen \(d_j, j=0,\ldots,i\), d.h.

\[\langle r_{i+1}, d_j \rangle \ = \ 0, \quad \text{ für alle } j=0,\ldots,i.\]

Proof. Aus Bemerkung 3.11 wissen wir, dass wir den Fehler \(e_{i+1}\) nach \(i\) Iterationen des Abstiegsverfahrens angeben können als

\[e_{i+1} \ = \ \sum_{k=i+1}^{n-1} \delta_k d_k.\]

Wir können beide Seiten der Gleichung mit einem Zeilenvektor \(-d_j^\top A \in \R^n\) für einen Index \(0 \leq j \leq i\) von links multiplizieren und erhalten damit:

\[-\langle d_j, Ae_{i+1} \rangle \ = \ - \sum_{k=i+1}^{n-1} \delta_k \underbrace{d_j Ad_k}_{=~0} \quad \Rightarrow \ \langle d_j, r_{i+1} \rangle \ = \ 0 \quad \text{ für alle } \ 0 \leq j \leq i.\]

◻

Man beachte, dass die Eigenschaft optimal bezüglich aller vorherigen Abstiegsrichtungen nur für den Fall von konjugierten Richtungen funktioniert und nicht im Fall des Gradientenabstiegsverfahren, wie wir in Motivation gesehen haben. Hier war man nur optimal bezüglich der letzten Abstiegsrichtung und nicht bezüglich aller vorherigen Richtungen. Das resultiert in dem typischen Zickzack-Pfad beim Abstieg, wie wir ihn in fig:zig-zag gesehen haben.

3.3.5. Konjugierte Gradienten

Wir haben in Konjugierte Abstiegsrichtungen gesehen, dass wir ein iteratives Abstiegsverfahren mit konjugierten Abstiegsrichtungen \(\lbrace d_j \rbrace_{j=0,\ldots,n-1}\) verwenden können, um in \(n\) Iterationen die eindeutige Lösung des quadratischen Minimierungsproblems (3.23) und somit die Lösung des linearen Gleichungssystems \(Ax = b\) zu erhalten. Bisher sind wir jedoch davon ausgegangen, dass wir die Menge der konjugierten Vektoren \(\lbrace d_0, \ldots, d_{n-1} \rbrace\) bereits kennen. Um einen Algorithmus angeben zu können müssen wir also noch ergründen, wie sich diese Menge mit möglichst geringen numerischen Aufwand finden lässt.

Eine naheliegende Idee wäre es das Gram-Schmidtsche Orthogonalisierungsverfahren so umzugestalten, dass wir eine Menge von linear unabhängigen Vektoren \(\lbrace u_0,\ldots, u_{n-1} \rbrace\) mit \(u_k \in \R^n, k=0,\ldots,n-1\) konjugieren bezüglich der Matrix \(A\). Hierzu würde man die erste Abstiegsrichtung \(d_0 \in \R^n\) des Abstiegsverfahrens mit konjugierten Richtungen als den ersten Vektor der Menge wählen, d.h., wir setzen \(d_0 = u_0\). Anschließend konstruieren wir die nächste Abstiegsrichtung \(d_1\) indem wir alle Komponenten von \(u_1\) entfernen, die nicht \(A\)-orthogonal zu \(d_0\) sind. Für die nächste Abstiegsrichtung \(d_2\) gehen wir analog vor, nur müssen wir darauf achten alle Komponenten von \(u_2\) zu entfernen, die nicht \(A\)-orthogonal zu \(d_0\) und \(d_1\) sind. Dieses Vorgehen lässt sich iterativ bis zum Vektor \(d_{n-1}\) fortführen und man erhält eine Menge von konjugierten Vektoren \(\lbrace d_0, \ldots, d_{n-1} \rbrace\). Diese lassen sich in geschlossener Form angeben als:

(3.35)\[d_i \ = \ u_i + \sum_{k=0}^{i-1} \beta_{i,k} d_k, \quad i=1,\ldots,n-1.\]

Wir müssen jedoch die Koeffizienten \(\beta_{i,k}\) so bestimmen, dass die Vektoren \(d_i\) konjugiert zu allen vorherigen Richtungsvektoren \(d_j \in \R^n, 0 \leq j < i\) sind. Um diese Koeffizienten zu bestimmen multiplizieren wir (3.35) wieder von links mit einem Zeilenvektor \(d_j^\top A \in \R^n\) für ein \(j \in \lbrace 0,\ldots,i-1\rbrace\) und erhalten

\[\begin{split}\begin{split} &\langle d_j, A d_i \rangle \ = \ \langle d_j, A u_i \rangle + \sum_{k=0}^{i-1} \beta_{i,k} \langle d_j, A d_k \rangle \\ \Rightarrow \quad & 0 \ = \ \langle d_j, Au_i \rangle + \beta_{i,j}\langle d_j, A d_j \rangle \\\ \Rightarrow \quad & \beta_{i,j} \ = \ - \frac{\langle d_j, Au_i \rangle}{\langle d_j, Ad_j \rangle}. \end{split}\end{split}\]

Der Ausdruck für die Koeffizienten \(\beta_{i,j}\) ist wohldefiniert, da wir angenommen haben, dass \(A\) eine symmetrische, positiv definite Matrix ist. Eigentlich könnten wir jetzt zufrieden sein, da wir ein Verfahren angeben können mit dem sich ein Abstiegsverfahren mit konjugierten Richtungen konstruieren lässt.

Leider haben wir durch die Verwendung des Gram-Schmidtschen Orthogonalisierungsverfahrens nichts gewonnen, da der numerische Aufwand zur Berechnung der unbekannten Koeffizienten im besten Fall in \(\mathcal{O}(n^3)\) liegt, was genau so teuer ist wie eine Invertierung der Matrix \(A\), zum Beispiel mit dem Eliminationsverfahren von Gauss [TR].

Glücklicherweise gibt es eine Möglichkeit eine Menge von konjugierten Abstiegsrichtungen im Laufe des Iterationsverfahren (3.31) zu konstruieren ohne den numerischen Rechenaufwand des modifizierten Gram-Schmidtschen Orthogonalisierungsverfahren zu benötigen. Hierzu werden wir ähnlich mit einer Menge von initialen Richtungsvektoren \(\lbrace u_0, \ldots, u_{n-1}\rbrace\) beginnen und diese geeignet anpassen. Es stellt sich nämlich heraus, dass in jeder Iteration ein Vektor existiert, der bereits \(A\)-orthogonal zu allen vorherigen Abstiegsrichtungen ist mit Ausnahme der letzten. Diese Aussage wird durch folgendes Lemma präzisiert.

Lemma 3.4 (Eigenschaften des Residuums)

Für \(i \in \lbrace 1,\ldots,n-1 \rbrace\) befinden wir uns im \((i+1)\)-ten Schritt des Abstiegsverfahrens mit konjugierten Richtungen in (3.31) und \(\lbrace d_0, \ldots, d_{i-1} \rbrace\) ist eine Menge von \(A\)-orthogonalen Vektoren und \(u_i = r_i\) sei ein initialer Richtungsvektor für den aktuellen Iterationsschritt. Dann gilt für \(r_{i+1} = b - A x_{i+1} = b - A (x_i + \tau_i r_i)\):

\[\langle r_{i+1}, r_j \rangle \ = \ 0, \quad \text{ für alle } \ j=0,\ldots,i,\]

und außerdem auch

\[\langle r_{i+1}, Ad_j \rangle \ = \ 0, \quad \text{ für alle } \ j=0,\ldots,i-1.\]

Proof. Wir definieren zuerst die lineare Hülle, die durch die \(A\)-orthogonalen Vektoren aufgespannt wird durch:

\[\mathcal{D}_i \ \coloneqq \ \operatorname{span}\lbrace d_0, \ldots, d_{i-1} \rbrace \subset \R^n.\]

Wir gehen in diesem Beweis konstruktiv vor und wollen in jedem Schritt den unbekannten \(A\)-orthogonalen Vektor \(d_i \in \R^n\) aus dem aktuellen Residuum \(r_i \in \R^n\) und den vorigen Richtungsvektoren \(\lbrace d_0, \dots, d_{i-1}\rbrace\) konstruieren.

Wir wählen als erste Abstiegsrichtung \(d_0 = r_0\) und erhalten damit:

\[\operatorname{span}\lbrace r_0 \rbrace \ = \ \operatorname{span}\lbrace d_0 \rbrace \ = \ \mathcal{D}_1.\]

Aus Satz 3.10 wissen wir, dass \(r_1\) orthogonal zu \(d_0\) steht und somit folgt auch schon:

\[\langle r_0, r_1 \rangle \ = \ \langle d_0, r_1 \rangle \ = \ 0.\]

Da wir den nächsten \(A\)-orthogonalen Richtungsvektor \(d_1 \in \R^n\) aus dem aktuellen Residuum \(r_1\) und dem Unterraum \(\mathcal{D}_1\) konstruieren wollen, können wir damit folgern:

\[\mathcal{D}_2 \ = \ \operatorname{span}\lbrace \mathcal{D}_1, d_1 \rbrace \ = \ \operatorname{span}\lbrace \mathcal{D}_1, r_1 \rbrace \ = \ \operatorname{span}\lbrace r_0, r_1 \rbrace.\]

Analog können wir nun für beliebiges \(i \in \lbrace 1,\ldots,n-1\rbrace\) folgern, dass

\[\mathcal{D}_i \ = \ \operatorname{span}\lbrace r_0,\ldots, r_{i-1} \rbrace.\]

Aus Satz 3.10 wissen wir, dass \(r_{i+1} \perp \mathcal{D}_{i+1}\) und damit wissen wir schon, dass die erste Aussage des Satzes gilt:

\[\langle r_{i+1}, r_j \rangle \ = \ 0, \quad \text{ für alle } j=0,\ldots,i.\]

Für die zweite Aussage des Lemmas drücken wir nun das Residuum \(r_i\) durch den Fehler \(e_i\) aus und erhalten:

\[r_i \ = \ -Ae_i \ = \ -A(e_{i-1} + \tau_{i-1}d_{i-1}) \ = \ \underbrace{r_{i-1}}_{\in \mathcal{D}_i} - \underbrace{\tau_{i-1}Ad_{i-1}}_{\in A\mathcal{D}_i}.\]

Wir sehen also, dass \(r_i \in \operatorname{span}\lbrace r_{i-1}, A d_i \rbrace\). Außerdem wissen wir durch unsere Folgerungen oben, dass \(r_{i-1} \in \mathcal{D}_i\) und \(A d_{i-1} \in A\mathcal{D}_i\) gilt. Damit gilt aber schon

\[\mathcal{D}_{i+1} \ = \ \operatorname{span}\lbrace \mathcal{D}_i, r_i \rbrace \ = \ \operatorname{span}\lbrace \mathcal{D}_i, A \mathcal{D}_i \rbrace.\]

Wenn wir dies rekursiv entwickeln sehen wir interessanterweise ein, dass

\[\mathcal{D}_i \ = \ \operatorname{span}\lbrace d_0, Ad_0, \ldots, A^{i-1}d_0\rbrace\]

Nach Satz 3.10 wissen wir jedoch auch, dass \(r_{i+1} \perp \mathcal{D}_{i+1}\) und somit muss gelten auch für den Unterraum gelten \(r_{i+1} \perp A\mathcal{D}_i\). Und damit haben wir die zweite Aussage des Satzes gezeigt, nämlich dass \(r_{i+1} \perp_A \mathcal{D}_i\) oder

\[\langle r_{i+1}, A d_j \rangle, \quad \text{ für alle } j=0,\ldots,i-1.\]

◻

Lemma 3.4 sagt aus, dass das Residuum \(r_i\) ein guter Ausgangspunkt für einen weiteren \(A\)-orthogonalen Richtungsvektor \(d_i \in \R^n\) im Punkt \(x_i\) ist, da es bereits zu allen bisherigen \(A\)-orthogonalen Richtungen \(d_0,\ldots,d_{i-2}\) konjugiert bezüglich der Matrix \(A\) ist. Wir müssen also nur noch dafür sorgen, dass der initiale Richtungsvektor \(r_i \in \R^n\) \(A\)-orthogonal zur letzten Suchrichtung \(d_{i-1}\) wird. Dies ist numerisch wesentlich günstiger als einen beliebig gewählten Richtungsvektor \(u_i \in \mathbb{R}^n \setminus \lbrace 0 \rbrace\) \(A\)-orthogonal zu machen bezüglich aller vorigen Richtungsvektoren \(d_0,\ldots,d_{i-1}\).

Wir wollen also im Folgenden das vollständige Abstiegsverfahren mit konjugierten Richtungen angeben. Da die initialen Richtungen nun als \(r_i = -\nabla f(x_i)\) für alle Iterationen \(i=0,\ldots, n-1\) gewählt werden, nennt man dieses Verfahren auch das Abstiegsverfahren der konjugierten Gradienten.

Satz 3.11 (Konjugation der Residuen bezüglich A)

Wir befinden uns im \((i+1)\)-ten Schritt des Verfahrens der konjugierten Gradienten für \(i=0,\ldots, n-1\) in einem Punkt \(x_{i+1} \in \mathbb{R}^n\) und wir wählen als initialen Richtungsvektor das aktuelle Residuum \(r_{i+1} \in \R^n\), welches nach Lemma 3.4 bereits \(A\)-orthogonal zu fast allen vorherigen Richtungsvektoren \(\lbrace d_0, \ldots, d_{i-1}\rbrace\) ist. Indem wir den neuen Richtungsvektor \(d_{i+1} \in \R^n\) definieren als

(3.36)\[d_{i+1} \ \coloneqq \ r_{i+1} + \beta_{i+1} d_i, \qquad \beta_{i+1} \ \coloneqq \ \frac{\langle r_{i+1}, r_{i+1} \rangle}{\langle r_i, r_i \rangle},\]

erhalten wir die Eigenschaft, dass dieser Richtungsvektor nun \(A\)-orthogonal zu allen bisherigen Richtungsvektoren des Verfahrens ist, d.h.,

\[d_{i+1} \: \perp_A \: d_j, \quad j = 0,\ldots,i.\]

Proof. Wir müssen den initialen Richtungsvektor \(u_{i+1} \in \mathbb{R}^n \setminus \lbrace 0 \rbrace\) so modifizieren, dass der resultierende Vektor \(d_{i+1}\) konjugiert zu \(d_i\) bezüglich der Matrix \(A\) ist. Mit dem modifizierten Gram-Schmidtschen Orthogonalisierungsverfahren erhalten wir die Form

\[d_{i+1} \ = \ r_{i+1} + \beta_{i+1} d_i, \qquad \beta_{i+1} \ = \ - \frac{\langle r_{i+1}, Ad_i\rangle}{\langle d_i, Ad_i \rangle}.\]

Die neue Abstiegsrichtung \(d_{i+1} \in \mathbb{R}^n \setminus \lbrace 0 \rbrace\) ist nach Konstruktion \(A\)-orthogonal zu allen vorherigen Abstiegsrichtungen \(\lbrace d_0, \ldots, d_i \rbrace\), jedoch wollen wir den Koeffizienten \(\beta_{i+1}\) noch näher charakterisieren im Folgenden.

Aus dem Beweis von Lemma 3.4 wissen wir, dass wir das aktuelle Residuum ausdrücken können als

\[r_{i+1} \ = \ r_i - \tau_iAd_i.\]

Wir multiplizieren diese Gleichung von links mit dem Zeilenvektor \(r_{i+1}^\top \in \R^n\) und erhalten

\[\langle r_{i+1}, r_{i+1} \rangle \ = \ \langle r_{i+ 1}, r_i \rangle - \tau_i\langle r_{i+1}, Ad_i \rangle.\]

Wir wissen aus Lemma 3.4 jedoch auch, dass \(\langle r_{i+1}, r_i \rangle = 0\) gilt und damit erhalten wir den Ausdruck

\[-\frac{1}{\tau_i}\langle r_{i+1}, r_{i+1} \rangle \ = \ \langle r_{i+1}, Ad_i \rangle.\]

Wenn wir nun die optimale Schrittweite \(\tau_i = \frac{\langle r_i, d_i \rangle}{\langle d_i, Ad_i \rangle}\) aus (3.31) einsetzen erhalten wir für den Koeffizienten \(\beta_{i+1}\):

\[\begin{split}\begin{split} & \ \langle r_{i+1}, Ad_i \rangle \ = \ -\frac{1}{\tau_i}\langle r_{i+1}, r_{i+1} \rangle \ = \ -\frac{\langle d_i, Ad_i \rangle}{\langle r_i, d_i \rangle} \langle r_{i+1}, r_{i+1} \rangle \\ \Rightarrow & \ \frac{\langle r_{i+1}, r_{i+1}\rangle }{\langle r_i, d_i \rangle} \ = \ -\frac{\langle r_{i+1}, Ad_i \rangle}{\langle d_i, Ad_i \rangle} \ = \ \beta_{i+1}. \end{split}\end{split}\]

Schlussendlich können wir den Nenner in diesem Ausdruck noch weiter umschreiben, da der letzte Richtungsvektor \(d_i \in \R^n\) mit dem modifizierten Gram-Schmidt Orthogonalisierungsverfahren auch ausgedrückt werden kann als

\[\langle r_i, d_i \rangle \ = \ \langle r_i, r_i + \beta_i d_{i-1} \rangle \ = \ \langle r_i, r_i \rangle + \beta_i \underbrace{\langle r_i, d_{i-1}\rangle}_{=~0} \ = \ \langle r_i, r_i \rangle.\]

Das Skalarprodukt in obiger Gleichung verschwindet auf Grund von der Aussage von Satz 3.10 und somit erhalten wir schlussendlich für den Koeffizienten \(\beta_{i+1}\) den simplen Ausdruck:

\[\beta_{i+1} \ = \ \frac{\langle r_{i+1}, r_{i+1}\rangle}{\langle r_i, r_i \rangle}.\]

◻

Mit der Herleitung von (3.36) können wir nun einen Algorithmus für das Abstiegsverfahren mit konjugierten Gradienten zum Lösen eines Gleichungssystems \(Ax = b\) angeben.

Algorithmus 3.4 (Lineares konjugierte Gradientenverfahren)

function \(x^*=\)conjugateGradient\((A, b, x_0)\)   # Initialisierung \(d_0 = r_0 = b - Ax_0\)   # Führe genau \(n\) Schritte durch   # Berechne Schrittweite \(\tau_k \ = \ \frac{\langle r_k, d_k\rangle}{\langle d_k, A d_k\rangle}\)   # Führe Abstiegsschritt durch \(x_{k+1} \ = \ x_k + \tau_k d_k\)     # Berechne effizient neues Residuum \(r_{k+1} \ = \ r_k - \tau_k Ad_k\)   # Berechne Koeffizienten \(\beta_{k+1} \ = \ \frac{\langle r_{k+1}, r_{k+1}\rangle}{\langle r_k, r_k \rangle}\)   # Berechne neue Abstiegsrichung mit Gram-Schmidt \(d_{k+1} = r_{k+1} + \beta_{k+1} d_k\)     # Ausgabe des letzten Punktes \(x^* = x_{k+1}\)

3.3.6. Verallgemeinerung für nichtlineare Optimierung

Wir haben festgestellt, dass das konjugierte Gradientenverfahren in Algorithmus 3.4 als Minimierung der konvexen quadratischen Funktion \(f(x) \coloneqq \langle x, Ax \rangle - \langle x, b \rangle + c\) konzipiert ist, um das äquivalente lineare Gleichungssystem \(Ax = b\) zu lösen. Es ist natürlich zu fragen, ob wir diesen Ansatz anpassen können, um allgemeine konvexe Zielfunktionen oder sogar allgemeine nichtlineare Zielfunktionen \(f\) zu minimieren.

Um eine nichtlineare Variante des konjugierte Gradientenverfahrens zu erhalten, müssen wir einige Anpassungen vornehmen. Fletcher und Reeves [] zeigten, dass man anstelle der explizit berechenbaren Schrittweite \(\tau_k\) für das quadratische Probleme in (3.30) zunächst einmal eine Schrittweite wählen muss, die ein approximatives Minimum der nichtlinearen Zielfunktion \(f\) entlang der Suchrichtung \(d_k \in \R^n\) identifiziert. Darüber hinaus muss das bisher verwendete Residuum \(r_k = b - Ax_x \in \R^n\) für den Einsatz in der nichtlinearen Optimierung durch den Gradienten der Zielfunktion \(f\) ersetzt werden.

Diese beiden Änderungen führen zu dem folgenden Algorithmus für die nichtlineare Optimierung.

Algorithmus 3.5 (Nichtlineares konjugierte Gradientenverfahren)

function \(x^*=\)nonlinearConjugateGradient\((f, \nabla f, x_0, \tau_0, \sigma)\)   # Initialisierung \(d_0 = - \nabla f(x_0)\)     # Verringere Schrittweite um Faktor \(\sigma\) \(\tau_k = \sigma\tau_k\)   # Führe Abstiegsschritt durch \(x_{k+1} \ = \ x_k + \tau_k d_k\)   # Berechne Koeffizienten \(\beta_{k+1} \ = \ \frac{\langle \nabla f(x_{k+1}), \nabla f(x_{k+1})\rangle}{\langle \nabla f(x_{k}), \nabla f(x_{k}) \rangle}\)   # Berechne neue Abstiegsrichung mit Gram-Schmidt \(d_{k+1} = -\nabla f(x_{k+1}) + \beta_{k+1} d_k\)     # Ausgabe des letzten Punktes \(x^* = x_{k+1}\)

Man beachte, dass man für die Implementierung von Algorithmus 3.5 keinerlei Matrix-Vektor Multiplikationen benötigt und man lediglich Auswertungen der Zielfunktion \(f\) und ihres Gradienten \(\nabla f\) braucht. Für den Fall, dass es sich bei der Zielfunktion \(f\) um eine strikt konvexe, quadratische Funktion handelt und wir in jedem Schritt die optimale Schrittweite \(\tau_k > 0\) bestimmten können, so entspricht Algorithmus 3.5 dem linearen konjugierte Gradientenverfahren in Algorithmus 3.4.

In der Literatur hat sich unter Anderem eine Modifikation des nichtlinearen konjugierte Gradientenverfahrens von Fletcher und Reeves etabliert, die sich in numerischen Experimenten häufig als robuster und effizienter herausgestellt hat.

Bemerkung 3.12 (Polak-Ribière Variante)

Bei der sogenannten Polak-Ribière Variante des Algorithmus 3.5 unterscheidet sich hauptsächlich die Berechnung des Parameters \(\beta_{k+1} \in \R\) in (3.36) für die Anpassung der nächsten Abstiegsrichtung. Hierbei wird der Faktor \(\beta_{k+1}\) nämlich wie folgt berechnet

\[\beta^{PR}_{k+1} \ \coloneqq \ \frac{\langle \nabla f(x_{k+1}), \nabla f(x_{k+1}) - \nabla f(x_{k}) \rangle}{\langle \nabla f(x_{k}), \nabla f(x_{k}) \rangle}.\]

Man sieht ein, dass im Falle einer strikt konvexen, quadratischen Zielfunktion \(f\) mit optimaler Schrittweitenwahl für die \(\tau_k > 0\) die Orthogonalitätsbedingung für sukzessive Gradienten aus Lemma 3.2 gilt mit

\[\langle \nabla f(x_{k+1}), \nabla f(x_k) \rangle \ = \ 0.\]

In diesem Fall stimmt der Faktor \(\beta^{PR}_{k+1}\) mit dem Faktor \(\beta_{k+1}\) aus Algorithmus 3.5 überein. In allen anderen Fällen hingegen unterscheiden sich die Faktoren im Allgemeinen und führen so zu einen signifikant unterschiedlichen Konvergenzverhalten der jeweiligen Abstiegsverfahren.