imported>ConGreif: /* Erweiterung auf unsymmetrische Matrizen */

2026-04-30T10:33:44Z

Erweiterung auf unsymmetrische Matrizen

Neue Seite

[[Datei:Conjugate gradient illustration.svg|mini|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 '''CG-Verfahren''' (von engl. '''c'''onjugate '''g'''radients oder auch '''Verfahren der konjugierten Gradienten''') ist eine effiziente [[Numerik|numerische]] Methode zur Lösung von großen [[Lineares Gleichungssystem|linearen Gleichungssystemen]] der Form <math>Ax=b</math> mit [[Symmetrische Matrix|symmetrischer]], [[Definitheit|positiv semidefiniter]] Systemmatrix <math>A</math>.

Das Verfahren liefert, in exakter Arithmetik, nach spätestens <math>m</math> Schritten die exakte Lösung, wobei <math>m</math> die Größe der quadratischen Matrix <math>A\in\R^{m\times m}</math> ist. In der Praxis hat es sich allerdings für große Probleme als [[iteratives Verfahren]] etabliert, denn meist reichen wenige Iterationen aus für eine gute Annäherung.
Das CG-Verfahren ist der prominenteste Vertreter der Klasse der [[Krylow-Unterraum-Verfahren]].

Es wurde zuerst 1952 von [[Eduard Stiefel]] und [[Magnus Hestenes]] vorgeschlagen.<ref>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]]</ref> Ein für bestimmte Gleichungssysteme äquivalentes Verfahren schlug auch [[Cornelius Lanczos]] Anfang der 1950er Jahre mit dem [[Lanczos-Verfahren]] vor.

== Idee des CG-Verfahrens ==
Die Idee des CG-Verfahrens besteht darin, dass für symmetrisches und positiv definites <math>A</math> das Minimieren der [[Quadratische Form|quadratischen Form]]
:<math>E(x):=\frac12\langle Ax,x\rangle - \langle b,x\rangle</math>
äquivalent zum Lösen von <math>Ax=b</math> ist. Hierbei bezeichnet <math>\langle \cdot,\cdot \rangle</math> das [[Standardskalarprodukt]].

Der [[Gradient (Mathematik)|Gradient]] von <math>E</math> an der Stelle <math>x_k</math> ist gerade <math>\left. \nabla E\right|_{x_k}=Ax_k-b=-r_k</math> und somit bei großen, [[Dünnbesetzte Matrix|dünn besetzten Matrizen]] schnell zu berechnen. Die Idee des CG-Verfahrens ist es nun, anstelle in Richtung des [[Residuum (Numerische Mathematik)|Residuums]] <math>r_k</math> wie beim [[Gradientenverfahren]] in eine andere Richtung <math>d_k</math> die Funktion <math>E</math> über einen Unterraum zu minimieren. Die Richtungen <math>d_k</math> sind dabei alle <math>A</math>-[[Konjugierte Richtungen|konjugiert]], das heißt, es gilt
:<math>\langle Ad_i,d_j\rangle=0\qquad\forall i\neq j</math>.
Die Iterierten <math>x_k</math> des CG-Verfahrens werden dann so gewählt, dass sie das Minimum von <math>E</math> in dem [[Affiner Raum|affinen Raum]] <math>V_k</math>, der durch die Vektoren <math>d_0,\ldots,d_k</math> aufgespannt und um <math>x_0</math> verschoben wird, bilden:
:<math>V_k:=x_0+\operatorname{span}\{d_0,\ldots,d_{k-1}\}.</math>
Es lässt sich zeigen, dass ebenfalls gilt:
:<math>V_k = x_0+\operatorname{span}\{r_0, Ar_0\ldots,A^{k-1}r_0\}.</math>
Der letzte Teil zeigt, dass die Suchrichtungen den [[Krylowraum]] zu ''A'' und <math>r_0</math> aufspannen. Das CG-Verfahren lässt sich deswegen alternativ direkt als Krylow-Unterraum-Verfahren definieren.

