Gaußsches Eliminationsverfahren

Das gaußsche Eliminationsverfahren oder einfach Gauß-Verfahren (nach Carl Friedrich Gauß) ist ein Algorithmus aus den mathematischen Teilgebieten der linearen Algebra und der Numerik. Es ist ein wichtiges Verfahren zum Lösen von linearen Gleichungssystemen und beruht darauf, dass elementare Umformungen zwar das Gleichungssystem ändern, aber die Lösung erhalten. Dies erlaubt es, jedes eindeutig lösbare Gleichungssystem auf Stufenform zu bringen, an der die Lösung durch sukzessive Elimination der Unbekannten leicht ermittelt oder die Lösungsmenge abgelesen werden kann.

Die Anzahl der benötigten Operationen ist bei einer $n\times n$ -Matrix von der Größenordnung n^3 . In seiner Grundform ist der Algorithmus aus numerischer Sicht anfällig für Rundungsfehler, aber mit kleinen Modifikationen (Pivotisierung) stellt er für allgemeine lineare Gleichungssysteme das Standardlösungsverfahren dar und ist Teil aller wesentlichen Programmbibliotheken für numerische lineare Algebra wie NAG, IMSL und LAPACK.

Erklärung

Ein lineares Gleichungssystem Ax=b mit drei Gleichungen und drei Unbekannten x = (x_1,~x_2,~x_3)^T und rechter Seite b = (b_1,~b_2,~b_3)^T hat die Form:

$\begin{array}{rcrcrcl} a_{11}x_1 & + & a_{12}x_2 & + & a_{13}x_3 & = & b_1, \\ a_{21}x_1 & + & a_{22}x_2 & + & a_{23}x_3 & = & b_2, \\ a_{31}x_1 & + & a_{32}x_2 & + & a_{33}x_3 & = & b_3. \end{array}$

Der Algorithmus zur Berechnung der Variablen $x_{1}$ , $x_{2}$ und $x_{3}$ lässt sich in zwei Etappen einteilen:

Vorwärtselimination,
Rückwärtseinsetzen (Rücksubstitution).

Im ersten Schritt wird das Gleichungssystem auf Stufenform gebracht. Stufenform heißt, dass pro Zeile mindestens eine Variable weniger auftritt, also mindestens eine Variable eliminiert wird. Im obigen Gleichungssystem würde man $a_{21}$ , $a_{31}$ und $a_{32}$ eliminieren, in der dritten Zeile ist dann nur noch die Variable $x_{3}$ :

$\begin{array}{rcrcrcl} \tilde{a}_{11}x_1 &+& \tilde{a}_{12}x_2 &+& \tilde{a}_{13}x_3 &=& \tilde{b}_1 \\ & & \tilde{a}_{22}x_2 &+& \tilde{a}_{23}x_3 &=& \tilde{b}_2 \\ & & & & \tilde{a}_{33}x_3 &=& \tilde{b}_3. \end{array}$

Zum Erreichen der Stufenform werden elementare Zeilenumformungen benutzt, mit Hilfe derer das Gleichungssystem in ein neues transformiert wird, welches aber dieselbe Lösungsmenge besitzt. Ausreichend sind zwei Arten von elementaren Zeilenumformungen:

Eine Zeile oder das Vielfache einer Zeile zu einer anderen Zeile addieren.
Zwei Zeilen vertauschen.

Das Verfahren besteht dann darin, angefangen in der ersten Spalte mit Umformungen der ersten Art durch geschicktes Dazuaddieren der ersten Zeile alle Einträge bis auf den ersten zu Null zu machen. Dies wird dann in der so modifizierten zweiten Spalte fortgesetzt, wobei diesmal Vielfache der zweiten Zeile zu den folgenden Zeilen addiert werden und so weiter. Dieser Schritt funktioniert nur, wenn das Diagonalelement der aktuellen Spalte nicht Null ist. In so einem Fall ist die zweite Art der Zeilenumformung nötig, da durch eine Zeilenvertauschung ein Nichtnulleintrag auf der Diagonale erzeugt werden kann. Mit Hilfe dieser beiden Arten von Umformungen ist es möglich, jedes lineare Gleichungssystem auf Stufenform zu bringen.

