# 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

## 3. Umformung in eine Iteration

Dieser Teil folgt in Kürze.

Wenn wir die Matrix A so in das Produkt L·U zerlegen, dass L eine untere Dreiecksmatrix mit Einsen auf der Hauptdiagonalen ist und U eine obere Dreiecksmatrix, können wir die relevaten Informationen in eine einzige Matrix packen:

\def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \pmatrix{\green{u_{11}} & \green{u_{12}} & \green{u_{13}} & \ldots & \green{u_{1n}} \cr \pink{l_{21}} & \green{u_{22}} & \green{u_{23}} & \ldots & \green{u_{2n}} \cr \pink{l_{31}} & \pink{l_{32}} & \green{u_{33}} & \ldots & \green{u_{2n}} \cr \vdots & \vdots & \vdots & \ddots & \vdots \cr \pink{l_{n1}} & \pink{l_{n2}} & \pink{l_{n3}} & \ldots & \green{u_{nn}}}

Daraus lassen sich L und U gewinnen:

- L, indem die Diagonalelemente (u11 bis unn) durch 1 und die übrigen grün unterlegten Elemente durch Null ersetzt werden,
- U, indem alle rötlich unterlegten Elemente durch Null ersetzt werden.

Wir entwerfen nun einen Algorithmus, der die Ausgangsmatrix A direkt in diese Kompaktdarstellung umformt.
Im letzten Abschnitt haben wir die Matrix so zerlegt:

\def\yellow{\bbox[#ffffa0,2pt]} \def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \def\blue{\bbox[#e8f0ff,2pt;border:1px #8090ff solid]} \def\white{\bbox[border:1px #8090ff solid]} A = \pmatrix{a_{11} & \textbf w \cr \textbf v & A'} = \pmatrix{ 1 & \textbf 0 \cr \pink{\frac{1}{a_{11}} \textbf{v}} & I_{n-1}} \pmatrix{\green {a_{11}} & \green{\textbf w} \cr \textbf 0 & \yellow{A' - \frac{1}{a_{11}} \textbf{v w}}}

Wenn wir die relevanten Informationen zusammenpacken, erhalten wir:

\def\yellow{\bbox[#ffffa0,2pt]} \def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \pmatrix{\green {a_{11}} & \green{\textbf w} \cr \pink{\frac{1}{a_{11}} \textbf{v}} & \yellow{A' - \frac{1}{a_{11}} \textbf{v w}}}

Im den rötlichen und grünen Bereichen stehen schon Teile der Dreicksmatrizen, der gelb unterlegte Bereich muss weiter zerlegt werden.
Dies bekommen wir durch folgende Operationen aus der Ausgangsmatrix (machen Sie sich dazu bitte klar, wie das Produkt vw gebildet wird):

- Die erste Zeile bleibt unverändert
- Der Rest der ersten Spalte wird durch a11 dividiert
- Im Rest der Matrix wird von jedem Element das Produkt aus dem links davon stehenden Element der ersten Spalte (das schon durch a11 dividiert ist!) und dem darüber stehenden Element der ersten Zeile subtrahiert.

In mathGUIde (Matrizen werden hier ab Null indiziert!) sieht das so aus:

    for i in fromTo(1, n-1):
        A[i,0] /= A[0,0]
    for i in fromTo(1, n-1):
        for j in fromTo(1, n-1):
            A[i,j] -= A[i,0] * A[0,j]

Jetzt ist es nur noch ein kleiner Schritt, dieses Verfahren auf den Rest der Matrix zu erweitern. Dazu führen wir eine äußere Schleifenvariable `k` für den linken (und oberen) Index der zu bearbeitenden Restmatrix ein:

    def luIter(A):
        n = A.height()
        for k in fromTo(0, n-2):
            assert A[k,k] != 0 # sonst gescheitert!
            for i in fromTo(k+1, n-1):
                A[i,k] /= A[k,k]
                for j in fromTo(k+1,n-1):
                    A[i,j] -= A[i,k] * A[k,j]

Wir testen die Funktion am Beispiel aus dem ersten Abschnitt:

\def\yellow{\bbox[#ffffa0,2pt]} \def\pink{\bbox[#ffe0e0,2pt]} \def\green{\bbox[#c0ffc0,2pt]} \yellow{\pmatrix{1 & 2 & 3 \cr 2 & 9 & 10 \cr 3 & 26 & 24 }} = \pink{\pmatrix{1 & 0 & 0 \cr 2 & 1 & 0 \cr 3 & 4 & 1 }} \green{\pmatrix{1 & 2 & 3 \cr 0 & 5 & 4 \cr 0 & 0 & -1 }}

    A = Matrix.fromString("1, 2, 3;  2, 9, 10;  3, 26, 24")
    luIter(A)
    print(A)

Bitte überzeugen Sie sich, dass nun in A die gesuchte linke und die rechte Dreiecksmatrix stecken!

Fortsetzung: 4. Beseitigung des Schönheitsfehlers: LUP-Zerlegung