GMRES-Verfahren

Das GMRES-Verfahren (für Generalized minimal residual method) ist ein iteratives numerisches Verfahren zur Lösung großer, dünnbesetzter linearer Gleichungssysteme. Das Verfahren ist aus der Klasse der Krylow-Unterraum-Verfahren und insbesondere auch für nicht-symmetrische Matrizen geeignet. In exakter Arithmetik, also wenn ohne Rundungsfehler gerechnet wird, liefert das Verfahren nach endlich vielen Schritten die exakte Lösung. Interessanter ist es jedoch als näherungsweises Verfahren, da es mit einer geeigneten Vorkonditionierung auch Gleichungssysteme mit Millionen Unbekannten in wenigen Iterationen mit befriedigender Genauigkeit lösen kann. Damit stellt es eine Art Black-Box-Löser für dünnbesetzte lineare Gleichungssysteme dar. Es wurde 1986 von Yousef Saad und Martin H. Schultz entwickelt.

Das Verfahren

Gegeben sei das lineare Gleichungssystem $Ax=b$ mit einer reellen $n\times n$ Matrix A. Das Gleichungssystem sei eindeutig lösbar, A habe also vollen Rang. Gegeben sei außerdem eine Startnäherung $x_{0}$, etwa einfach die rechte Seite b. Dann wird das GMRES-Verfahren dadurch definiert, dass im m-ten Schritt die euklidische Norm des Residuums $\|Ax-b\|_{2}$ über den affinen Krylow-Unterraum $x_{0}+{\mathcal {K}}_{m}(A,r_{0})=x_{0}+{\mbox{span}}\{r_{0},Ar_{0},\ldots ,A^{m-1}r_{0}\}$ minimiert wird, mit dem Fehler $r_{0}=b-Ax_{0}$ der Startnäherung.

Hierzu wird eine orthonormale Basis $\{v_{1},\dots ,v_{m}\}$ des Raumes mit Hilfe der Arnoldi-Prozedur iterativ berechnet. Diese erlaubt eine Darstellung der von den Basisvektoren gebildeten Matrizen $V_{m}\in \mathbb {R} ^{n\times m}$ und $V_{m+1}\in \mathbb {R} ^{n\times (m+1)}$ über eine Matrix $H_{m}\in \mathbb {R} ^{(m+1)\times m}$, die eine obere Hessenbergmatrix ist, an die eine Zeile angehängt wurde, in der nur der letzte Eintrag nicht Null ist, als

$AV_{m}=V_{m+1}H_{m}$.

Mit dem Ansatz $x_{m}=x_{0}+V_{m}y$ ergibt sich eine effizient berechenbare Form der Norm des Residuums, da $V_{m+1}$ die Euklidische Norm erhält:

$\|Ax_{m}-b\|_{2}=\|\|r_{0}\|_{2}v_{1}-V_{m+1}H_{m}y\|_{2}=\|\|r_{0}\|_{2}e_{1}-H_{m}y\|_{2}$.

Hierbei bezeichnet $e_{1}\in \mathbb {R} ^{m+1}$ den ersten Einheitsvektor. Die Hessenbergmatrix H wird in jedem Schritt aufdatiert und dann durch eine zusammengesetzte orthogonale Transformation $Q_{m}$, meist durch Givens-Rotationen wie im unten angegebenen Pseudo-Code, auf eine rechte obere Dreiecksmatrix ${\bar {R}}_{m}\in \mathbb {R} ^{(m+1)\times m}$ mit Nullen in der letzten Zeile, gebracht. Hier sind nur m-1 Rotationen notwendig, da jede ein Element auf der unteren Nebendiagonalen auf Null setzen kann. In manchen Fällen verlieren die berechneten Vektoren aufgrund von Rundungsfehlern ihre Orthogonalität. Dies kann meist durch Verwendung der aufwändigeren Householder-Spiegelungen statt der Drehungen behoben werden. Anwendung von $Q_{m}$ liefert in beiden Fällen