Eine weitere Art der elementaren Umformung ist das Vertauschen von Spalten. Diese wird zur Durchführung des Algorithmus nicht benötigt, aber manchmal in Computerprogrammen aus Stabilitätsgründen eingesetzt. Dabei wird die Position der Variablen im Gleichungssystem geändert. Beim Rechnen per Kopf ist manchmal noch die Multiplikation einer Zeile mit einer Zahl nützlich, etwa um komplizierte Brüche zu vermeiden. Dies verursacht zusätzlichen Rechenaufwand und ist deswegen in Computerprogrammen keine Option und ändert ferner die Determinante der Koeffizientenmatrix, was theoretische Nachteile mit sich bringt.

Im zweiten Schritt des Verfahrens, dem Rückwärtseinsetzen, werden ausgehend von der letzten Zeile, in der nur noch eine Variable auftaucht, die Variablen ausgerechnet und in die darüberliegende Zeile eingesetzt.

Eine Alternative hierzu ist der Gauß-Jordan-Algorithmus, bei dem nicht nur die unteren Teile eliminiert werden, sondern auch die oberen, so dass eine Diagonalform entsteht, bei der dann der oben genannte zweite Schritt entfällt.

Beispiel

, hier: $a_{11} = 1$ , $a_{12} = 2$ , $a_{13} = 3$ und

Zur besseren Übersichtlichkeit werden die Koeffizienten $a_{ij}$ des linearen Gleichungssystems in die mit erweiterte Koeffizientenmatrix geschrieben:

$\left( \begin{array}{c c c | c} 1 & 2 & 3 & 2\\ 1 & 1 & 1 & 2\\ 3 & 3 & 1 & 0 \end{array} \right)$

Jetzt wird so umgeformt, dass $a_{21}$ und $a_{31}$ Null werden, indem man geeignete Vielfache der ersten Gleichung zur zweiten und dritten Gleichung addiert. Den entsprechenden Multiplikator erhält man, indem man das zu eliminierende Element (als erstes $a_{21}=1$ ) durch das Pivotelement $a_{11}$ teilt (hier: $\tfrac{1}{1}=1$ und $\tfrac{3}{1}=3$ ). Da die beiden Elemente $a_{21}$ und $a_{31}$ Null werden sollen, werden die beiden Multiplikatoren jeweils mit (-1) multipliziert.

Zur zweiten Zeile wird also das (-1) -fache und zur dritten Zeile das (-3) -fache der ersten Zeile addiert. Damit $a_{32}$ Null wird, wird ein Vielfaches der zweiten Zeile zur dritten Zeile addiert, in diesem Fall das (-3) -fache:

Falls die Zahl, durch die zur Berechnung des Multiplikators dividiert wird (hier für die ersten beiden Zeilen die Zahl , beim dritten Mal die Zahl ), Null ist, wird diese Zeile mit einer weiter unten liegenden vertauscht. Die letzte Zeile bedeutet

Diese Gleichung ist einfach lösbar und liefert x_3 = 3 . Damit ergibt sich für die zweite Zeile

, mit

also

und weiter x_1 = 5 . Damit sind alle Variablen berechnet:

und

Kontrolle durch Zeilensumme

Die Umformungen können durch das Berechnen der Zeilensumme kontrolliert werden.

Hier wurde in der letzten Spalte die Summe aller Elemente der jeweiligen Zeile addiert. Für die erste Zeile ist die Zeilensumme 1 + 2 + 3 + 2 = 8 . Da an der ersten Zeile keine Umformungen durchgeführt werden, ändert sich ihre Zeilensumme nicht. Bei der ersten Umformung dieses Gleichungssystems wird zur zweiten Zeile das (-1) -fache der ersten addiert. Macht man das auch für die Zeilensumme, so gilt $5 + (-1) \cdot 8 = -3$ . Dieses Ergebnis ist die Zeilensumme der umgeformten zweiten Zeile -1 -2 + 0 = -3 . Zur Überprüfung der Rechnungen kann man also die Umformungen an der Zeilensumme durchführen. Sind alle Rechnungen korrekt, muss sich die Zeilensumme der umgeformten Zeile ergeben.