Da die Vektoren <math>d_k</math> alle <math>A</math>-konjugiert sind, ist die [[Dimension (Mathematik)|Dimension]] von <math>V_k</math> gerade <math>k</math>, falls die Vektoren <math>d_k\neq 0</math> sind. Man kann zeigen, dass <math>r_k= 0</math> ist, wenn <math>d_k= 0</math> ist. Ist also <math>A</math> eine <math>m\times m</math>-Matrix, so terminiert das Verfahren nach spätestens <math>m</math> Schritten, falls exakt gerechnet wird. Numerische Fehler können durch weitere Iterationen eliminiert werden. Hierzu betrachtet man den Gradienten <math>r_k</math>, der das Residuum angibt. Unterschreitet die [[Norm (Mathematik)|Norm]] dieses Residuums einen gewissen Schwellenwert, wird das Verfahren abgebrochen.

Das Verfahren baut sukzessive eine <math>A</math>-orthogonale Basis für den <math>\mathbb R^m</math> auf und minimiert in die jeweilige Richtung bestmöglich.

Das Problem bei dem iterativen Verfahren ist das Finden der optimalen Schrittweite. Um die Güte eines Punktes zu bestimmen, ist jeweils eine vollständige [[Matrixmultiplikation]] notwendig, welche nebenbei gleich einen neuen Gradienten liefert. Ist die Schrittweite entlang eines vorgegebenen Gradienten zu ungenau, entspricht die Methode eher einem einfachen [[Bergsteigeralgorithmus]].

== CG-Verfahren ohne Vorkonditionierung ==

Zunächst wählt man ein <math>x_0 \in \mathbb{R}^m</math> beliebig und berechnet:

:<math>r_0 = b - A x_0</math>
:<math>d_0 = r_0 </math>

Für <math>k = 0,1,...</math> führt man aus:
* Speichere [[Matrix-Vektor-Produkt]], um es nur einmal auszurechnen
::<math>\begin{align}
z&=Ad_k
\end{align}</math>
* Finde von <math>x_k</math> in Richtung <math>d_k</math> den Ort <math>x_{k+1}</math> des Minimums der Funktion <math>E</math> und aktualisiere den Gradienten bzw. das Residuum
::<math>\begin{align}
\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{align}</math>
* Korrigiere die Suchrichtung <math>d_{k+1}</math> mit Hilfe von <math>d_k</math> und <math>r_{k+1}</math>
::<math>\begin{align}
\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{align}</math>
bis das Residuum in der A-Norm kleiner als eine Toleranz ist (<math>\|r_{k+1}\| _A<\text{tol}</math>).