$\|\|r_{0}\|_{2}e_{1}-H_{m}y\|_{2}=\|Q_{m}(\|r_{0}\|_{2}e_{1}-H_{m}y)\|_{2}=\|{\bar {g}}_{m}-{\bar {R}}_{m}y\|_{2}={\sqrt {|\gamma _{m+1}|^{2}+\|g_{m}-R_{m}y\|_{2}^{2}}}$,

wobei $g_{m}$ und $R_{m}$ aus ihren Pendants durch Weglassen der letzten Zeile erhalten werden. Hier ist nun ersichtlich, an welcher Stelle das Residuum minimal wird, nämlich für den eindeutig bestimmten Vektor y, der $g_{m}=R_{m}y$ erfüllt. Das Residuum im m-ten Schritt ist damit genau $|\gamma _{m+1}|$.

Eine Besonderheit des Verfahrens ist, dass die aktuelle Näherung $x_{m}$ im Laufe der Iteration zunächst nicht berechnet wird, sondern nur der Hilfsvektor y. Stattdessen liefert das Verfahren in jedem Schritt die Norm des Residuums. Ist diese kleiner als die gewünschte Genauigkeit wird das Verfahren üblicherweise abgebrochen. Dann wird die aktuelle Näherung als Linearkombination der Basisvektoren berechnet. Hierbei sind die Komponenten von y einfach die Koeffizienten der Basisdarstellung.

Alternativ ist die Lösung des obigen Minimierungsproblems gegeben als der Vektor $x_{m}$ des affinen Krylow-Unterraumes $x_{0}+{\mathcal {K}}_{m}(A,r_{0})$, dessen Residuum $b-Ax_{m}$ senkrecht auf dem Raum $A{\mathcal {K}}_{m}(A,r_{0})$ steht. Damit ist GMRES eine schiefe Projektionsmethode.

Pseudocode

Gegeben $x_{0}\in \mathbb {R} ^{n}$ und eine Abbruchtoleranz TOL für die Norm des Residuums, berechne $r_{0}=b-Ax_{0}$.

If $\|r_{0}\|_{2}\leq TOL$, then END.

$v_{1}={\frac {r_{0}}{\|r_{0}\|_{2}}}$.

$\gamma _{1}=\|r_{0}\|_{2}$.

For $j=1,\dots ,n$

$q=Av_{j}$

For $i=1,\dots ,j$ do $h_{ij}=v_{i}^{T}q$.

$w_{j}=q-\sum _{i=1}^{j}h_{ij}v_{i},\quad h_{j+1,j}=\|w_{j}\|_{2}.$

For $i=1,\dots ,j-1$ do ${\begin{pmatrix}h_{ij}\\h_{i+1,j}\end{pmatrix}}={\begin{pmatrix}c_{i+1}&s_{i+1}\\-s_{i+1}&c_{i+1}\end{pmatrix}}{\begin{pmatrix}h_{ij}\\h_{i+1,j}\end{pmatrix}}$.

$\beta ={\sqrt {h_{jj}^{2}+h_{j+1,j}^{2}}};\quad s_{j+1}={\frac {h_{j+1,j}}{\beta }}$.

$c_{j+1}={\frac {h_{jj}}{\beta }};\quad h_{jj}=\beta$.

$\gamma _{j+1}=-s_{j+1}\gamma _{j};\quad \gamma _{j}=c_{j+1}\gamma _{j}$.

if $|\gamma _{j+1}|\geq TOL$, $v_{j+1}={\frac {w_{j}}{h_{j+1,j}}}.$

else

for $i=j,\dots ,1$ do $y_{i}={\frac {1}{h_{ii}}}\left(\gamma _{i}-\sum _{k=i+1}^{j}h_{ik}y_{k}\right)$.

$x=x_{0}+\sum _{i=1}^{j}y_{i}v_{i}$.

END.

Konvergenzresultate