Pivotisierung

Das gaußsche Eliminationsverfahren ist im Allgemeinen nicht ohne Zeilenvertauschungen durchführbar. Ersetzt man im obigen Beispiel $a_{11}=1$ durch $a_{11}=0$ , so kann der Algorithmus ohne Zeilenvertauschung gar nicht starten. Zur Abhilfe wählt man ein Element der ersten Spalte der Koeffizientenmatrix, das so genannte Pivotelement welches ungleich 0 ist.

Wahl des Pivot
$x_{1}$	$x_{2}$	$x_{3}$
0	2	3	4
1	1	1	2
3	3	1	0

Danach vertauscht man die erste Zeile mit der Pivotzeile:

Zeilenvertauschung
$x_{1}$	$x_{2}$	$x_{3}$
1	1	1	2
0	2	3	4
3	3	1	0

Für die Rechnung per Hand ist es hilfreich, eine 1 oder minus 1 als Pivotelement zu wählen, damit im weiteren Verlauf des Verfahrens keine Brüche entstehen. Für die Berechnung mit Hilfe eines Computers ist es sinnvoll, das betragsgrößte Element zu wählen, um einen möglichst stabilen Algorithmus zu erhalten. Wählt man das Pivotelement in der aktuellen Spalte, spricht man von Spaltenpivotisierung. Alternativ kann man das Pivot auch in der aktuellen Zeile wählen.

Wahl des Pivot
$x_{1}$	$x_{2}$	$x_{3}$
0	2	3	4
1	1	1	2
3	3	1	0

In diesem Fall werden entsprechend die Spalten getauscht.

Spaltenvertauschung
$x_{2}$	$x_{1}$	$x_{3}$
2	0	3	4
1	1	1	2
3	3	1	0

Beim Rückwärtseinsetzen ist dabei zu beachten, dass die Variablen ihre Position im Gleichungssystem geändert haben. Wählt man als Pivot das betragsgrößte Element der gesamten Restmatrix, so spricht man von vollständiger Pivotisierung beziehungsweise Totalpivotisierung. Dafür sind im Allgemeinen sowohl Zeilen- als auch Spaltenvertauschungen notwendig.

Pivotisierung ist ohne nennenswerten Zusatzaufwand durchführbar, wenn nicht die Einträge der Matrix und der rechten Seite vertauscht, sondern die Vertauschungen in einem Indexvektor gespeichert werden.

LR-Zerlegung

Will man das Lösen eines quadratischen eindeutig lösbaren Gleichungssystems Ax=b als Computerprogramm umsetzen, bietet es sich an, den Gaußalgorithmus als LR-Zerlegung (auch LU-Zerlegung oder Dreieckszerlegung genannt) zu interpretieren. Dies ist eine Zerlegung der regulären Matrix in das Produkt einer linken unteren Dreiecksmatrix (links, bzw. engl. „lower“) und einer rechten oberen Dreiecksmatrix (rechts, auch mit bezeichnet, von engl. „upper“). Das folgende Beispiel zeigt dies:

$A= \begin{pmatrix} 1 & 2 & 3 \\ 1 & 1 & 1 \\ 3 & 3 & 1 \\ \end{pmatrix} = \begin{pmatrix} 1 & 0 & 0 \\ 1 & 1 & 0 \\ 3 & 3 & 1 \\ \end{pmatrix} \cdot \begin{pmatrix} 1 & 2 & 3 \\ 0 & -1 & -2 \\ 0 & 0 & -2 \\ \end{pmatrix} =L\cdot R$

Dabei hat die oben erwähnte Stufenform und die Matrix dient dem Speichern der benötigten Umformungsschritte, die Multiplikationen mit Frobeniusmatrizen entsprechen. Das zeigt die Existenz der Zerlegung. Um Eindeutigkeit zu erreichen, werden die Diagonalelemente der Matrix als 1 festgelegt. Die Umformungsschritte zu speichern hat den Vorteil, dass für verschiedene „rechte Seiten“ das Gleichungssystem effizient durch Vorwärts- und Rückwärtseinsetzen gelöst werden kann.

