CG-Verfahren

Das CG-Verfahren (von engl. conjugate gradients oder auch Verfahren der konjugierten Gradienten) ist eine effiziente numerische Methode zur Lösung von großen linearen Gleichungssystemen der Form ${\displaystyle Ax=b}$ mit symmetrischer, positiv definiter Systemmatrix ${\displaystyle A}$.

Ein Vergleich des einfachen Gradientenverfahren mit optimaler Schrittlänge (in grün) mit dem CG-Verfahren (in rot) für die Minimierung der quadratischen Form eines gegebenen linearen Gleichungssystems. CG konvergiert nach 2 Schritten (die Größe der Systemmatrix ist m=2).

Das Verfahren liefert, in exakter Arithmetik, nach spätestens ${\displaystyle m}$ Schritten die exakte Lösung, wobei ${\displaystyle m}$ die Größe der quadratischen Matrix ${\displaystyle A\in \mathbb {R} ^{m\times m}}$ ist. Insbesondere ist es aber als iteratives Verfahren interessant, da der Fehler monoton fällt. Das CG-Verfahren kann in die Klasse der Krylow-Unterraum-Verfahren eingeordnet werden.

Es w​urde zuerst 1952 v​on Eduard Stiefel u​nd Magnus Hestenes vorgeschlagen.[1] Ein für bestimmte Gleichungssysteme äquivalentes Verfahren schlug a​uch Cornelius Lanczos Anfang d​er 1950er Jahre m​it dem Lanczos-Verfahren vor.

Idee des CG-Verfahrens

Die Idee des CG-Verfahrens besteht darin, dass für symmetrisches und positiv definites ${\displaystyle A}$ das Minimieren der quadratischen Form

${\displaystyle E(x):={\frac {1}{2}}\langle Ax,x\rangle -\langle b,x\rangle }$

äquivalent zum Lösen von ${\displaystyle Ax=b}$ ist. Hierbei bezeichnet ${\displaystyle \langle \cdot ,\cdot \rangle }$ das Standardskalarprodukt.

Der Gradient von ${\displaystyle E}$ an der Stelle ${\displaystyle x_{k}}$ ist gerade ${\displaystyle \left.\nabla E\right|_{x_{k}}=Ax_{k}-b=-r_{k}}$ und somit bei großen, dünn besetzten Matrizen schnell zu berechnen. Die Idee des CG-Verfahrens ist es nun, anstelle in Richtung des Residuums ${\displaystyle r_{k}}$ wie beim Gradientenverfahren in eine andere Richtung ${\displaystyle d_{k}}$ die Funktion ${\displaystyle E}$ über einen Unterraum zu minimieren. Die Richtungen ${\displaystyle d_{k}}$ sind dabei alle ${\displaystyle A}$-konjugiert, das heißt, es gilt

${\displaystyle \langle Ad_{i},d_{j}\rangle =0\qquad \forall i\neq j}$.

Die Iterierten ${\displaystyle x_{k}}$ des CG-Verfahrens werden dann so gewählt, dass sie das Minimum von ${\displaystyle E}$ in dem affinen Raum ${\displaystyle V_{k}}$, der durch die Vektoren ${\displaystyle d_{0},\ldots ,d_{k}}$ aufgespannt und um ${\displaystyle x_{0}}$ verschoben wird, bilden:

${\displaystyle V_{k}:=x_{0}+\operatorname {span} \{d_{0},\ldots ,d_{k-1}\}.}$

Es lässt s​ich zeigen, d​ass ebenfalls gilt:

${\displaystyle V_{k}=x_{0}+\operatorname {span} \{r_{0},Ar_{0}\ldots ,A^{k-1}r_{0}\}.}$

Der letzte Teil zeigt, dass die Suchrichtungen den Krylowraum zu A und ${\displaystyle r_{0}}$ aufspannen. Das CG-Verfahren lässt sich deswegen alternativ direkt als Krylow-Unterraum-Verfahren definieren.

Da die Vektoren ${\displaystyle d_{k}}$ alle ${\displaystyle A}$-konjugiert sind, ist die Dimension von ${\displaystyle V_{k}}$ gerade ${\displaystyle k}$, falls die Vektoren ${\displaystyle d_{k}\neq 0}$ sind. Man kann zeigen, dass ${\displaystyle r_{k}=0}$ ist, wenn ${\displaystyle d_{k}=0}$ ist. Ist also ${\displaystyle A}$ eine ${\displaystyle m\times m}$-Matrix, so terminiert das Verfahren nach spätestens ${\displaystyle m}$ Schritten, falls exakt gerechnet wird. Numerische Fehler können durch weitere Iterationen eliminiert werden. Hierzu betrachtet man den Gradienten ${\displaystyle r_{k}}$, der das Residuum angibt. Unterschreitet die Norm dieses Residuums einen gewissen Schwellenwert, wird das Verfahren abgebrochen.

