# LUP-Zerlegung

1. Lösung linearer Gleichungen mit LU-Zerlegung
2. Ein einfacher rekursiver Algorithmus mit einem Schönheitsfehler
3. Umformung in eine Iteration
4. Beseitigung des Schönheitsfehlers: LUP-Zerlegung

Diese Präsentation ist unfertig. Sie zeigt einige Möglichkeiten des TEX-Formelsatzes in mathGUIde (unter Verwendung von jsMath).

## 2. Ein einfacher rekursiver Algorithmus mit einem Schönheitsfehler

Wir gliedern die quadratische n·n-Matrix

A = \pmatrix{{a_{11}} & a_{12} & \ldots & a_{1n} \cr a_{21} & a_{22} & \ldots & a_{2n} \cr \vdots & \vdots & \ddots & \vdots \cr a_{n1} & a_{n2} & \ldots & a_{nn}}
in vier Teilmatrizen:
\def\blue{\bbox[#e8f0ff,2pt;border:1px #8090ff solid]} \def\white{\bbox[border:1px #8090ff solid]} A = \pmatrix{ \white {\matrix{a_{11}}} & \blue {\matrix{a_{12} & \ldots & a_{1n}}} \cr \blue {\matrix{a_{21} \cr \vdots \cr a_{n1}}} & \white {\matrix{a_{22} & \ldots & a_{2n} \cr \vdots & \ddots & \vdots \cr a_{n2} & \ldots & a_{nn}}}} = \pmatrix{a_{11} & \textbf w \cr \textbf v & A'}

Die Elemente dieser Zerlegung (statt A' schreiben wir `A1`) können wir in mathGUIde so schreiben:

    v = A.submatrix(1,0, n-1,1)
    w = A.submatrix(0,1, 1,n-1)
    A1 = A.submatrix(1,1, n-1,n-1)

Zum Verständnis: Die Matrizen-Indizes zählen in mathGUIde nicht ab 1 sondern ab Null (um mit den Python-Listen kompatibel zu bleiben). Der `submatrix`-Methode übergibt man in den ersten beiden Parametern die Zeile und Spalte des linken oberen Elements der Untermatrix, dann die Anzahl der zu übernehmenden Zeilen und Spalten.

Durch Ausmultiplizieren sieht man, dass diese Matrix in folgendes Produkt zerlegt werden kann, vorausgesetzt, dass a11 von Null verschieden ist:

\def\yellow{\bbox[#ffffa0,2pt]} A = \pmatrix{1 & \textbf 0 \cr \frac{1}{a_{11}} \textbf{v} & I_{n-1}} \pmatrix{a_{11} & \textbf w \cr \textbf 0 & \yellow{A' - \frac{1}{a_{11}} \textbf{v w}}}
Wir ziehen uns jetzt an den eigenen Haaren aus dem Sumpf, indem wir das Problem auf eine kleinere Matrix abwälzen, nämlich auf die gelb unterlegte (n-1)·(n-1)-Matrix. Nehmen wir also einmal an, dass wir die als Produkt einer unteren (linken) und einer oberen (rechten) Dreiecksmatrix darstellen können:
\def\yellow{\bbox[#ffffa0,2pt]} B = \yellow{A' - \frac{1}{a_{11}} \textbf{v w}} = L' U'
Damit erhalten wir
\def\yellow{\bbox[#ffffa0,2pt]} \def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} A = \pmatrix{1 & \textbf 0 \cr \frac{1}{a_{11}} \textbf v & I_{n-1}} \pmatrix{a_{11} & \textbf w \cr \textbf 0 & \yellow{L'U'}} = \pink{\pmatrix{1 & \textbf 0 \cr \frac{1}{a_{11}} \textbf{v} & L'}} \green{\pmatrix{a_{11} & \textbf w \cr \textbf 0 & U'}}
Man sieht leicht, dass die rötlich unterlegte Matrix eine untere Dreiecksmatrix ist und die grün unterlegte eine obere Dreiecksmatrix.

Vorausgesetzt, dass wir nicht auf eine Division durch Null stoßen (das ist ein wunder Punkt, den wir noch klären müssen), können wir also durch rekursive Anwendung dieser Zerlegung eine n·n-Matrix in eine untere und eine obere Dreiecksmatrix zerlegen.
Wir müssen nur dafür sorgen, dass die Rekursion abbricht. Dazu behandeln wir die 1·1-Matrizen separat: Jede 1·1-Matrix ist ja immer gleichzeitig eine untere und auch eine obere Dreiecksmatrix, da sie nur aus der Hauptdiagonale besteht.
Wir können dann einfach schreiben:

\def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \pmatrix{a_{11}} = \pink{\pmatrix{1}} \green{\pmatrix{a_{11}}}

Diese Rekursion können wir nun fast direkt in eine mathGUIde-Funktion übernehmen. Die farbig unterlegten Programmzeilen entsprechen den farbigen Unterlegungen in den Formeln:

    def luRecursive(A):
        n = A.height()
        if n == 1: # Rekursionsabbruch: 1·1-Matrix
            return Matrix.identity(1), A.copy()
        v = A.submatrix(1,0, n-1,1)
        w = A.submatrix(0,1, 1,n-1)
        A1 = A.submatrix(1,1, n-1,n-1)
        B = A1 - (1/A[0,0]) * v * w
        L1, U1 = luRecursive(B) # Rekursiver Aufruf für (n-1)·(n-1)-Matrix
        L_top = Matrix.identity(1).joinRight(Matrix.null(1,n-1))
        L_btm = ((1/A[0,0]) * v).joinRight(L1)
        L = L_top.joinBottom(L_btm)
        U_top = A.submatrix(0,0, 1,1).joinRight(w)
        U_btm = Matrix.null(n-1,1).joinRight(U1)
        U = U_top.joinBottom(U_btm)
        return L,U

Wir testen den Algorithmus mit einer 3·3-Matrix:

    A = Matrix.fromString("5,4,2; 1,9,7; 3,0,6")
    L, U = luRecursive(A)
    print(L, U)

Zur Kontrolle vergleichen wir nun die Originalmatrix mit dem Produkt der beiden Dreiecksmatrizen:

    print(A, L*U)

Fortsetzung: 3. Umformung in eine Iteration