Die im Allgemeinen benötigten Zeilenvertauschungen können durch eine Permutationsmatrix beschrieben werden.

$P\cdot A=L\cdot R$

Vorwärtseinsetzen

Beim Vorwärtseinsetzen berechnet man eine Lösung des linearen Gleichungssystems Ly = b , beziehungsweise bei Rechnung mit Pivotisierung von Ly = Pb . Diese steht über die Gleichung y = Rx mit der Lösung des ursprünglichen Gleichungssystems in Beziehung.

Ausgeschrieben hat das Gleichungssystem Ly = b folgende Gestalt:

$\begin{pmatrix} l_{11} & 0 & 0 & \ldots & 0\\ l_{21} & l_{22} & 0 & & \vdots\\ l_{31} & l_{32} & l_{33} & \ddots & \vdots\\ \vdots & \vdots & \vdots & \ddots & 0\\ l_{n1} & l_{n2} & l_{n3} & \ldots & l_{n,n} \end{pmatrix} \cdot \begin{pmatrix}y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_n\end{pmatrix} = \begin{pmatrix}b_1 \\ b_2 \\ b_3 \\ \vdots \\ b_n\end{pmatrix}$

Für die Komponenten $y_{i}$ gilt dann die folgende Formel:

$y_i = \frac{1}{l_{ii}} \left( b_i - \sum_{k=1}^{i-1} l_{ik} \cdot y_k \right).$

Beginnend mit $y_1 = \frac{b_1}{l_{11}}$ können nacheinander $y_1, y_2, \ldots,y_n$ ausgerechnet werden, indem jeweils die schon bekannten $y_{i}$ eingesetzt werden.

Rückwärtseinsetzen

Beim Rückwärtseinsetzen berechnet man die Lösung des ursprünglichen Gleichungssystems, indem man Rx = y ähnlich wie beim Vorwärtseinsetzen löst. Der Unterschied besteht darin, dass man bei $x_n = \frac{y_n}{r_{nn}}$ beginnt und dann nacheinander die Werte von $x_{n-1},x_{n-2},\ldots,x_1$ berechnet. Die entsprechende Formel lautet

$x_i = \frac{1}{r_{ii}} \left( y_i - \sum_{k=i+1}^{n} r_{ik} \cdot x_k \right)$

Algorithmus in Pseudocode

Der folgende Algorithmus führt eine LR-Zerlegung der Matrix A ohne Pivotisierung aus, indem er simultan L und R anhand von A erzeugt:

   Eingabe: Matrix A

   // Initialisierung
   R := A
   L := E_n

   // n-1 Iterationsschritte 
   for i := 1 to n-1
     // Zeilen der Restmatrix werden durchlaufen
     for k := i+1 to n
       // Berechnung von L
       L(k,i) := R(k,i) / R(i,i) // Achtung: vorher Prüfung auf Nullwerte notwendig
       // Spalten der Restmatrix werden durchlaufen
       for j := i to n
         // Berechnung von R 
         R(k,j) := R(k,j) - L(k,i) * R(i,j)

   Ausgabe: Matrix L, Matrix R

Final beschreiben L und R gerade die Faktorisierung von A in gewünschten Form, d.h. A=LR für eine linke untere Dreiecksmatrix L und eine rechte obere Dreiecksmatrix R.

Alternativ ist (aus möglichem Interesse an Speichereffizienz) eine simultane Entwicklung von L und R direkt in A möglich, welcher durch folgenden Algorithmus beschrieben wird:

   Eingabe: Matrix A

   // n-1 Iterationsschritte 
   for i := 1 to n-1
     // Zeilen der Restmatrix werden durchlaufen
     for k := i+1 to n
       // Berechnung von L
       A(k,i) := A(k,i) / A(i,i) // Achtung: vorher Prüfung auf Nullwerte notwendig
       // Spalten der Restmatrix werden durchlaufen
       for j := i+1 to n
         // Berechnung von R 
         A(k,j) := A(k,j) - A(k,i) * A(i,j)

   Ausgabe: Matrix A (in modifizierter Form)

