MINRES-Verfahren

Ein Vergleich der Norm des Fehlers sowie des Residuums beim CG-Verfahren (blau) und dem MINRES-Verfahren (grün). Die verwendete Matrix stammt aus einem 2d-Randwertproblem.

Das MINRES-Verfahren (von englisch minimum residual „minimales Residuum“) ist ein Krylow-Unterraum-Verfahren zur iterativen Lösung von symmetrischen linearen Gleichungssystemen. Es wurde 1975 von den Mathematikern Christopher Conway Paige und Michael Alan Saunders vorgeschlagen.

Im Unterschied zum CG-Verfahren wird beim MINRES-Verfahren nicht vorausgesetzt, dass die Matrix positiv definit ist, nur die Symmetrie der Matrix wird zwingend vorausgesetzt.

Eigenschaften des MINRES-Verfahrens

Das MINRES-Verfahren berechnet iterativ eine Näherungslösung eines linearen Gleichungssystems der Form

$Ax=b$.

Dabei ist $A\in \mathbb {R} ^{n\times n}$ eine symmetrische Matrix und $b\in \mathbb {R} ^{n}$ die rechte Seite.

Hierzu wird die Norm des Residuums $r(x):=b-Ax$ im $k$-dimensionalen Krylowraum

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

minimiert. Dabei ist $x_{0}\in \mathbb {R} ^{n}$ ein Startwert und $r_{0}:=r(x_{0})$.

Genauer definieren wir die Näherungslösungen $x_{k}$ durch

$x_{k}:=\mathrm {argmin} _{x\in V_{k}}\|r(x)\|$.

Dabei ist $\|\cdot \|$ die euklidische Norm im $\mathbb {R} ^{n}$.

Wegen der Symmetrie von $A$ ist es dabei im Gegensatz zum GMRES-Verfahren möglich, diesen Minimierungsprozess rekursiv durchzuführen. Im $k$-ten Schritt ist jeweils nur der Rückgriff auf die Iterierten aus den letzten beiden Schritten nötig (kurze Rekursion).

Konstruktionsidee von MINRES

Krylov-Unterraum-Verfahren beschränken sich darauf, nur Vektoren aus Matrix-Vektor-Produkten mit der Systemmatrix zu benutzen. Das hat Vorteile, weil die Matrix dazu nicht explizit, sondern nur als Funktion für das Matrix-Vektor-Produkt verfügbar sein muss. Zu Beginn sind die einzigen bekannten Vektoren die aktuelle Näherungslösung $x_{0}$ (zu Beginn meist ein Nullvektor), die rechte Seite $b$ und das Residuum $r_{0}=b-A\cdot x_{0}$.

Man kopiert das Residuum in einen Vektor $p_{0}$ und nimmt diesen aus obigem Grund als Korrekturrichtung für die Näherungslösung. Dazu berechnet man sein Bild $s_{0}=A\cdot p_{0}$. Dieses Bild will man optimal an das Residuum heranaddieren, sodass dessen Länge kleinstmöglich wird (daher der Verfahren-Name). Dazu rechnet man $r_{1}=r_{0}-\xi \cdot s_{0}$, mit $r_{1}\perp s_{0}$. Dazu muss $\xi =\langle s_{0},r_{0}\rangle /\langle s_{0},s_{0}\rangle$ sein (Gram-Schmidt). Die zugehörige Näherungslösung $x_{1}$ für dieses Residuum kennt man: $x_{1}=x_{0}+\xi \cdot p_{0}$.

Für das neue Residuum erstellt man wieder eine Kopie $p_{1}$ und berechnet wieder das Bild $s_{1}=A\cdot p_{1}$. Um durch Wiederholung dieses Prinzips das Residuum immer weiter zu verkleinern, möchte man im nächsten Schritt ein Residuum $r_{2}$ erzeugen, dass auf $s_{0}$ und $s_{1}$ senkrecht steht. Da $s_{1}$ einen Richtungsanteil von $s_{0}$ enthalten könnte, muss $s_{1}$ auf $s_{0}$ orthogonalisiert und $p_{1}$ analog angepasst werden, damit danach weiterhin $s_{1}=A\cdot p_{1}$ gilt. Für $s_{1}:=s_{1}-\tau \cdot s_{0}$ wird $p_{1}:=p_{1}-\tau \cdot p_{0}$. So setzt man das über viele Iterationen fort.