Aufgrund der Definition des Verfahrens über das Minimierungsproblem fällt die euklidische Norm der Residuen monoton. In exakter Arithmetik ist GMRES sogar ein direktes Lösungsverfahren, welches spätestens nach n Schritten die exakte Lösung liefert. Wird die Dimension des Krylow-Unterraums in jedem Schritt um eins erhöht, ist diese Aussage klar, da dann im letzten Schritt über den kompletten $\mathbb {R} ^{n}$ minimiert wird. Ist dies nicht der Fall, so kommt es vorher zu einem Abbruch des Algorithmus, allerdings mit der exakten Lösung.

Für allgemeine Matrizen ist dies auch das stärkste Ergebnis das möglich ist, denn nach einem Satz von Greenbaum, Pták und Strakoš gibt es zu jeder monoton fallenden Folge eine Matrix, so dass die Folge der durch GMRES erzeugten Residuen der gegebenen Folge entspricht. Insbesondere ist es also möglich, dass die Residuen konstant bleiben und erst im allerletzten Schritt auf Null fallen.

Für spezielle Matrizen gibt es schärfere Konvergenzresultate. Ist die Matrix diagonalisierbar, so existiert eine reguläre Matrix V und eine Diagonalmatrix $\Lambda$ mit $A=V\Lambda V^{-1}$. Dann gilt für jedes Polynom vom Grad k mit $p(0)=1$:

${\frac {\|r_{k}\|_{2}}{\|r_{0}\|_{2}}}\leq \kappa _{2}(V)\max _{z\in \sigma (A)}|p_{k}(z)|.$

Hierbei bezeichnet $\kappa _{2}$ die Konditionszahl der Matrix in euklidischer Norm und $\sigma (A)$ das Spektrum, also die Menge der Eigenwerte. Für eine normale Matrix ist $\kappa _{2}(V)=1$. Aus der Ungleichung folgt insbesondere, dass die Vorkonditionierung die Eigenwerte zu Clustern zusammenführen sollte.

Ist die Matrix positiv definit (nicht notwendigerweise symmetrisch) so gilt:

$\|r_{m}\|_{2}\leq \left(1-{\frac {\lambda _{\text{min}}^{2}({\frac {A^{T}+A}{2}})}{\lambda _{\text{max}}(A^{T}A)}}\right)^{m/2}\|r_{0}\|_{2}$,

wobei $\lambda _{\text{min}}$ und $\lambda _{\text{max}}$ den größten beziehungsweise kleinsten Eigenwert einer Matrix bezeichnen.

Ist die Matrix A nicht nur positiv definit, sondern auch symmetrisch, so gilt sogar:

$\|r_{m}\|_{2}\leq \left({\frac {\kappa _{2}^{2}(A)-1}{\kappa _{2}^{2}(A)}}\right)^{m/2}\|r_{0}\|_{2}$.

All diese Aussagen gelten nur für die Residuen und geben damit keine Auskunft über den tatsächlichen Fehler, also den Abstand der aktuellen Näherung zur exakten Lösung. Zu diesem sind keine Aussagen bekannt.

Aufwand und Restarted GMRES

GMRES benötigt pro Iteration eine Matrix-Vektor-Multiplikation und eine Reihe von Skalarprodukten, deren Anzahl um eine pro Iterationsschritt steigt, ebenso wie die Anzahl der (vollbesetzten!) zu speichernden Basisvektoren. Dies liegt daran, dass das Verfahren nicht durch eine kurze Rekursion gegeben ist, sondern auf alle Basisvektoren in jedem Schritt zugegriffen wird.

Da der Aufwand und der Speicherplatz linear mit der Iterationszahl steigen, ist es üblich, nach k Schritten die berechnete Basis zu verwerfen und die Iteration mit der aktuellen Näherungslösung neu zu starten (=Restart). Dieses Verfahren wird GMRES(k) genannt, übliche Restart-Längen sind 20 bis 40. Hier lässt sich allerdings nur noch für Spezialfälle Konvergenz beweisen, und es lassen sich Matrizen angeben, so dass ein Restart nicht zu Konvergenz führt.