Unvollständige Zerlegungen

Die LR-Zerlegung hat den Nachteil, dass sie auch bei dünnbesetzten Matrizen häufig vollbesetzt ist. Werden dann statt aller Einträge nur jene in einem vorgegebenen Besetzungsmuster berechnet, spricht man von einer unvollständigen LU-Zerlegung. Diese liefert eine günstige Approximation an die Matrix und kann somit als Vorkonditionierer bei der iterativen Lösung linearer Gleichungssysteme eingesetzt werden. Im Fall symmetrisch positiv definiter Matrizen spricht man von einer unvollständigen Cholesky-Zerlegung.

Eigenschaften des Verfahrens

Rechenaufwand und Speicherplatzbedarf

Die Anzahl arithmetischer Operationen für die LR-Zerlegung ist bei einer $n\times n$ -Matrix ca. $\tfrac 23n^3$ . Der Aufwand für das Vorwärts- und Rückwärtseinsetzen ist quadratisch ( O(n^2) ) und daher insgesamt vernachlässigbar. Da der Aufwand kubisch mit der Dimension der Matrix wächst, kann man anhand der Rechenzeit für eine Matrix die Rechenzeit für eine andere Matrix abschätzen. Benötigt der Algorithmus also auf einem bestimmten Rechner für eine Matrix der Dimension n = 10000 etwa 10 Sekunden, so benötigt er für eine Matrix der Dimension $n = 100000 = 10 \cdot 10000$ auf demselben Rechner ungefähr $10^3 \cdot 10$ Sekunden also ca. 3 Stunden. Damit ist das gaußsche Eliminationsverfahren ein schnelles direktes Verfahren zur Lösung linearer Gleichungssysteme, für eine QR-Zerlegung benötigt man mindestens doppelt so viele Rechenoperationen. Dennoch sollte der Algorithmus nur für Gleichungssysteme kleiner bis mittlerer Dimension verwendet werden (bis etwa n = 10000 ). Für Matrizen höherer Dimension sind iterative Verfahren oft besser. Diese nähern die Lösung schrittweise an und benötigen in jedem Schritt für eine vollbesetzte Matrix O(n^2) Rechenoperationen. Die Konvergenzgeschwindigkeit solcher Verfahren hängt stark von den Eigenschaften der Matrix ab und man kann die konkret benötigte Rechenzeit nur schwer vorhersagen.

Die Rechnung kann auf dem Speicher der Matrix durchgeführt werden, so dass außer der Speicherung von selbst kein zusätzlicher Speicherbedarf entsteht. Für eine vollbesetzte Matrix der Dimension n = 1000 müsste man eine Million Koeffizienten abspeichern. Dies entspricht im IEEE 754-Format double in etwa 8 Megabyte. Bei iterativen Verfahren, die mit Matrix-Vektor-Multiplikationen arbeiten, kann allerdings eine explizite Speicherung von selbst nicht nötig sein, so dass diese Verfahren ggf. vorzuziehen sind.

Für Spezialfälle lassen sich Aufwand und Speicherplatz deutlich reduzieren, indem spezielle Eigenschaften der Matrix und ihrer LR-Zerlegung ausgenutzt werden können. So benötigt die Cholesky-Zerlegung für symmetrische positiv definite Matrizen nur die Hälfte an Rechenoperationen und Speicher. Ein anderes Beispiel sind Bandmatrizen mit fester Bandbreite , da hier die LR-Zerlegung die Bandstruktur erhält und sich so der Aufwand auf O(nm^2) reduziert. Für wenige spezielle dünnbesetzte Matrizen ist es möglich, die Besetzungsstruktur auszunutzen, so dass die LR-Zerlegung ebenfalls dünnbesetzt bleibt. Beides geht einher mit einem verringerten Speicherbedarf.

Genauigkeit