=== Varianten ===
Es existieren verschiedene Varianten des Verfahrens, neben der ersten von [[Roger Fletcher (Mathematiker)|Roger Fletcher]] und [[Colin Reeves]] z. B. von [[Magnus Hestenes]] und [[Eduard Stiefel]], von [[William Davidon]], Fletcher und [[Michael J. D. Powell]] oder von Elijah Polak und Gerard Ribière. Diese sind für quadratische Formen (wie oben definiert) identisch, da die weiteren Terme aufgrund der Orthogonalität der Residuen verschwinden. Verwendet man das CG-Verfahren aber, um eine durch eine quadratische Form angenäherte Funktion zu minimieren, so zeigen diese Varianten oft besseres Konvergenzverhalten als die ursprüngliche Formulierung von Fletcher und Reeves.
* <math>\beta_{k} = \frac{r_{k+1}^T r_{k+1}}{r_k^T r_k}</math>   (Fletcher-Reeves)
* <math>\beta_{k} = \frac{r_{k+1}^T (r_{k+1}-r_k)}{r_k^T r_k}</math>   (Polak-Ribière)
* <math>\beta_{k} = \frac{r_{k+1}^T (r_{k+1}-r_k)}{d_k^T (r_{k+1}-r_k)}</math>   (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 [[Vorkonditionierung|Vorkonditionierer]] berücksichtigen. Bei einer symmetrischen Vorkonditionierung wird das Gleichungssystem <math>Ax=b</math> mit Hilfe einer Vorkonditionierer-Matrix <math>C=KK^T\approx A^{-1}</math> zu <math>K^TAKy=K^Tb</math> mit <math>y=K^{-1}x</math> transformiert, und darauf das CG-Verfahren angewandt.

Die Matrix <math>K^TAK</math> ist symmetrisch, da <math>A</math> symmetrisch ist. Sie ist ferner positiv definit, da nach dem [[Trägheitssatz von Sylvester]] <math>A</math> und <math>K^TAK</math> die gleichen Anzahlen positiver und negativer [[Eigenwert]]e besitzen.

Das resultierende Verfahren ist das sogenannte PCG-Verfahren (von engl. '''P'''reconditioned '''C'''onjugate '''G'''radient):

Zunächst wählt man ein <math>x_0 \in \mathbb{R}^m</math> beliebig und berechnet:
:<math>r_0 = b - A x_0</math>
:<math>h_0 = C r_0</math>
:<math>d_0 = h_0</math>

Für <math>k = 0,1,\dotsc</math> setzt man:
* Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen
::<math>z=Ad_k</math>
* Finde von <math>x_k</math> in Richtung <math>d_k</math> das Minimum <math>x_{k+1}</math> und aktualisiere Gradienten und vorkonditionierten Gradienten
::<math>\alpha_k=\frac{r_k^T h_k}{d_k^T z}</math>
::<math>x_{k+1}=x_k+\alpha_k d_k</math>
::<math>r_{k+1}=r_k-\alpha_k z</math> ([[Residuum (Numerische Mathematik)|Residuum]])
::<math>h_{k+1}=C r_{k+1}</math>

* Korrigiere die Suchrichtung <math>d_{k+1}</math>
::<math>\beta_k=\frac{r_{k+1}^T h_{k+1}}{r_k^T h_k}</math>
::<math>d_{k+1}=h_{k+1}+\beta_k d_k</math>

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

[[Datei:ICCG-CG-comparison.png|mini|Vergleich von ICCG mit CG anhand der 2D-[[Poisson-Gleichung]]]]
Ein häufiger Vorkonditionierer im Zusammenhang mit CG ist die [[unvollständige Cholesky-Zerlegung]]. Diese Kombination wird auch als ICCG bezeichnet und wurde in den 1970ern von Meijerink und [[Henk van der Vorst|van der Vorst]] eingeführt.

Zwei weitere für das PCG-Verfahren zulässige Vorkonditionierer sind der [[Jacobi-Verfahren|Jacobi]]-Vorkonditionierer <math>C=D^{-1}</math>, wobei <math>D</math> die [[Hauptdiagonale]] von <math>A</math> ist, und der [[SSOR-Verfahren|SSOR]]-Vorkonditionierer
:<math>
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}
</math>
mit <math>\omega \in (0, \,2)</math>, wobei <math>D</math> die Hauptdiagonale und <math>L</math> die strikte untere Dreiecksmatrix von <math>A</math> ist.

== Konvergenzrate des CG-Verfahrens ==

Man kann zeigen, dass die Konvergenzgeschwindigkeit des CG-Verfahrens durch
:<math>\|x_k-x\|_A \le 2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^k\|x_{0}-x\|_A</math>
beschrieben wird. Hierbei ist <math>\kappa(A)</math> die [[Kondition (Mathematik)|Kondition]] der Matrix <math>A</math> bezüglich der [[Spektralnorm]], also der von der euklidischen Norm erzeugten Matrixnorm, sowie
<math>\|x\|_A = \sqrt{x^T A x}</math> die Energienorm von <math>A</math>. Der Ausdruck <math>\sqrt{\kappa(A)}-1</math> ist nicht negativ, da die Konditionszahl (bzgl. einer von einer Vektornorm erzeugten Matrixnorm) einer Matrix immer größer oder gleich 1 ist. Da <math>A</math> symmetrisch und positiv definit ist, gilt
:<math>\kappa(A) = \frac{\lambda_\mathrm{max}(A)}{\lambda_\mathrm{min}(A)}</math>.

