Gaußsches Eliminationsverfahren

Grundidee und Ablauf

Das gaußsche Eliminationsverfahren (Gauß-Verfahren) ist ein Algorithmus der linearen Algebra und Numerik zum Lösen linearer Gleichungssysteme. Für ein System Ax=b werden Äquivalenzumformungen verwendet: Sie verändern die Gleichungen, erhalten aber die Lösungsmenge. Ein eindeutig lösbares System lässt sich so in Stufenform überführen und anschließend lösen.

Das Verfahren hat zwei Schritte. Bei der Vorwärtselimination werden unterhalb der Diagonale Einträge zu null gemacht, sodass von Zeile zu Zeile mindestens eine Variable weniger vorkommt. Dazu darf man ein Vielfaches einer Zeile zu einer anderen addieren und Zeilen vertauschen. Ist das aktuelle Diagonalelement null, muss ein Nichtnulleintrag durch Zeilenvertauschung auf diese Stelle gebracht werden. Danach beginnt das Rückwärtseinsetzen: Aus der letzten Gleichung wird die letzte Variable berechnet und schrittweise in die darüberliegenden Gleichungen eingesetzt.

Spaltenvertauschungen sind für den Ablauf nicht nötig, werden aber in Programmen teils aus Stabilitätsgründen eingesetzt; dabei ändert sich die Reihenfolge der Variablen. Eine Multiplikation einer Zeile mit einer Zahl kann beim Rechnen per Hand Brüche vermeiden, verursacht aber zusätzlichen Aufwand und verändert die Determinante. Beim Gauß-Jordan-Algorithmus werden zusätzlich die Einträge oberhalb der Diagonale eliminiert; es entsteht eine Diagonalform und das Rückwärtseinsetzen entfällt.

Beispiel und Kontrolle

Für das Gleichungssystem x₁+2x₂+3x₃=2, x₁+x₂+x₃=2 und 3x₁+3x₂+x₃=0 lautet die erweiterte Koeffizientenmatrix

(1 2 3 | 2; 1 1 1 | 2; 3 3 1 | 0).

Mit der ersten Zeile als Pivotzeile werden zur zweiten Zeile das (−1)-fache und zur dritten Zeile das (−3)-fache der ersten Zeile addiert. Anschließend wird zur neuen dritten Zeile das (−3)-fache der neuen zweiten Zeile addiert. Die letzte Gleichung ist −2x₃=−6, also x₃=3. Aus der zweiten Gleichung −x₂−2x₃=0 folgt x₂=−6; schließlich erhält man x₁=5.

Eine Kontrolle ist über Zeilensummen möglich. Werden Zeilen umgeformt, muss dieselbe Umformung auf ihre Summen angewandt werden. Im Beispiel hat die erste erweiterte Zeile die Summe 1+2+3+2=8. Addiert man zu einer Zeile das (−1)-fache der ersten, wird auch ihre Summe entsprechend verändert; stimmt das Ergebnis mit der Summe der neuen Zeile überein, kontrolliert dies die Rechnung.

Pivotisierung und LR-Zerlegung

Pivotisierung bedeutet, vor einem Eliminationsschritt ein geeignetes Pivotelement ungleich 0 auszuwählen. Ist etwa a₁₁=0, kann ohne Zeilenvertauschung nicht begonnen werden. Beim Rechnen per Hand sind 1 oder −1 als Pivot günstig. Für Computer wird meist ein betragsgrößtes Element gewählt, um die Stabilität zu verbessern. Liegt das Pivot in der aktuellen Spalte, heißt dies Spaltenpivotisierung; alternativ kann es in der aktuellen Zeile gewählt werden. Bei vollständiger Pivotisierung oder Totalpivotisierung wird das betragsgrößte Element der gesamten Restmatrix gewählt, wofür im Allgemeinen Zeilen- und Spaltenvertauschungen nötig sind. Vertauschungen können in einem Indexvektor gespeichert werden.

Die LR-Zerlegung interpretiert das Verfahren für eine reguläre quadratische Matrix als A=L·R. L ist eine untere normierte Dreiecksmatrix, also mit Diagonalelementen 1, und R eine obere Dreiecksmatrix in Stufenform. Im Beispiel gilt

A=(1 2 3; 1 1 1; 3 3 1)=(1 0 0; 1 1 0; 3 3 1)·(1 2 3; 0 −1 −2; 0 0 −2)=L·R.

L speichert die Umformungsschritte. Mit nötigen Zeilenvertauschungen lautet die Zerlegung P·A=L·R, wobei P eine Permutationsmatrix ist: Sie entsteht aus der Einheitsmatrix durch Zeilenvertauschungen und enthält weiterhin nur Nullen und Einsen. Für jede reguläre Matrix A∈ℝⁿ×ⁿ existieren P, L und R dieser Art.