Damit die Berechnung von x ausreichend genau ist, darf zum einen die Kondition der Matrix nicht zu schlecht und die verwendete Maschinengenauigkeit nicht zu gering sein. Zum anderen benötigt man ein Lösungsverfahren, das ausreichend stabil ist. Ein guter Algorithmus zeichnet sich also durch eine hohe Stabilität aus.

Im Allgemeinen ist das Verfahren ohne Pivotisierung instabil. Daher wird meist Spaltenpivotisierung zur Lösung verwendet. Damit ist das Verfahren für die meisten Matrizen stabil durchführbar, wie insbesondere durch die Arbeiten von James H. Wilkinson nach dem Zweiten Weltkrieg klar wurde. Es lassen sich allerdings Matrizen angeben, für welche die Stabilitätskonstante exponentiell mit der Dimension der Matrix wächst. Mit vollständiger Pivotisierung lässt sich die Stabilität noch verbessern, allerdings steigt dann auch der Aufwand für die Pivotsuche auf $O(n^{3})$ , daher wird diese selten verwendet. Generell bessere Stabilität haben QR-Zerlegungen, die allerdings auch aufwändiger zu berechnen sind.

Bei strikt diagonaldominanten oder positiv definiten Matrizen (siehe auch Cholesky-Zerlegung) ist das Gauß-Verfahren stabil und ohne Pivotisierung durchführbar, es treten also keine Nullen auf der Diagonale auf.

Nachiteration

Ein praktischer Ansatz zum Ausgleich dieser Rechenungenauigkeiten besteht aus einer Nachiteration mittels Splitting-Verfahren, da über die LR-Zerlegung eine gute Näherung der Matrix A zur Verfügung steht, die leicht invertierbar ist. Dazu startet man mit der berechneten Lösung $x_{0}=x$ und berechnet in jedem Schritt das Residuum

Danach berechnet man unter Verwendung der LR-Zerlegung die Lösung z_k des Gleichungssystems

und setzt

$x_{k+1}=x_k+z_k.$

Da es meistens nur um kleine Korrekturen geht, reichen oft wenige Iterationsschritte. Im Allgemeinen ist für die Berechnung des Residuums $r_{k}$ allerdings eine höhere Genauigkeit notwendig. Reicht auch die Nachiteration nicht aus, um auf die gewünschte Genauigkeit zu kommen, bleibt nur die Wahl eines anderen Verfahrens oder eine Umformung des Problems, um eine günstigere Matrix zu erhalten, etwa eine mit kleinerer Kondition.

Die Nachiteration wird beispielsweise in der LAPACK-Routine DSGESV angewandt. In dieser Routine wird die LR-Zerlegung in einfacher Genauigkeit ermittelt und die doppelte Genauigkeit der Lösung durch Nachiteration mit doppeltgenau berechnetem Residuum erreicht.

Das Gauß-Verfahren als theoretisches Hilfsmittel

Das Gauß-Verfahren ist neben seiner Bedeutung zur numerischen Behandlung von eindeutig lösbaren linearen Gleichungssystemen auch ein wichtiges Hilfsmittel in der theoretischen linearen Algebra.

Aussagen zur Lösbarkeit des linearen Gleichungssystems

Ein lineares Gleichungssystem kann eine, mehrere oder keine Lösung haben. Bei Verwendung von vollständiger Pivotisierung bringt das Gauß-Verfahren jede Matrix auf Dreiecksstufenform, bei der alle Einträge unterhalb einer gewissen Zeile Null sind und auf der Diagonale keine Nulleinträge auftauchen. Der Rang der Matrix ergibt sich dann als Zahl der Nichtnullzeilen der Matrix. Die Lösbarkeit ergibt sich dann aus dem Zusammenspiel mit der rechten Seite: Gehören zu den Nullzeilen Nichtnulleinträge der rechten Seite, ist das Gleichungssystem unlösbar, ansonsten lösbar, wobei die Dimension der Lösungsmenge der Anzahl der Unbekannten minus dem Rang entspricht.

Beispiel

$\begin{array}{rcrcl} x & + & 4y & = & 8, \\ 3x & + & 12y & = & 24. \end{array}$