Der Gesamtaufwand von GMRES ist wie bei allen Krylow-Unterraum-Verfahren bei dünnbesetzten Matrizen O(n) mit einer hohen Konstanten, wenn deutlich weniger Iterationen durchgeführt werden, als es Unbekannte gibt.

Vergleich mit anderen Lösern

Für symmetrische Matrizen fällt das Arnoldi-Verfahren zur Berechnung der orthogonalen Basis mit dem Lanczos-Verfahren zusammen. Das entsprechende Krylow-Unterraum-Verfahren ist das MINRES-Verfahren (für Minimal Residual Method) von Paige und Saunders. Dieses kommt im Gegensatz zur verallgemeinerten Variante mit einer Dreitermrekursion aus. Es lässt sich zeigen, dass es für allgemeine Matrizen kein Krylow-Unterraum-Verfahren gibt, welches mit kurzen Rekursionen arbeitet, aber gleichzeitig wie das GMRES-Verfahren eine Optimalitätsbedingung bezüglich der Norm des Residuums erfüllt.

Eine andere Klasse von Verfahren baut auf dem unsymmetrischen Lanczos-Verfahren auf, insbesondere das BiCG-Verfahren. Solche Verfahren zeichnen sich durch eine Dreitermrekursion aus, allerdings haben sie aufgrund der fehlenden Optimalität keine monotone Konvergenzhistorie mehr. Darüber hinaus liefern sie zwar im Konvergenzfall die exakte Lösung, haben allerdings keine garantierte Konvergenz mehr.

Die dritte Variante sind Verfahren wie CGS und BiCGSTAB. Diese arbeiten ebenfalls mit Dreitermrekursionen ohne Optimalität und können ebenfalls vorzeitig ohne Konvergenz abbrechen. Die Idee bei diesen Verfahren ist es, die generierenden Polynome der Iterationssequenz geschickt zu wählen.

Keine der drei Gruppen ist für alle Matrizen besser, es gibt jeweils Beispiele, wo eine Klasse die anderen übertrumpft. In der Praxis werden deswegen mehrere Löser ausprobiert, um für das gegebene Problem Erfahrungswerte zu sammeln.

Für den Fall symmetrischer, positiv-definiter Matrizen ist das von NVIDIA bereitgestellte CHOLMOD eine Hochleistungsbibliothek für die spärliche Cholesky-Faktorisierung. Die Berechnung kann auf einfache Weise auch auf die GPU ausgelagert werden, wodurch eine wesentliche zusätzliche Beschleunigung der Berechnung erzielt werden kann. CHOLMOD ist Teil des SuiteSparse-Pakets für lineare Algebra, das von Prof. Tim Davis von der Texas A&M University entwickelt wurde. SuiteSparse und CHOLMOD werden seit langem in der Industrie und im akademischen Bereich eingesetzt – insbesondere als Löser linearer Systeme, die mit dem Matlab-Operator „Backslash“ aufgerufen werden (wie in x = A\b).

Vorkonditionierung

Weniger entscheidend als die Auswahl des tatsächlichen Lösers ist die Wahl des Vorkonditionierers, durch den entscheidende Geschwindigkeitsverbesserungen erzielt werden können. Für sequentielle Codes bietet sich hier eine ILU-Zerlegung an, aber je nach Problem können sich auch andere Vorkonditionierer als sehr effizient erweisen. Da ILU nicht parallelisierbar ist, werden in diesem Falle andere eingesetzt, beispielsweise Schwarz-Gebietszerlegungs-Verfahren.

Beispiel Code

Reguläres GMRES (Julia)

using LinearAlgebra

