Beispiel

1. ${\displaystyle x_{1}+2x_{2}+3x_{3}=2}$, hier: ${\displaystyle a_{11}=1}$, ${\displaystyle a_{12}=2}$, ${\displaystyle a_{13}=3}$ und ${\displaystyle b_{1}=2}$
2. ${\displaystyle x_{1}+x_{2}+x_{3}=2}$
3. ${\displaystyle 3x_{1}+3x_{2}+x_{3}=0}$

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

${\displaystyle \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 ${\displaystyle a_{21}}$ und ${\displaystyle 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 ${\displaystyle a_{21}=1}$) durch das Pivotelement ${\displaystyle a_{11}}$ teilt (hier: ${\displaystyle {\tfrac {1}{1}}=1}$ und ${\displaystyle {\tfrac {3}{1}}=3}$). Da die beiden Elemente ${\displaystyle a_{21}}$ und ${\displaystyle a_{31}}$ Null werden sollen, werden die beiden Multiplikatoren jeweils mit ${\displaystyle (-1)}$ multipliziert.

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

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

${\displaystyle -2x_{3}=-6}$.

Diese Gleichung ist einfach lösbar und liefert ${\displaystyle x_{3}=3}$. Damit ergibt sich für die zweite Zeile

${\displaystyle -1x_{2}-2x_{3}=0}$, mit ${\displaystyle x_{3}=3}$ also ${\displaystyle x_{2}=-6}$

und weiter ${\displaystyle x_{1}=5}$. Damit sind alle Variablen berechnet:

${\displaystyle x_{1}=5}$, ${\displaystyle x_{2}=-6}$ und ${\displaystyle x_{3}=3}$.

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 angeschrieben. Für die erste Zeile ist die Zeilensumme ${\displaystyle 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 ${\displaystyle (-1)}$-fache der ersten addiert. Macht man das auch für die Zeilensumme, so gilt ${\displaystyle 5+(-1)\cdot 8=-3}$. Dieses Ergebnis ist die Zeilensumme der umgeformten zweiten Zeile ${\displaystyle -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.

Einführendes Beispiel

Will man das Lösen eines quadratischen eindeutig lösbaren Gleichungssystems ${\displaystyle 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 ${\displaystyle A}$ in das Produkt einer linken unteren, normierten Dreiecksmatrix ${\displaystyle L}$ (links, bzw. Englisch „left“, oder auch „lower“) und einer rechten oberen Dreiecksmatrix ${\displaystyle R}$ (rechts, bzw. Englisch „right“, oder auch „upper“, und dann mit ${\displaystyle U}$ bezeichnet). Das folgende Beispiel zeigt dies:

${\displaystyle 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 dient die Matrix ${\displaystyle L}$ dem Speichern der benötigten Umformungsschritte, die Multiplikationen mit Frobeniusmatrizen entsprechen, und ${\displaystyle R}$ hat die oben erwähnte Stufenform. Das zeigt die Existenz der Zerlegung. Um Eindeutigkeit zu erreichen, werden die Diagonalelemente der Matrix ${\displaystyle L}$ als 1 festgelegt. Die Umformungsschritte zu speichern hat den Vorteil, dass für verschiedene „rechte Seiten“ ${\displaystyle b}$ 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 ${\displaystyle P}$ beschrieben werden:

${\displaystyle P\cdot A=L\cdot R.}$

Existenzsatz

Für jede reguläre Matrix ${\displaystyle A\in \mathbb {R} ^{n\times n}}$ existiert eine Permutationsmatrix ${\displaystyle P\in \mathbb {R} ^{n\times n}}$, eine untere, normierte Dreiecksmatrix ${\displaystyle L\in \mathbb {R} ^{n\times n}}$ und eine obere Dreiecksmatrix ${\displaystyle R\in \mathbb {R} ^{n\times n}}$, sodass gilt:

${\displaystyle P\cdot A=L\cdot R}$.

Eine Permutationsmatrix ${\displaystyle P}$ ist eine Matrix, die aus der Einheitsmatrix durch eine beliebige Anzahl an Zeilenvertauschungen entsteht und somit weiterhin nur aus Nullen und Einsen besteht.

Algorithmus

Der Algorithmus zur Berechnung der Matrizen ${\displaystyle P,L,R}$ für ein vorgegebenes ${\displaystyle A\in \mathbb {R} ^{n\times n}}$ lautet wie folgt.

Es werden ${\displaystyle n}$ Matrixumformungen vollzogen ( ${\displaystyle k=0,\ldots ,n-1}$). Dabei führt man die Umformungsmatrizen ${\displaystyle A^{(k)}}$ ein:

${\displaystyle A^{(k)}=\left(a_{ij}^{(k)}\right)={\begin{cases}A&,\quad k=0\\(I-L^{(k)})P^{(k-1)}A^{(k-1)}&,\quad k=1,\ldots ,n-1.\end{cases}}}$

Dabei wurden neue Hilfsmatrizen ${\displaystyle L^{(k)},P^{(k)}}$ eingeführt:

${\displaystyle {\begin{aligned}\left(L^{(k)}\right)_{ij}&={\begin{cases}{\frac {a_{ik}^{(k-1)}}{a_{kk}^{(k-1)}}}&,\quad j=k\;\wedge \;i>k\\0&,\quad {\text{sonst}}\end{cases}}\\P^{(k)}&=(e_{1},\ldots ,e_{k-1},e_{\hat {k}},e_{k+1},\ldots ,e_{{\hat {k}}-1},e_{k},e_{{\hat {k}}+1},\ldots ,e_{n})\\{\hat {k}}&\in \{k,\ldots ,n\}\quad {\text{sodass gilt}}\quad \left|a_{{\hat {k}}k}^{(k)}\right|=\max {\left\lbrace \left|a_{ik}^{(k)}\right|\colon i=k,\ldots ,n\right\rbrace }.\end{aligned}}}$

Beachte:

• ${\displaystyle e_{i}\in \mathbb {R} ^{n}}$ stellt den ${\displaystyle i}$-ten Basisvektor dar
• in ${\displaystyle P^{(k)}}$ wurde nur eine Zeile der Einheitsmatrix mit einer anderen vertauscht
• ${\displaystyle {\hat {k}}}$ muss so gewählt werden, dass ${\displaystyle a_{{\hat {k}}k}^{(k)}}$ den betragsmäßig größten Wert für alle Elemente aus der Teilspalte ${\displaystyle k}$ hat

Man benötigt noch weitere Hilfsmatrizen ${\displaystyle Q^{(k)}}$:

${\displaystyle Q^{(k)}={\begin{cases}I&,\quad k=n\\P^{(n-1)}\cdot \ldots \cdot P^{(k)}&,\quad k<n.\end{cases}}}$

Nun können die gewünschten Matrizen angegeben werden:

${\displaystyle {\begin{aligned}R&=A^{(n)}\\P&=Q^{(1)}\\L&=\left(I+\sum _{k=1}^{n-1}Q^{(k+1)}L^{(k)}\right).\end{aligned}}}$

Algorithmus in Pseudocode

   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

   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)

   Eingabe: Matrix 
      ${\displaystyle A}$

   // n-1 Iterationsschritte
   for k := 1 to n-1
     // Pivotisierung
     // finde betragsmäßig größtes Element in k-ter Spalte
     
      ${\displaystyle {\hat {k}}}$ = k
     
      ${\displaystyle {\hat {a}}=|a_{{\hat {k}}k}^{(k-1)}|}$
     for i := k+1 to n
       if < 
      ${\displaystyle \left(|a_{ik}^{(k-1)}|>{\hat {a}}\right)}$ >
         
      ${\displaystyle {\hat {k}}}$ = i
         
      ${\displaystyle {\hat {a}}=|a_{ik}^{(k-1)}|}$
     // vertausche Zeilen
     <Erstelle 
      ${\displaystyle P^{(k)}}$>
     <Vertausche Zeilen k,
      ${\displaystyle {\hat {k}}}$>
     // Rest analog zu obigen Algorithmen
     <Rest 
      ${\displaystyle \ldots }$>
   
      ${\displaystyle P=P^{(n)}\cdots P^{(1)}}$

   Ausgabe: Matrix L, Matrix R, Matrix P

Lösen eines linearen GleichungssystemS

Das ursprüngliche LGS ${\displaystyle Ax=b}$ wird mittels der LR-Zerlegung nun wie folgt vereinfacht:

${\displaystyle {\begin{aligned}Ax=b\quad &{\text{und}}\quad PA=LR\\\Rightarrow PAx&=Pb\\\Leftrightarrow LRx&=Pb.\end{aligned}}}$

Nun definiert man die folgenden Hilfsvariablen

${\displaystyle y:=Rx\quad {\text{und}}\quad {\hat {b}}:=Pb.}$

Folglich hat sich das LGS ${\displaystyle Ax=b}$ in eine vereinfachte Struktur gewandelt:

${\displaystyle Ly={\hat {b}}\quad {\text{und}}\quad Rx=y.}$

Diese können leicht durch Vorwärts- bzw. Rückwärtseinsetzen gelöst werden.

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 ${\displaystyle A}$ 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.

Rechenaufwand und Speicherplatzbedarf

Die Anzahl arithmetischer Operationen für die LR-Zerlegung ist bei einer ${\displaystyle n\times n}$-Matrix ca. ${\displaystyle {\tfrac {2}{3}}n^{3}}$. Der Aufwand für das Vorwärts- und Rückwärtseinsetzen ist quadratisch ( ${\displaystyle {\mathcal {O}}(n^{2})}$) und daher insgesamt vernachlässigbar. Das gaußsche Eliminationsverfahren ist 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 ${\displaystyle 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 ${\displaystyle {\mathcal {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 ${\displaystyle A}$ durchgeführt werden, so dass außer der Speicherung von ${\displaystyle A}$ selbst kein zusätzlicher Speicherbedarf entsteht. Für eine vollbesetzte Matrix der Dimension ${\displaystyle 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 ${\displaystyle A}$ 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 ${\displaystyle m}$, da hier die LR-Zerlegung die Bandstruktur erhält und sich so der Aufwand auf ${\displaystyle {\mathcal {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

Aussagen zur Lösbarkeit des linearen Gleichungssystems

Ein lineares Gleichungssystem kann keine Lösung (unlösbar), genau eine Lösung (eindeutig lösbar) oder unendlich viele Lösungen haben. Bei Verwendung von vollständiger Pivotisierung bringt das Gauß-Verfahren jede Koeffizientenmatrix auf eine reduzierte Stufenform. Der Rang der (ursprünglich gegebenen) Koeffizientenmatrix ist gleich der Anzahl der Nichtnullzeilen der in reduzierte Stufenform gebrachten Matrix. Die Lösbarkeit ergibt sich dann aus dem Zusammenspiel mit der rechten Seite: Gehören zu den Nullzeilen der in reduzierte Stufenform gebrachten Matrix Nichtnulleinträge der rechten Seite, ist das Gleichungssystem unlösbar, ansonsten lösbar. Die Anzahl der freien Parameter in der Lösungsmenge ist gleich der Anzahl der Unbekannten minus dem Rang. Das alles ergibt sich aus dem Satz von Kronecker-Capelli.

Beispiel

${\displaystyle {\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, der dann die Rolle eines freien Parameters spielt, angeben:

${\displaystyle {\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.