# 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

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

Wir kommen nun auf das Problem der Division durch Null (vom Ende des 1. Abschnitts) zurück. Bei der folgenden Matrix versagt unser Algorithmus, obwohl sie invertierbar ist:

A = \pmatrix{0 & 5 & 2 \cr 2 & 6 & 4 \cr 2 & 1 & 1 }

    A = Matrix.fromString("0, 5, 2;  2, 6, 4;  2, 1, 1")
    luIter(A)

Wir können die Division durch Null hier verhindern, indem wir die ersten beiden Zeilen der Matrix vertauschen:

A' = \pmatrix{2 & 6 & 4 \cr 0 & 5 & 2 \cr 2 & 1 & 1 }

    A1 = Matrix.fromString("2, 6, 4;  0, 5, 2;  2, 1, 1")
    luIter(A1)
    print(A1)

Die Vertauschung der der ersten beiden Zeilen lässt sich erreichen, indem man A von links mit der entsprechenden Permutationsmatrix multipliziert. Diese entsteht aus der Einheitsmatrix durch Vertauschung der der entsprechenden Zeilen:

A' = \pmatrix{2 & 6 & 4 \cr 0 & 5 & 2 \cr 2 & 1 & 1 } = \pmatrix{0 & 1 & 0 \cr 1 & 0 & 0 \cr 0 & 0 & 1 } \pmatrix{0 & 5 & 2 \cr 2 & 6 & 4 \cr 2 & 1 & 1 } = PA

In einer invertierbaren Matrix lassen sich die Zeilen immer so permutieren, dass unser Algortihmus anschließend erfolgreich ist. Es gilt nämlich:

Satz

Zu jeder invertierbaren Matrix A gibt es

- eine Permutationsmatrix P,
- eine untere Dreiecksmatrix L (nur Einsen auf der Hauptdiagonalen),
- eine obere Dreiecksmatrix U (keine Null auf der Hauptdiagonalen),

so dass PA = LU.

Wir können das lineare Gleichungssystem auf beiden Seiten mit P multiplizieren, und dann mit der Zerlegung der der permutierten Matrix arbeiten:

\def\green{\bbox[#d0ffd0,2pt]} A\textbf x = \textbf b \quad \Rightarrow\quad \green{PA}\textbf x = P\textbf b \quad \Rightarrow\quad \green{LU}\textbf x = P\textbf b

Dieses Gleichungssystem lässt sich in zwei Schritten lösen:

#### 1. Schritt:

Mit dem Ansatz U x = y wird die Gleichung zu L y = P b. Diese Gleichung nach y auflösen!

#### 2. Schritt:

Mit dem nun bekannten y kann der Ansatz U x = y nach x aufgelöst werden.

In der nun folgenden Implementierung vertauschen wir nicht nur, um Divisionen durch Null zu verhindern, sondern bringen immer das dem Betrage nach größte Element an die Divisionsposition. Damit erhöht sich die numerische Stabilität beim Rechnen mit Gleitkommazahlen.

    def lup(A):
        A = A.copy() # Kopie, damit übergebene Matrix unverändert bleibt
        n = A.height()
        P = Matrix.identity(n) # Ausgangsbasis für Permutationsmatrix
        for k in fromTo(0,n-2):
            # Suche das dem Betrag nach größte Element
            # der k-ten Spalte ab der k-ten Zeile:
            absMax, kMax = max([(abs(A[i,k]), i) for i in fromTo(k,n-1)])
            assert absMax > 0

            # Vertausche die Zeile des größten Elements mit der k-ten Zeile...
            P[k], P[kMax] = P[kMax], P[k]  # ... in der Permutation
            A[k], A[kMax] = A[kMax], A[k]  # ... und in der Matrix

            for i in fromTo(k+1,n-1):
                # Dividiere die restlichen Elemente
                # der k-ten Spalte durch das „Kopfelement“:
                A[i,k] /= A[k,k]
                # Subtrahiere in der Restmatrix überall das Produkt
                # aus den „Randelementen“ links und oberhalb:
                for j in fromTo(k+1,n-1):
                    A[i,j] -= A[i,k] * A[k,j]

        # A in linke und rechte Dreiecksmatrix zerlegen
        L = Matrix.identity(n)
        U = Matrix.null(n)
        for i in range(n):
            for k in range(n):
                if k < i:  L[i,k] = A[i,k]
                else:      U[i,k] = A[i,k]
        return L, U, P

Ein Test mit der Matrix von oben:

    A = Matrix.fromString("0, 5, 2;  2, 6, 4;  2, 1, 1")
    L,U,P = lup(A)
    print(L,U,P)

Und die Kontrolle:

    print(P * A)
    print(L * U)

Die Funktion `lup` ist in mathGUIde bereits eingebaut, allerdings nicht als globale Funktion sondern als Matrix-Methode:

    A = Matrix.fromString("0, 5, 2;  2, 6, 4;  2, 1, 1")
    L,U,P = A.lup()
    print(L,U,P)