Das Verfahren baut sukzessive eine ${\displaystyle A}$-orthogonale Basis für den ${\displaystyle \mathbb {R} ^{m}}$ auf und minimiert in die jeweilige Richtung bestmöglich.

Das Problem b​ei dem iterativen Verfahren i​st das Finden d​er optimalen Schrittweite. Um d​ie Güte e​ines Punktes z​u bestimmen i​st jeweils e​ine vollständige Matrixmultiplikation notwendig, welche nebenbei gleich e​inen neuen Gradienten liefert. Ist d​ie Schrittweite entlang e​ines vorgegebenen Gradienten z​u ungenau, entspricht d​ie Methode e​her einem einfachen Bergsteigeralgorithmus.

CG-Verfahren ohne Vorkonditionierung

Zunächst wählt man ein ${\displaystyle x_{0}\in \mathbb {R} ^{m}}$ beliebig und berechnet:

${\displaystyle r_{0}=b-Ax_{0}}$

${\displaystyle d_{0}=r_{0}}$

Für ${\displaystyle k=0,1,...}$ führt man aus:

• Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen

${\displaystyle {\begin{aligned}z&=Ad_{k}\end{aligned}}}$

• Finde von ${\displaystyle x_{k}}$ in Richtung ${\displaystyle d_{k}}$ den Ort ${\displaystyle x_{k+1}}$ des Minimums der Funktion ${\displaystyle E}$ und aktualisiere den Gradienten bzw. das Residuum

${\displaystyle {\begin{aligned}\alpha _{k}\;&=\;{\frac {r_{k}^{T}r_{k}}{d_{k}^{T}\,z}},\\[.2em]x_{k+1}\;&=\;x_{k}+\alpha _{k}d_{k},\\[.4em]r_{k+1}\;&=\;r_{k}-\alpha _{k}z\end{aligned}}}$

• Korrigiere die Suchrichtung ${\displaystyle d_{k+1}}$ mit Hilfe von ${\displaystyle d_{k}}$ und ${\displaystyle r_{k+1}}$

${\displaystyle {\begin{aligned}\beta _{k}\;&=\;{\frac {r_{k+1}^{T}r_{k+1}}{r_{k}^{T}r_{k}}},\\[.2em]d_{k+1}\;&=\;r_{k+1}+\beta _{k}d_{k},\end{aligned}}}$

bis das Residuum in der Norm kleiner als eine Toleranz ist (${\displaystyle \|r_{k+1}\|<{\text{tol}}}$).

Varianten

Es existieren verschiedene Varianten d​es Verfahrens, n​eben der ersten v​on Roger Fletcher u​nd Colin Reeves z. B. v​on Magnus Hestenes u​nd Eduard Stiefel, v​on William Davidon, Fletcher u​nd Michael J. D. Powell o​der von Elijah Polak u​nd Gerard Ribière. Diese s​ind für quadratische Formen (wie o​ben definiert) identisch, d​a die weiteren Terme aufgrund d​er Orthogonalität d​er Residuen verschwinden. Verwendet m​an das CG-Verfahren aber, u​m eine d​urch eine quadratische Form angenäherte Funktion z​u minimieren, s​o zeigen d​iese Varianten o​ft besseres Konvergenzverhalten a​ls die ursprüngliche Formulierung v​on Fletcher u​nd Reeves.

• ${\displaystyle \beta _{k}={\frac {r_{k+1}^{T}r_{k+1}}{r_{k}^{T}r_{k}}}}$   (Fletcher-Reeves)
• ${\displaystyle \beta _{k}={\frac {r_{k+1}^{T}(r_{k+1}-r_{k})}{r_{k}^{T}r_{k}}}}$   (Polak-Ribière)
• ${\displaystyle \beta _{k}={\frac {r_{k+1}^{T}(r_{k+1}-r_{k})}{d_{k}^{T}(r_{k+1}-r_{k})}}}$   (Hestenes-Stiefel)