function my_gmres(A::Matrix, b::Vector, x0::Vector, num_iter::Int, tol::Float64)
    n = length(b)
    num_iter = min(num_iter, n)
    
    if x0 == zeros(n)
        r_0 = b
    else
        r_0 = b - A * x0
    end

    b_norm = norm(b)
    r_norm = norm(r_0)

    # initialisiere Matrix V und H
    H = zeros(num_iter+1, num_iter)
    V = zeros(n, num_iter+1)
    V[:, 1] = r_0/r_norm

    # initialisiere Givens-Drehvektoren
    cn = zeros(num_iter)
    sn = zeros(num_iter)
    
    
    β = zeros(num_iter+1)
    β[1] = r_norm

    error = r_norm / b_norm
    e = [error]

    last_iter = 0

    # beginne Iteration
    for k in 1:num_iter

        # beginne Arnoldi-Verfahren
        q_k = A * V[:, k]
        for i in 1:k
            H[i, k] = dot(V[:, i], q_k)
            q_k = q_k - H[i, k] * V[:, i]
        end
        H[k+1, k] = norm(q_k)
        V[:, k+1] = q_k / H[k+1, k]
        # beende Arnoldi-Verfahren

        # Givens-Drehung anwenden
        for i in 1:k-1
            temp     =  cn[i]*H[i,k] + sn[i]*H[i+1,k]
            H[i+1,k] = -sn[i]*H[i,k] + cn[i]*H[i+1,k]
            H[i,k]   =  temp
        end

        # berechne die neuen cos- und sin-Werte für die neue Givens-Drehung
        # nachfolgend eine numerisch stabilere Version des auskommentierten Codes
        #t = sqrt(H[k,k]^2 + H[k+1, k]^2)
        #c[k] = H[k,k]/t
        #s[k] = H[k+1, k]/t

        if H[k+1, k] <= abs(H[k,k])

            t = H[k+1, k]/abs(H[k,k])
            cn[k] = sign(H[k,k])/sqrt(1 + t^2)
            sn[k] = t/sqrt(1 + t^2)
        else

            t = H[k,k]/H[k+1, k]
            cn[k] = t/sqrt(1 + t^2)
            sn[k] = 1/sqrt(1 + t^2)
        end    

        # Givens Drehung nochmals mit neuen cos- und sin-Werten anwenden
        H[k,k] = cn[k]*H[k,k] + sn[k]*H[k+1,k]
        # letzte Givens Drehung H[k+1,k] = -sn[i]*H[i,k] + cn[i]*H[i+1,k] , sollte ohnehin null sein
        H[k+1,k] = 0.0

        # Ende der Givens-Drehung Anwendung

        # Residuum aktualisieren
        β[k+1] = -sn[k]*β[k]
        β[k] = cn[k]*β[k]
        error = abs(β[k+1]) / b_norm

        e = [e; error]

        if error <= tol
            last_iter = k
            break
        end

    end
    # löse das Problem der kleinsten Quadrate
    # weil H eine obere Dreiecksmatrix ist, können wir Rückwärtssubstitution anwenden
    # (H̃ ist eine obere Hessenbergmatrix, wurde aber Givens-gedreht nach H)

    # löse y durch Rückwärtssubstitution
    y = H[1:last_iter, 1:last_iter] \ β[1:last_iter]

    x = x0 + V[:, 1:last_iter] * y

    return x
end

# teste die Funktion
dim = 10
A = I + 1e-1*rand(dim, dim) # sicherstellen, dass A gut konditioniert ist
b = rand(dim)
x0 = zeros(dim)
num_iter = 10
tol = 1e-5

x = my_gmres(A, b, x0, num_iter, tol)
println(norm(A*x - b))

Literatur

• C. T. Kelley: Iterative Methods for Linear and Nonlinear Equations. Society for Industrial and Applied Mathematics SIAM, Philadelphia PA 1995, ISBN 0-89871-352-8, (Frontiers in applied mathematics 16).
• Andreas Meister: Numerik linearer Gleichungssysteme. Eine Einführung in moderne Verfahren. 2. überarbeitete Auflage. Vieweg, Wiesbaden 2005, ISBN 3-528-13135-7.
• Yousef Saad: Iterative Methods for Sparse Linear Systems. 2nd edition. Society for Industrial and Applied Mathematics SIAM, Philadelphia PA 2003, ISBN 0-89871-534-2.
• Yousef Saad, Martin H. Schultz: GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems. In: SIAM Journal on Scientific and Statistical Computing Bd. 7, 1986, ISSN 0196-5204, S. 856–869.