$Webrelaunch 2020$

Fakultät für Mathematik

Startseite

Institute

Institut für Angewandte und Numerische Mathematik

Arbeitsgruppe 3: Wissenschaftliches Rechnen

Warum funktioniert das CG-Verfahren? - S. Ritterbusch

Sei $A$ eine symmetrische positiv definite $n\times n$ -Matrix, d.h. $A^\top=A$ und für alle $x\in V\backslash\{0\}$ gilt $x^\top A x>0$ .

Wir suchen eine Lösung $x\in V$ im $n$ -dimensionalen Vektorraums $V$ des Gleichungssystems

$Ax=c\qquad\qquad (1)$

für ein $c\in V$ .

Diese Lösung kann man mit Hilfe des CG-Verfahrens (oder conjugate gradients method, bzw. konjugierte Gradienten Verfahren) bestimmen. Häufig wird das Verfahren aus der (namensgebenden) Sicht einer Minimierung einer quadratischen Form beschrieben, im Gegensatz dazu wird es hier aus der Sicht der Basisdarstellung beschrieben, die im Algorithmus schrittweise gebildet wird. Dieser Zugang hat den Vorteil, dass die einzelnen Schritte und die Eigenschaften konstruktiv leichter ersichtlich werden, statt wie sonst oft üblich einen gegebenen Algorithmus im Nachhinein nur zu belegen.

Wie wird die Lösung konstruiert?
Warum funktioniert der Ansatz überhaupt?
Woher kommt die Effizienz des Verfahrens?
Der vollständige Algorithmus.
Woher kommt der Name?
Weitere Darstellungen.

Eine Basisdarstellung

Die Grundidee des Verfahrens ist, dass wir ausnutzen, dass eine positiv definite symmetrische Matrix ein Skalarprodukt auf dem Vektorraum definiert:

$\langle x, y\rangle_A:=x^\top A y\,,\qquad x\perp_A y\ :\Leftrightarrow\ x^\top Ay=0$

Mit Hilfe des Gram-Schmidtschen Orthogonalisierungsverfahren kann man somit eine $A$ -orthogonale Basis $B=\left(b_1,\ldots,b_n\right)$ , d.h. $b_k^\top Ab_l=0$ wenn $k\not=l$ , des $n$ -dimensionalen Vektorraums $V$ erzeugen und die Lösung in dieser Basis mit Koeffizienten $\lambda_j\in\mathbb{R}$ darstellen:

$x=\sum_{j=1}^n \lambda_j b_j$

Multipliziert man die zu lösende Gleichung $(1)$ von links mit den Basisvektoren so erhalten wir die Koeffizienten sehr schnell dank der $A$ -Orthogonalität der Basis $B$ :

$b_k^\top A x=b_k^\top A \sum_{j=1}^n \lambda_j b_j= \lambda_k b_k^\top A b_k\stackrel!=b_k^\top c$

Die Koeffizienten sind also durch $\lambda_k=\frac{b_k^\top c}{b_k^\top A b_k}$ bestimmt. Im CG-Verfahren berechnet man die $\lambda_k$ und orthogonale $b_k$ schrittweise und erhält jeweils eine neue Näherung $x_k$ an die Lösung:

$x_k:=\sum_{j=1}^k \lambda_k b_k$

Spätestens im $n$ -ten Schritt ist die Lösung $x=x_n$ erreicht und wegen der $A$ -Orthogonalität der Basis, sind die Näherungen als $A$ -orthogonale Projektionen der Lösung $x$ auf die Untervektorräume

$K_k:=\mathrm{span}\{b_1,\ldots,b_k\}$

in jedem Schritt in diesem Sinne optimal. Es bleibt nur die Frage, wie man die Basis $B$ sinnvoll konstruiert, um Arbeit zu sparen, denn eine vollständige Gram-Schmidtsche Orthogonalisierung ist genauso aufwendig, wie das LGS direkt zu lösen.

Der Residuenansatz