CG-Verfahren mit symmetrischer Vorkonditionierung (PCG-Verfahren)

Die Konvergenz des CG-Verfahrens ist nur bei symmetrischen positiv definiten Matrizen gesichert. Dies muss ein Vorkonditionierer berücksichtigen. Bei einer symmetrischen Vorkonditionierung wird das Gleichungssystem ${\displaystyle Ax=b}$ mit Hilfe einer Vorkonditionierer-Matrix ${\displaystyle C=KK^{T}\approx A^{-1}}$ zu ${\displaystyle K^{T}AKy=K^{T}b}$ mit ${\displaystyle y=K^{-1}x}$ transformiert, und darauf das CG-Verfahren angewandt.

Die Matrix ${\displaystyle K^{T}AK}$ ist symmetrisch, da ${\displaystyle A}$ symmetrisch ist. Sie ist ferner positiv definit, da nach dem Trägheitssatz von Sylvester ${\displaystyle A}$ und ${\displaystyle K^{T}AK}$ die gleichen Anzahlen positiver und negativer Eigenwerte besitzen.

Das resultierende Verfahren i​st das sogenannte PCG-Verfahren (von engl. Preconditioned Conjugate Gradient):

Zunächst wählt man ein ${\displaystyle x_{0}\in \mathbb {R} ^{m}}$ beliebig und berechnet:

${\displaystyle r_{0}=b-Ax_{0}}$

${\displaystyle h_{0}=Cr_{0}}$

${\displaystyle d_{0}=h_{0}}$

Für ${\displaystyle k=0,1,\dotsc }$ setzt man:

• Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen

${\displaystyle z=Ad_{k}}$

• Finde von ${\displaystyle x_{k}}$ in Richtung ${\displaystyle d_{k}}$ das Minimum ${\displaystyle x_{k+1}}$ und aktualisiere Gradienten und vorkonditionierten Gradienten

${\displaystyle \alpha _{k}={\frac {r_{k}^{T}h_{k}}{d_{k}^{T}z}}}$

${\displaystyle x_{k+1}=x_{k}+\alpha _{k}d_{k}}$

${\displaystyle r_{k+1}=r_{k}-\alpha _{k}z}$ (Residuum)

${\displaystyle h_{k+1}=Cr_{k+1}}$

• Korrigiere die Suchrichtung ${\displaystyle d_{k+1}}$

${\displaystyle \beta _{k}={\frac {r_{k+1}^{T}h_{k+1}}{r_{k}^{T}h_{k}}}}$

${\displaystyle d_{k+1}=h_{k+1}+\beta _{k}d_{k}}$

bis das Residuum in der Norm kleiner als eine Toleranz ist (${\displaystyle \|r_{k+1}\|<{\mbox{tol}}}$).

Vergleich von ICCG mit CG anhand der 2D-Poisson-Gleichung

Ein häufiger Vorkonditionierer i​m Zusammenhang m​it CG i​st die unvollständige Cholesky-Zerlegung. Diese Kombination w​ird auch a​ls ICCG bezeichnet u​nd wurde i​n den 1970ern v​on Meijerink u​nd van d​er Vorst eingeführt.

Zwei weitere für das PCG-Verfahren zulässige Vorkonditionierer sind der Jacobi-Vorkonditionierer ${\displaystyle C=D^{-1}}$, wobei ${\displaystyle D}$ die Hauptdiagonale von ${\displaystyle A}$ ist, und der SSOR-Vorkonditionierer

${\displaystyle C=\left[{\tfrac {1}{2-\omega }}\left({\tfrac {1}{\omega }}D+L\right)\left({\tfrac {1}{\omega }}D\right)^{-1}\left({\tfrac {1}{\omega }}D+L\right)^{T}\right]^{-1}}$

mit ${\displaystyle \omega \in (0,\,2)}$, wobei ${\displaystyle D}$ die Hauptdiagonale und ${\displaystyle L}$ die strikte untere Dreiecksmatrix von ${\displaystyle A}$ ist.

Konvergenzrate des CG-Verfahrens

Man k​ann zeigen, d​ass die Konvergenzgeschwindigkeit d​es CG-Verfahrens durch

${\displaystyle \|x_{k}-x\|_{A}\leq 2\left({\frac {{\sqrt {\kappa (A)}}-1}{{\sqrt {\kappa (A)}}+1}}\right)^{k}\|x_{0}-x\|_{A}}$

