Home | english  |  Impressum  |  Datenschutz  |  Sitemap  |  Intranet  |  KIT
Arbeitsgruppe: Numerische Simulation, Optimierung und Hochleistungsrechnen

Sekretariat
Kollegiengebäude Mathematik (20.30)
Zimmer 3.039

Adresse
Karlsruher Institut für Technologie (KIT)
Institut für Angewandte und Numerische Mathematik
Englerstrasse 2
76131 Karlsruhe

Öffnungszeiten:

Tel.: +49 721 608 - 42062

Fax.: +49 721 608 - 44178

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.

  1. Wie wird die Lösung konstruiert?
  2. Warum funktioniert der Ansatz überhaupt?
  3. Woher kommt die Effizienz des Verfahrens?
  4. Der vollständige Algorithmus.
  5. Woher kommt der Name?
  6. 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

Auf 64 Seiten wird das Verfahren aus dem Blick der Minimierungsaufgabe anschaulich beschrieben.
Die Darstellung hier ist kurz und bündig- zu beachten ist, dass hier per Definition eine positiv definite Matrix symmetrisch ist.
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.
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.
Der Artikel liefert eine Übersicht über das Verfahren und alternative Ansätze.