Da die zweite Gleichung ein Vielfaches der ersten Gleichung ist, hat das Gleichungssystem unendlich viele Lösungen. Bei der Elimination von x in der zweiten Gleichung verschwindet diese vollständig, übrig bleibt nur die erste Gleichung. Löst man diese nach x auf, kann man die Lösungsmenge in Abhängigkeit von y angeben:

$\begin{array}{rcl}x & = & 8 - 4y \\ L & = & \{(8 - 4y,y)^T|y\} \end{array}$

Determinante

Ferner liefert das Gauß-Verfahren eine Möglichkeit, die Determinante einer Matrix zu berechnen. Da die elementaren Zeilenumformungen die Determinante 1 haben, bis auf Zeilenvertauschungen, deren Determinante -1 ist (dies ändert jedoch nur das Vorzeichen und lässt sich daher leicht korrigieren), hat die sich ergebende obere Dreiecksmatrix dieselbe Determinante wie die ursprüngliche Matrix, kann aber wesentlich einfacher berechnet werden: Sie ist das Produkt der Diagonalelemente.

Berechnung der Inversen

Eine weitere Möglichkeit der Anwendung des Gauß-Verfahrens besteht in der Berechnung der Inversen der Matrix. Hierzu wird der Algorithmus auf ein von rechts durch eine Einheitsmatrix erweitertes Schema angewandt und nach der ersten Phase fortgesetzt, bis links eine Einheitsmatrix erreicht ist. Im rechten Teil steht dann die inverse Matrix. Dieses Verfahren ist numerisch nicht zu empfehlen und die explizite Berechnung der Inversen kann meist umgangen werden.

Geschichte

Bereits im chinesischen Mathematikbuch Jiu Zhang Suanshu (dt. Neun Bücher arithmetischer Technik), das zwischen 200 vor und 100 nach Christus verfasst wurde, findet sich eine beispielhafte, aber klare Demonstration des Algorithmus anhand der Lösung eines Systems mit drei Unbekannten. 263 veröffentlichte Liu Hui einen umfassenden Kommentar zu dem Buch, der daraufhin in den Textkorpus einging. Das Jiu Zhang Suanshu war bis ins 16. Jahrhundert eine wesentliche Quelle der mathematischen Bildung in China und umliegenden Ländern.

In Europa wurde erst 1759 von Joseph-Louis Lagrange ein Verfahren publiziert, das die grundlegenden Elemente enthält. Carl Friedrich Gauß beschäftigte sich im Rahmen seiner Entwicklung und Anwendung der Methode der kleinsten Quadrate mit linearen Gleichungssystemen, den dort auftretenden Normalgleichungen. Seine erste Veröffentlichung zu dem Thema stammt von 1810 (Disquisitio de elementis ellipticis Palladis), allerdings erwähnt er bereits 1798 in seinen Tagebüchern kryptisch, er habe das Problem der Elimination gelöst. Sicher ist, dass er das Verfahren zur Berechnung der Bahn des Asteroiden Pallas zwischen 1803 und 1809 nutzte. In den 1820ern beschrieb er das erste Mal etwas wie eine LR-Zerlegung. Das Eliminationsverfahren wurde in der Folgezeit vor allem in der Geodäsie eingesetzt, und so ist der zweite Namensgeber des Gauß-Jordan-Verfahrens nicht etwa der Mathematiker Camille Jordan, sondern der Geodät Wilhelm Jordan.

Im und nach dem Zweiten Weltkrieg gewann die Untersuchung numerischer Verfahren an Bedeutung und das Gauß-Verfahren wurde nun auch vermehrt auf Probleme unabhängig von der Methode der kleinsten Quadrate angewandt. John von Neumann und Alan Turing definierten die LR-Zerlegung in der heute üblichen Form und untersuchten das Phänomen der Rundungsfehler. Befriedigend gelöst wurden diese Fragen erst in den 1960ern durch James Hardy Wilkinson, der zeigte, dass das Verfahren mit Pivotisierung rückwärtsstabil ist.

Basierend auf einem Artikel in:

Wikipedia.de