CG-Verfahren

aus Wikipedia, der freien Enzyklopädie
Wechseln zu: Navigation, Suche
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. conjugate gradients oder auch Verfahren der konjugierten Gradienten) ist eine effiziente numerische Methode zur Lösung von großen linearen Gleichungssystemen der Form Ax=b mit symmetrischer, positiv definiter Systemmatrix A.

Das Verfahren liefert, in exakter Arithmetik, nach spätestens m Schritten die exakte Lösung, wobei m die Größe der quadratischen Matrix A\in\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 wurde zuerst 1952 von Eduard Stiefel und Magnus Hestenes vorgeschlagen.[1] Ein für bestimmte Gleichungssysteme äquivalentes Verfahren schlug auch Cornelius Lanczos Anfang der 1950er Jahre mit dem Lanczos-Verfahren vor.

Idee des CG-Verfahrens[Bearbeiten]

Die Idee des CG-Verfahrens besteht darin, dass das Minimieren der quadratischen Form

E(x):=\frac12\langle Ax,x\rangle - \langle b,x\rangle

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

Der Gradient von E an der Stelle x_k ist gerade \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 r_k wie beim Gradientenverfahren in eine andere Richtung d_k die Funktion E über einen Unterraum zu minimieren. Die Richtungen d_k sind dabei alle A-konjugiert, das heißt, es gilt

\langle Ad_i,d_j\rangle=0\qquad\forall i\neq j.

Die Iterierten x_k des CG-Verfahrens werden dann so gewählt, dass sie das Minimum von E in dem affinen Raum V_k, der durch die Vektoren d_0,\ldots,d_k aufgespannt und um x_0 verschoben wird, bilden:

V_k:=x_0+\operatorname{span}\{d_0,\ldots,d_{k-1}\}.

Es lässt sich zeigen, dass ebenfalls gilt:

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 r_0 aufspannen. Das CG-Verfahren lässt sich deswegen alternativ direkt als Krylow-Unterraum-Verfahren definieren.

Da die Vektoren d_k alle A-konjugiert sind, ist die Dimension von V_k gerade k. Ist also A eine m\times m-Matrix, so terminiert das Verfahren nach spätestens m Schritten, falls exakt gerechnet wird. Numerische Fehler können durch weitere Iterationen eliminiert werden. Hierzu betrachtet man den Gradienten r_k, der das Residuum angibt. Unterschreitet die Norm dieses Residuums einen gewissen Schwellenwert, wird das Verfahren abgebrochen.

Das Verfahren baut sukzessive eine A-orthogonale Basis für den \mathbb R^m 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[Bearbeiten]

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

r_0 = b - A x_0
d_0 = r_0

Für k = 0,1,... führt man aus:

\begin{align}
z&=Ad_k
\end{align}
  • Finde von x_k in Richtung d_k den Ort x_{k+1} des Minimums der Funktion E und aktualisiere den Gradienten bzw. das Residuum
\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}
  • Korrigiere die Suchrichtung d_{k+1} mit Hilfe von d_k und r_{k+1}
\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}

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

Varianten[Bearbeiten]

Es existieren verschiedene Varianten des Verfahrens, neben der ersten von Roger Fletcher und Colin Reeves z. B. von Hestenes und Stiefel, von Davidon, Fletcher und Powell oder von Polak und 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.

  • \beta_{k} = \frac{(r_{k+1})^T r_{k+1}}{(r_k)^T r_k} (Fletcher-Reeves)
  • \beta_{k} = \frac{(r_{k+1})^T (r_{k+1}-r_k)}{(r_k)^T r_k} (Polak-Ribière)
  • \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)[Bearbeiten]

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 Ax=b mit Hilfe einer Vorkonditionierer-Matrix C=KK^T\approx A^{-1} zu K^TAKy=K^Tb mit y=K^{-1}x transformiert, und darauf das CG-Verfahren angewandt.

Die Matrix K^TAK ist symmetrisch, da A symmetrisch ist. Sie ist ferner positiv definit, da nach dem Trägheitssatz von Sylvester A und K^TAK die gleichen Anzahlen positiver und negativer Eigenwerte besitzen.

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

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

r_0 = b - A x_0
h_0 = C r_0
d_0 = h_0

Für k = 0,1,... setzt man:

  • Speichere Matrix-Vektor-Produkt, um es nur einmal auszurechnen
z=Ad_k
  • Finde von x_k in Richtung d_k das Minimum x_{k+1} und aktualisiere Gradienten und vorkonditionierten Gradienten
\alpha_k=\frac{\langle r_k, h_k\rangle}{\langle d_k, z\rangle}
x_{k+1}=x_k+\alpha_k d_k
r_{k+1}=r_k-\alpha_k z (Residuum)
h_{k+1}=C r_{k+1}
  • Korrigiere die Suchrichtung d_{k+1}
\beta_k=\frac{\langle r_{k+1}, h_{k+1}\rangle}{\langle r_k, h_k\rangle}
d_{k+1}=h_{k+1}+\beta_k d_k

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

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 van der Vorst eingeführt.

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


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 \omega \in (0, \,2), wobei D die Hauptdiagonale und L die strikte untere Dreiecksmatrix von A ist.

Konvergenzrate des CG-Verfahrens[Bearbeiten]

Man kann zeigen, dass die Konvergenzgeschwindigkeit des CG-Algorithmus durch

\|x_k-x\|_A \le 2\left(\frac{\sqrt{\kappa(A)}-1}{\sqrt{\kappa(A)}+1}\right)^k\|x_{0}-x\|_A

beschrieben wird. Hierbei ist \kappa(A) die Kondition der Matrix A bzgl. der Spektralnorm, also der von der euklidischen Norm erzeugten Matrixnorm, sowie \|x\|_A = \sqrt{x^T A x} die Energienorm von A. Der Ausdruck \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 A symmetrisch und positiv definit ist, gilt

\kappa(A) = \frac{\lambda_{max}(A)}{\lambda_{min}(A)}.

Aus der Minimierungseigenschaft lässt sich ferner herleiten, dass

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

wobei p_k(z) ein beliebiges Polynom vom Grad k ist mit p_k(0)=1 und x^* die Lösung. Mit \sigma(A) ist das Spektrum, also die Menge der Eigenwerte der Matrix A gemeint. Daraus folgt, dass CG ein System zu einer Matrix mit nur k verschiedenen Eigenwerten in k Schritten löst und dass CG 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[Bearbeiten]

Ist die Systemmatrix A unsymmetrisch, aber regulär, so kann das CG-Verfahren auf die Normalgleichungen

A^TAx = A^Tb

angewendet werden, da A^TA für eine reguläre Matrix A symmetrisch und positiv definit ist. Dieses Verfahren nennt sich auch CGNR, da bei diesem Vorgehen die Norm des Residuums von b-Ax minimiert wird. Alternativ gibt es das Verfahren CGNE, welches

AA^Ty=b

löst mit x=A^Ty. Hierbei wird der Fehler minimiert.

Beide Verfahren haben den Nachteil, dass zum einen 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[Bearbeiten]

Einzelnachweise[Bearbeiten]

  1. Hestenes, Stiefel: Methods of conjugate gradients for solving linear systems, Journal of Research of the National Bureau of Standards, Bd. 49, 1952, S.409–436, online (pdf)