Lösen mit Dreiecksmatrizen

Aus Ax=b und PA=LR folgt LRx=Pb. Mit y:=Rx und b̂:=Pb entstehen zwei leicht lösbare Systeme: Ly=b̂ und Rx=y. Für mehrere rechte Seiten b ist dies besonders effizient, weil die LR-Zerlegung nur einmal berechnet werden muss.

Beim Vorwärtseinsetzen wird y von oben nach unten bestimmt. Für Ly=b gilt

yᵢ=(1/lᵢᵢ)(bᵢ−∑ₖ₌₁ⁱ⁻¹ lᵢₖ·yₖ),

beginnend mit y₁=b₁/l₁₁. Beim Rückwärtseinsetzen wird x von unten nach oben bestimmt, beginnend mit xₙ=yₙ/rₙₙ:

xᵢ=(1/rᵢᵢ)(yᵢ−∑ₖ₌ᵢ⁺¹ⁿ rᵢₖ·xₖ).

Bei dünnbesetzten Matrizen wird eine LR-Zerlegung häufig vollbesetzt. Eine unvollständige LU-Zerlegung berechnet deshalb nur Einträge eines vorgegebenen Besetzungsmusters. Sie approximiert A und kann als Vorkonditionierer für iterative Verfahren dienen; für symmetrische positiv definite Matrizen heißt der entsprechende Fall unvollständige Cholesky-Zerlegung.

Aufwand, Genauigkeit und weitere Anwendungen

Für eine n×n-Matrix benötigt die LR-Zerlegung circa 2/3 n³ arithmetische Operationen. Vorwärts- und Rückwärtseinsetzen haben Aufwand O(n²). Das Gauß-Verfahren ist damit ein schnelles direktes Verfahren; eine QR-Zerlegung benötigt mindestens doppelt so viele Operationen. Es wird für kleine bis mittlere Dimensionen, bis etwa n=10000, empfohlen. Bei höheren Dimensionen sind iterative Verfahren oft besser; bei vollbesetzten Matrizen benötigen sie je Schritt O(n²), ihre Konvergenzgeschwindigkeit hängt jedoch stark von der Matrix ab.

Die Rechnung kann direkt im Speicher von A erfolgen. Eine vollbesetzte Matrix mit n=1000 hat eine Million Koeffizienten und benötigt im IEEE-754-Format double etwa 8 Megabyte. Für symmetrische positiv definite Matrizen halbiert die Cholesky-Zerlegung Rechenaufwand und Speicher. Bei Bandmatrizen fester Bandbreite m bleibt die Bandstruktur erhalten; der Aufwand sinkt auf O(nm²).

Ohne Pivotisierung ist das Verfahren im Allgemeinen instabil. Mit Spaltenpivotisierung ist es für die meisten Matrizen stabil; vollständige Pivotisierung verbessert die Stabilität weiter, macht die Pivotsuche aber zu O(n³). QR-Zerlegungen sind generell stabiler, jedoch aufwändiger. Bei strikt diagonaldominanten oder positiv definiten Matrizen ist das Gauß-Verfahren auch ohne Pivotisierung stabil. Die Genauigkeit hängt außerdem von der Kondition der Matrix und der Maschinengenauigkeit ab.

Eine Nachiteration verbessert eine berechnete Lösung x₀=x. Man berechnet rₖ=b−Axₖ, löst Azₖ=rₖ mithilfe der LR-Zerlegung und setzt xₖ₊₁=xₖ+zₖ. Meist reichen wenige Schritte; für rₖ ist oft höhere Genauigkeit nötig. Die LAPACK-Routine DSGESV bestimmt etwa die LR-Zerlegung in einfacher Genauigkeit und erreicht durch Nachiteration mit doppeltgenau berechnetem Residuum doppelte Genauigkeit.

Theoretisch zeigt das Verfahren auch die Lösbarkeit: Nichtnullwerte der rechten Seite zu Nullzeilen der reduzierten Stufenform bedeuten Unlösbarkeit. Andernfalls ist das System lösbar; die Zahl freier Parameter ist Zahl der Unbekannten minus Rang. Die Determinante der oberen Dreiecksmatrix ist das Produkt ihrer Diagonalelemente; Zeilenvertauschungen ändern nur ihr Vorzeichen. Zur Inversenberechnung wird A rechts um die Einheitsmatrix erweitert und bis links die Einheitsmatrix steht; rechts steht dann A⁻¹. Dies ist numerisch nicht zu empfehlen.