Aus der Minimierungseigenschaft lässt sich ferner herleiten, dass

:<math>\frac{\|x_k-x^*\|_A}{\|x_0-x^*\|_A} \leq \max_{z \in \sigma(A)}|p_k(z)|</math>,

wobei <math>p_k(z)</math> ein beliebiges [[Polynom]] vom Grad <math>k</math> ist mit <math>p_k(0)=1</math> und <math>x^*</math> die Lösung. Mit <math>\sigma(A)</math> ist das [[Spektrum (lineare Algebra)|Spektrum]], also die Menge der Eigenwerte der Matrix <math>A</math> gemeint. Daraus folgt, dass das CG-Verfahren ein System zu einer Matrix mit nur <math>k</math> verschiedenen Eigenwerten in <math>k</math> 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 die Systemmatrix ''A'' unsymmetrisch, aber [[Reguläre Matrix|regulär]], so kann das CG-Verfahren auf die [[Normalgleichungen]]
:<math>A^TAx = A^Tb</math>
angewendet werden, da <math>A^TA</math> für eine reguläre Matrix ''A'' symmetrisch und positiv definit ist. Dieses Verfahren nennt sich auch CGNR (von engl. '''C'''onjugate '''G'''radients '''N'''ormal '''R'''esidual), da bei diesem Vorgehen die Norm des Residuums von <math>b-Ax</math> minimiert wird. Alternativ gibt es das Verfahren CGNE (von engl. '''C'''onjugate '''G'''radient Method on the '''N'''ormal '''E'''quations), welches
:<math>AA^Ty=b</math>
löst mit <math>x=A^Ty</math>. Hierbei wird der Fehler minimiert.

CGNR ist ein [[Krylow-Unterraum-Verfahren]] bezüglich <math>A^T A</math>, nicht mehr bezüglich der Systemmatrix <math>A</math>. Bei CGNE ist das ähnlich.
Beide Verfahren haben den Nachteil, dass zum einen <math>A^T</math> zur Verfügung stehen muss, was nicht immer gegeben ist. Zum anderen wird bei den Normalgleichungen die Konditionszahl quadriert, was die Konvergenz typischerweise verlangsamt.

=== Krylow-Unterraum-Verfahren für unsymmetrische Matrizen ===
Direkte Krylow-Unterraum-Verfahren für unsymmetrische Matrizen sind unter anderem [[BiCG-Verfahren|BiCG]] (von engl. '''bic'''onjugate '''g'''radients), [[CGS-Verfahren|CGS]] oder BiCGSTAB. BiCG verwendet gekoppelte Krylowräume zu <math>A</math> und <math>A^T</math>; CGS ist daraus abgeleitet aber vermeidet den expliziten Zugriff auf <math>A^T</math>, während BiCGSTAB eine stabilisierte Variante dieser Verfahrensfamilie ist.

== Literatur ==
* C. T. Kelley: ''Iterative Methods for Linear and Nonlinear Equations.'' SIAM, ISBN 0-89871-352-8. [https://archive.siam.org/books/textbooks/fr16_book.pdf (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: [https://www.cs.cmu.edu/~quake-papers/painless-conjugate-gradient.pdf ''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 ==
<references />

{{SORTIERUNG:Cgverfahren}}
[[Kategorie:Numerische lineare Algebra]]
[[Kategorie:Optimierungsalgorithmus]]

CG-Verfahren - Versionsgeschichte

imported>ConGreif: /* Erweiterung auf unsymmetrische Matrizen */