Auf diese Weise müsste in der $i$-ten Iteration die Richtung $s_{i}$ auf $i-1$ Vorgängern $s_{1},\dotsc ,s_{i-1}$ orthogonalisiert werden. Lanczos konnte jedoch zeigen, dass $s_{i}$ bereits auf all diesen Richtungen senkrecht steht, falls $s_{i}$ nur auf seinen beiden Vorgängern $s_{i-1},s_{i-2}$ orthogonalisiert (= senkrecht gestellt) wird. Dies liegt an der Symmetrie von $A$ (weshalb das Verfahren auch nur im symmetrischen Fall funktioniert).

Algorithmus des MINRES-Verfahrens

Anmerkung: Das MINRES-Verfahren ist vergleichsweise komplizierter als das algebraisch äquivalente Conjugate Residual Verfahren. Im Folgenden wurde daher ersatzweise das Conjugate Residual (CR) Verfahren niedergeschrieben. Es unterscheidet sich insoweit von MINRES, dass in CR nicht wie bei MINRES die Spalten einer Basis des Krylov-Raums (unten bezeichnet mit $p_{k}$), sondern deren Bilder (unten bezeichnet mit $s_{k}$) über die Lanczos-Rekursion orthogonalisiert werden. Es existieren effizientere und präkonditionierte Varianten mit weniger AXPYs. Vgl. dazu mit dem englischsprachigen Artikel.

Zunächst wählt man $x_{0}\in \mathbb {R} ^{n}$ beliebig und berechnet

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

$p_{0}=r_{0}$

$s_{0}=Ap_{0}$

Dann iterieren wir für $k=1,2,\dots$ die folgenden Schritte:

• Berechne die $x_{k},r_{k}$ durch

$\alpha _{k-1}={\frac {\langle r_{k-1},s_{k-1}\rangle }{\langle s_{k-1},s_{k-1}\rangle }}$

$x_{k}=x_{k-1}+\alpha _{k-1}p_{k-1}$

$r_{k}=r_{k-1}-\alpha _{k-1}s_{k-1}$

falls $\|r_{k}\|$ kleiner als eine vorgegebene Toleranz ist, bricht man an dieser Stelle den Algorithmus mit der Näherungslösung $x_{k}$ ab, ansonsten berechnet man eine neue Abstiegsrichtung $p_{k}$ mittels

$p_{k}\leftarrow s_{k-1}$

$s_{k}\leftarrow As_{k-1}$

• für $l=1,2$ (der Schritt $l=2$ wird erst ab dem zweiten Iterationsschritt durchgeführt) berechne:

$\beta _{k,l}={\frac {\langle s_{k},s_{k-l}\rangle }{\langle s_{k-l},s_{k-l}\rangle }}$

$p_{k}\leftarrow p_{k}-\beta _{k,l}p_{k-l}$

$s_{k}\leftarrow s_{k}-\beta _{k,l}s_{k-l}$

Konvergenzrate des MINRES-Verfahrens

Im Fall von positiv definiten Matrizen lässt sich die Konvergenzrate des MINRES-Verfahrens ähnlich wie beim CG-Verfahren abschätzen. Im Gegensatz zum CG-Verfahren gilt die Abschätzung allerdings nicht für die Fehler der Iterierten, sondern für das Residuum. Es gilt:

$\|r_{k}\|\leq 2\left({\frac {{\sqrt {\kappa (A)}}-1}{{\sqrt {\kappa (A)}}+1}}\right)^{k}\|r_{0}\|$.

Dabei ist $\kappa (A)$ die Konditionszahl der Matrix $A$.

Beispiel-Implementierung in GNU Octave / Matlab

function [x,r] = minres(A,b,x0,maxit,tol)
  x = x0;
  r = b - A*x0;
  p0 = r;
  s0 = A*p0;
  p1 = p0;
  s1 = s0;
  for iter=[1:maxit]
    p2 = p1;p1 = p0;
    s2 = s1;s1 = s0;
    alpha = r'*s1/(s1'*s1);
    x += alpha*p1;
    r -= alpha*s1;
    if (r'*r < tol^2)
      break
    end
    p0 = s1;
    s0 = A*s1;
    beta1 = s0'*s1/(s1'*s1);
    p0 -= beta1*p1;
    s0 -= beta1*s1;
    if iter > 1
      beta2 = s0'*s2/(s2'*s2);
      p0 -= beta2*p2;
      s0 -= beta2*s2;
    end
  end
end