Das (namensgebende, s.u.) Grundprinzip lautet die Basisvektoren iterativ orthogonalisiert aus den Residuen

$r_k:=Ax_k-c$

zu bestimmen. Startet man dabei mit $x_0=0$ erhält man damit als Lösungen die Projektionen in die sogenannten Krylow-Räume

$K_k=\mathrm{span}\{c,Ac,A^2c\ldots,A^{k-1}c\}\,.$

Die wichtigste Eigenschaft dieses Krylow-Ansatzes ist, dass für alle $y\in K_k$ und $k<n$ automatisch $Ay\in K_{k+1}$ gilt- dies werden wir im Beweis des zweiten Lemmas benötigen.

Wir können aber jetzt schon aus den bisherigen Überlegungen den naheliegenden $k$ -ten Schritt ( $k=1,\ldots,n$ ) des (noch ineffizienten) Algorithmus so formulieren:

$\begin{array}{lrcll}\textrm{Residuum:}&r_{k-1}&=&Ax_{k-1}-c\not=0 &\qquad (2)\textrm{Orthogonalisierter Basisvektor:}&b_k&=& \displaystyle r_{k-1}-\sum_{j=1}^{k-1} \frac{r_{k-1}^\top A b_j}{b_j^\top A b_j}b_j&\qquad (3)\textrm{Darstellung in Basis:}& \lambda_{k}&=&\displaystyle \frac{b_k^\top c}{b_k^\top A b_k}\textrm{N\"{a}chste N\"aherung:}&x_k&=&x_{k-1}+\lambda_k b_k\end{array}$

Wenn der Ansatz mit den Residuen funktioniert, erhält man so nach obigen Überlegungen nach spätestens $n$ Schritten die Lösung, oder wenn in Gleichung $(2)$ ein $r_{k-1}=0$ ist. Es ist aber noch offen, ob die Residuen wirklich ein geeigneter Ansatz für die Basisvektoren sind, denn nur wenn

$r_{k-1}\not\in\mathrm{span}\{b_1,\ldots,b_{k-1}\}=K_{k-1}$

ist, so wird man mit Gleichung $(3)$ einen neuen Basisvektor erhalten. Diese Aussage können wir induktiv beweisen:

Lemma 1: Im $k$ -ten Schritt ist $r_{k-1}$ entweder der Nullvektor, oder linear unabhängig zu den vorangegangenen Residuen bzw. den Basisvektoren.

Beweis: Im ersten Schritt ist die Aussage automatisch erfüllt. Sei nun $k\leq n$ und $r_{k-2}\not=0$ und linear unabhängig zu den Basisvektoren $b_1,\ldots,b_{k-2}$ . Angenommen es gelte $0\not=r_{k-1}\in K_{k-1}$ . Dann müsste auch $0\not=A^{k-1}c\in K_{k-1}$ sein (!), bzw. es gäbe $\mu_j$ mit

$A^{k-1}c=\sum_{j=0}^{k-1}\mu_j A^j c\,.$

Wendet man auf beiden Seiten $A^{-1}$ an, so erhalten wir entweder, dass $A^{-1}c\in K_{k-1}$ war, und somit $r_{k-1}=0$ gewesen sein muss, oder dass entgegen der Annahme $A^{k-2}c\in K_{k-2}$ war, das ist ein Widerspruch. Also ist $r_{k-1}=0$ oder $r_{k-1}\not\in K_{k-1}$ . $\square$

Damit ist sichergestellt, dass der Ansatz mit den Residuen zumindest funktioniert. Diese Aussage kann man auch sehr elegant aus dem Satz von Cayley und Hamilton ableiten: Nach dem Satz kann man die Inverse einer regulären Matrix durch ein Polynom $n-1$ -ten Grades der Matrix darstellen. Da aber die Lösung gerade $x=-A^{-1}c$ ist, heisst das, dass $x$ als Linearkombination von den Vektoren $c,Ac,\ldots,A^{n-1}c$ dargestellt werden kann, also in einem der Krylow-Räume liegen wird. Diese Vorgehensweise wird oft als Motivation für den Krylow-Ansatz dargestellt.

Ein-Schritt Gram-Schmidt

Die Orthogonalisierung in der Gleichung $(3)$ wird mit jedem Schritt und Basisvektor aufwendiger und insgesamt langsamer als sogar eine Lösung z.B. mit Hilfe einer $LR$ -Zerlegung. Der erstaunliche Vorteil des Residuenansatzes ist aber, dass statt dem in Gleichung $(3)$ beschriebenen Verfahren, nur gegenüber dem letzten Basisvektor orthogonalisiert werden muss, also überraschenderweise

$b_k=r_{k-1}-\sum_{j=1}^{k-1} \frac{r_{k-1}^\top A b_j}{b_j^\top A b_j}b_j=\underline{r_{k-1}-\frac{r_{k-1}^\top A b_{k-1}}{b_{k-1}^\top A b_{k-1}}b_{k-1}}\qquad (3')$

gilt! Damit bleibt der Aufwand für jeden Schritt gleich und als iteratives Verfahren erhalten wir fortlaufend verbesserte Näherungen an die Lösung- ein sehr großer Vorteil gegenüber anderen Verfahren, da man bei sehr großen Gleichungssystem sich schon mit der Näherung nach wenigen Iterationen zufrieden gibt.

Es bleibt zu zeigen, dass $r_k\perp_A K_{k-1}$ , bzw. $b_l^\top A r_k=0$ für alle $l<k$ gilt: Das ist die entscheidende Aussage für die Effizienz des Verfahrens:

Lemma 2: Ist $r_k\not=0$ , so ist $b_l^\top A r_k=0$ , wenn $l<k$ .

Beweis: Für $k=0,1$ ist die Aussage erfüllt. Für $k\geq 2$ nutzen wir die Darstellung

$r_k=Ax_k-c=A\left(\sum_{l=1}^k \lambda_l b_l\right)-c=-c+\sum_{l=1}^k\lambda_lAb_l$

und sehen, dass damit $b_m^\top r_k=0$ , bzw. $b_m\perp r_k$ , für alle $m\leq k$ gilt, denn dann ist wegen der $A$ -Orthogonalität der Basis

$b_m^\top r_k=-b_m^\top c+\sum_{l=1}^k\lambda_lb_m^\top Ab_l=-b_m^\top c +\underbrace{\lambda_m}_{=\frac{b_m^\top c}{b_m^\top A b_m}} b_m^\top A b_m=-b_m^\top c+b_m^\top c=0\,.$

In anderen Worten, $r_k$ ist (herkömmlich) orthogonal auf allen Räumen $K_m$ , $m\leq k$ . Aus der Konstruktion bzw. Eigenschaft der Krylovräume, dass mit $b_l\in K_l$ automatisch $Ab_l\in K_{l+1}$ ist, folgt $Ab_l\perp r_k$ für alle $l+1\leq k$ , bzw. $l<k$ , d.h. $r_k^\top A b_l=0$ für $l<k$ , was zu zeigen war. $\square$

Der Algorithmus

Nach den obigen Überlegungen kann man den Algorithmus so formulieren:

$x_0=0$

Für $k=1,\ldots,n$ :

$\begin{array}{crcll}\mathrm{(I)}\qquad&r_{k-1}&=&Ax_{k-1}-c &\qquad \textrm{Ende, wenn }r_{k-1}=0\,.\mathrm{(II)}\qquad& b_k&=& \displaystyle r_{k-1}- \frac{r_{k-1}^\top A b_{k-1}}{b_{k-1}^\top A b_{k-1}}b_{k-1}\mathrm{(III)}\qquad& \lambda_{k}&=&\displaystyle \frac{b_k^\top c}{b_k^\top A b_k}\mathrm{(IV)}\qquad& x_k&=&x_{k-1}+\lambda_k b_k&\qquad \textrm{Ende, wenn }k=n\,.\end{array}$

Es gibt verschiedene Optimierungsmöglichkeiten: Im Schritt $\mathrm{(I)}$ kann man z.B. die Identität $r_{k-1}=r_{k-2}+\lambda_{k-1}Ab_{k-1}$ verwenden, und Schritte $\mathrm{(II)}$ und $\mathrm{(III)}$ vereinfachen, wenn man grundsätzlich die Basis zusätzlich auch bezüglich $A$ normiert- damit kann man sich fast alle Matrix-Multiplikationen sparen. Es gibt auch unterschiedliche Ansätze für den Schritt $\mathrm{(III)}$ . Da die Konvergenzgeschwindigkeit in einem Schritt von der Kondition der Matrix $A$ abhängt, ist es sinnvoll Vorkonditionierer zu verwenden, die die Konditionszahl verbessern.

Allgemeiner ist das Verfahren auch auf nicht symmetrische Matrizen anwendbar: Hat man eine invertierbare Matrix $B$ und will das Problem $Bx=d$ lösen, so kann man es auf $B^\top B=B^\top d$ umformen, kann es also mit symmetrischen $A=B^\top B$ und $c=B^\top d$ auf geeignete Form für das CG-Verfahren bringen.

Die Namensgebung

Der Name des Verfahrens stammt daher, dass die Lösung von $Ax=c$ genau die Minimumsstelle der quadratischen Form

$F(x)=\frac12 x^\top A x-x^\top c$

entspricht, denn es ist die Nullstelle des Gradienten $\nabla F(x)=Ax-c$ von $F$ . Jetzt ist auch klar, um welche Gradienten es geht, denn die $r_k$ sind die Gradienten $r_k=\nabla F(x_k)$ , und der Begriff der Konjugation bezeichnet die Orthogonalisierung, die zu den Basisvektoren, oder geometrisch gesehen Abstiegsrichtungen, führt.

In dieser geometrischen Interpretation liegt der Vergleich mit anderen Abstiegsverfahren, wie z.B. dem steilsten Abstieg, nahe- doch hinkt der Vergleich auf den ersten Blick, da das CG-Verfahren im Gegensatz zu den anderen Methoden theoretisch nach endlichen Schritten die exakte Lösung liefert. Auf den zweiten Blick ist der Vergleich bei wirklich großen Problemen berechtigt, denn dann wird man das CG-Verfahren nie bis zum Ende iterieren, sondern sich mit Lösungen mit hinreichend kleinem Residuum begnügen.

Weitere Quellen

J. R. Shewchuk: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain (pdf). CMU, 1994.

Auf 64 Seiten wird das Verfahren aus dem Blick der Minimierungsaufgabe anschaulich beschrieben.

T. Arens, F. Hettlich, Ch. Karpfinger, U. Kockelkorn, K. Lichtenegger, H. Stachel: Mathematik. Spektrum Akademischer Verlag, Heidelberg, 2008.

Die Darstellung hier ist kurz und bündig- zu beachten ist, dass hier per Definition eine positiv definite Matrix symmetrisch ist.

P. Deufelhard, A. Hohmann: Numerische Mathematik I, Eine algorithmisch orientierte Einführung. De Gruyter, 2008.

Die sehr verständliche Darstellung wird durch das Fehlerminimierungsproblem innerhalb der Krylow-Räume motiviert und erklärt auch die Sichtweise des Iterationsschritts als Variationsproblem.

A. Kirsch: An Introduction to the Mathematical Theory of Inverse Problems. Springer, 1996.

Die CG-Methode wird hier für beschränkte, injektive, lineare Operatoren auf einem Hilbertraum vorgestellt, insbesondere auch die optimale Stopstrategie für gestörte Daten.

CG-Verfahren in der deutschen Wikipedia

Der Artikel liefert eine Übersicht über das Verfahren und alternative Ansätze.