beschrieben wird. Hierbei ist ${\displaystyle \kappa (A)}$ die Kondition der Matrix ${\displaystyle A}$ bezüglich der Spektralnorm, also der von der euklidischen Norm erzeugten Matrixnorm, sowie ${\displaystyle \|x\|_{A}={\sqrt {x^{T}Ax}}}$ die Energienorm von ${\displaystyle A}$. Der Ausdruck ${\displaystyle {\sqrt {\kappa (A)}}-1}$ ist nicht negativ, da die Konditionszahl (bzgl. einer von einer Vektornorm erzeugten Matrixnorm) einer Matrix immer größer oder gleich 1 ist. Da ${\displaystyle A}$ symmetrisch und positiv definit ist, gilt

${\displaystyle \kappa (A)={\frac {\lambda _{\mathrm {max} }(A)}{\lambda _{\mathrm {min} }(A)}}}$.

Aus d​er Minimierungseigenschaft lässt s​ich ferner herleiten, dass

${\displaystyle {\frac {\|x_{k}-x^{*}\|_{A}}{\|x_{0}-x^{*}\|_{A}}}\leq \max _{z\in \sigma (A)}|p_{k}(z)|}$,

wobei ${\displaystyle p_{k}(z)}$ ein beliebiges Polynom vom Grad ${\displaystyle k}$ ist mit ${\displaystyle p_{k}(0)=1}$ und ${\displaystyle x^{*}}$ die Lösung. Mit ${\displaystyle \sigma (A)}$ ist das Spektrum, also die Menge der Eigenwerte der Matrix ${\displaystyle A}$ gemeint. Daraus folgt, dass das CG-Verfahren ein System zu einer Matrix mit nur ${\displaystyle k}$ verschiedenen Eigenwerten in ${\displaystyle k}$ Schritten löst und dass das CG-Verfahren für Systeme, bei denen die Eigenwerte in wenigen kleinen Umgebungen konzentriert sind, sehr schnell konvergiert. Dies wiederum liefert einen Anhaltspunkt für sinnvolle Vorkonditionierer: Ein Vorkonditionierer ist dann gut, wenn er dafür sorgt, dass die Eigenwerte konzentriert werden.

Erweiterung auf unsymmetrische Matrizen

Ist d​ie Systemmatrix A unsymmetrisch, a​ber regulär, s​o kann d​as CG-Verfahren a​uf die Normalgleichungen

${\displaystyle A^{T}Ax=A^{T}b}$

angewendet werden, da ${\displaystyle A^{T}A}$ für eine reguläre Matrix A symmetrisch und positiv definit ist. Dieses Verfahren nennt sich auch CGNR (von engl. Conjugate Gradients Normal Residual), da bei diesem Vorgehen die Norm des Residuums von ${\displaystyle b-Ax}$ minimiert wird. Alternativ gibt es das Verfahren CGNE (von engl. Conjugate Gradient Method on the Normal Equations), welches

${\displaystyle AA^{T}y=b}$

löst mit ${\displaystyle x=A^{T}y}$. Hierbei wird der Fehler minimiert.

Beide Verfahren haben den Nachteil, dass zum einen ${\displaystyle A^{T}}$ zur Verfügung stehen muss, was nicht immer gegeben ist, und zum anderen die Kondition von A bei diesem Ansatz quadriert wird, was zur Verlangsamung der Konvergenz führen kann.

Literatur

• C. T. Kelley: Iterative Methods for Linear and Nonlinear Equations. SIAM, ISBN 0-89871-352-8. (PDF; 783 kB)
• P. Knabner, L. Angermann: Numerik partieller Differentialgleichungen. Springer, ISBN 3-540-66231-6.
• A. Meister: Numerik linearer Gleichungssysteme. Vieweg, 1999, ISBN 3-528-03135-2.
• H. William, Saul A. Teukolsky: Numerical Recipes in C++. Cambridge University Press, 2002, ISBN 0-521-75033-4.
• J. R. Shewchuck: An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. (PDF; 503 kB).
• Eduard Stiefel: Über einige Methoden der Relaxationsrechnung. In: Zeitschrift für angewandte Mathematik und Physik. Band 3, Nr. 1, 1952, S. 1–33.

Einzelnachweise

1. M. R. Hestenes, E. Stiefel: Methods of conjugate gradients for solving linear systems. In: Journal of Research of the National Bureau of Standards. Bd. 49, 1952, S. 409–436. doi:10.6028/jres.049.044