LU-Zerlegung: Vor- / Rückwärtseinsetzen¶

Ziel¶

In dieser Lektion implementieren wir das Lösungsverfahren \(A x = b\) unter der Annahme, dass eine \(LR\)-Zerlegung \(A = L R\) der Matrix \(A\) bereits gegeben ist.

Abgabe: Sie geben als Praktikumsbericht ihren Code Python ab. Für die Bewertung ist ausschlaggebend, dass Sie sich ersichtlich mit der Aufgabe auseinandergesetzt haben, unabhängig davon, ob der Code am Ende läuft oder nicht. Wir werden allerdings im nachfolgenden Praktikum darauf aufbauen.

Theorie¶

Lineare Gleichungssysteme \(Ax = b\) für \(A \in \mathbb{R}^{n \times n}, b \in \mathbb{R}^n\) können durch Vorwärts- bzw. Rückwärtseinsetzen gelöst werden, wenn die Matrix \(A=R\) eine obere bzw. \(A=L\) eine untere Dreiecksmatrix ist. Sie implementieren und testen in diesem Praktikum die beiden Algorithmen. Das Vorwärtseinsetzen ist Algorithmus 1 im Skript. Rückwärtseinsetzen funktioniert analog.

Je nach Kontext kann zudem eine effiziente Implementierung entscheidend sein, sowohl memory-efficient (kein unnötiges Belegen von Speicherplatz) wie auch time-efficient. Im einzelnen bedeutet das

Die LU-Zerlegung wird in-place ausgeführt, wobei die ursprüngliche Matrix \(A\) überschrieben wird. Die resultierende Matrix LR enthält im Dreieck links unten die Matrix \(L\), ohne die Diagonalelemente und im Dreieck oben rechts die Matrix R (inklusive Diagonale). Dies enspricht der Vorgehensweise im Unterricht.
Es werden keine Zeilenvertauschungen explizit durchgeführt, stattdessen wird mit einem Indexvektor gearbeitet.

Aufgaben¶

Aufgabe 1¶

Beschreiben Sie das Vorwärtseinsetzen als Pseudocode, analog zu Algorithmus \(1\) im Skript.
Beschreiben Sie in Pseudocode das Vor- und Rückwärtseinsetzen, wenn die LU-Zerlegung in einer Matrix \(LR\) gegeben ist. Beachten Sie, dass der Vektor \(x\) beim Rückwärtseinsetzen überschrieben werden kann. Drücken Sie die Operationen soweit möglich als Skalarprodukte aus.
Beschreiben Sie die in-place LU-Zerlegung ohne Zeilentausch als Pseudocode.

Aufgabe 2¶

Implementieren Sie eine Funktion fbSubs für das Vorwärts- und Rückwärtseinsetzen. Versuchen Sie, so wenig Schleifen wie möglich zu verwenden!

[ ]:

import numpy as np

"""
forwardBackwardSubs: Vorwärts- und Rückwärtseinsetzen
in:
 - Matrix LR, die das Ergebnis einer LU-Zerlegung enthält
 - Vektor b: rechte Seite des LGS
out: Lösung x des LGS
"""
def fbSubs(LR, b):
    # code

Aufgabe 3¶

Testen Sie fbSubs zuerst mit dem Zahlenbeispiel aus dem Unterricht (wo Sie alles schrittweise nachvollziehen können) und dann mit dem folgenden Testcode, den Sie auch variieren dürfen, bis Sie sicher sind, dass alles wie gewünscht funktioniert. Sie müssen sich im nachfolgenden Praktikum auf fbSubs verlassen können!

[ ]:

#Test LR
n = 7 # Dimension der Koeffizientenmatrix
for k in range(1000):                      # Testläufe
    LR = np.array( np.random.rand(n,n) )   # zufällige Matrix LR
    rhs = np.array(np.random.rand(n))    # zufällige rechte Seite des LGS
    x = fbSubs(LR, rhs)   # Aufruf der eigenen Funktion

    L,R = np.tril(LR,-1)+np.identity(n), np.triu(LR) # L und R extrahieren

    assert( np.linalg.norm(rhs - L@R@x) < 1e-10)     # Test, mit numerischer Toleranz

Implementierung Gauss-Algorithmus mit Spaltenpivotisierung¶

Ziel¶

In diesem Praktikum implementieren wir das Lösungsverfahren \(A x = b\) inklusive \(LU\)-Zerlegung \(P A = L R\) der Matrix \(A\) mit Spaltenpivotisierung.

Abgabe: Sie geben als Praktikumsbericht ihren Code in Python ab. Für die Bewertung ist ausschlaggebend, dass Sie sich ersichtlich mit der Aufgabe auseinandergesetzt haben, unabhängig davon, ob der Code am Ende läuft oder nicht. Sie können den Code direkt in das notebook schreiben oder in einer Entwicklungsumgebung Ihrer Wahl.

Theorie¶

Ein lineares Gleichungssystem \(Ax = b\) kann gelöst werden, indem die Matrix \(A\) in eine untere bzw. obere Dreiecksmatrix zerlegt wird: \(A = L R\), falls keine Zeilenvertauschungen notwendig sind bzw. \(P A = L R\) mit der zusätzlichen Permutationsmatrix \(P\), falls Zeilen vertauscht werden.

Diesen Prozess nennt man auf Deutsch \(LR\)-Zerlegung (englisch: \(LU\)-factorization, für lower, upper). Der Algorithmus ist eine kleine Erweiterung des Gauss-Algorithmus (Algorithmus 2 im Skript). \(R\) ist die gewohnte Zeilenstufenform, in \(l_{ij}\) merken wir uns zusätzlich, mit welchem Faktor Zeile \(i\) von Zeile \(j\) subtrahiert wurde. Die Beschreibung des Vorgehens entnehmen Sie dem Unterricht.

Sind Zeilenvertauschungen notwendig (oder auf Grund einer Pivotstrategie erwünscht), so führt die Permutationsmatrix \(P\) Buch über diese Vertauschungen: alle Zeilenvertauschungen in \(A\) werden auch in der Matrix \(P\) übernommen.

Je nach Kontext kann zudem eine effiziente Implementierung entscheidend sein, sowohl memory-efficient (kein unnötiges Belegen von Speicherplatz) wie auch time-efficient. Im einzelnen bedeutet das

Die LU-Zerlegung wird in-place ausgeführt, wobei die ursprüngliche Matrix \(A\) überschrieben wird. Die resultierende Matrix LR enthält im Dreieck links unten die Matrix \(L\), ohne die Diagonalelemente und im Dreieck oben rechts die Matrix R (inklusive Diagonale). Dies enspricht der Vorgehensweise im Unterricht.
Zeilenvertauschungen können formal durch eine Permutationsmatrix \(P\) erfasst werden. Es gilt dann \(L R = P A\) und man hat wegen \(A x = b\) <=> \(P A x = P b\) die beiden LGS

\(L \cdot y = P \cdot b\) und \(R \cdot x = y\)

zu lösen. Um die Matrixmultiplikation mit der Permutationsmatrix zu sparen arbeitet man stattdessen mit einem Indexvektor \(i\), der zu Beginn auf \(i = [0, 1, ..., n-1]\) initialisiert wird. Werden Zeilen in \(A\) getauscht, so tauscht man auch die entsprechenden Einträge in \(i\) und löst am Ende die beiden LGS

\(L \cdot y = b[i]\) und \(R \cdot x = y\)

Ziel dieses Praktikums ist eine Funktion linsolve, die LGS \(A \cdot x = b\) mit Hilfe einer in-place LU-Zerlegung mit Spaltenpivotisierung löst. Wir verwenden die Ergebnisse des letzten Prakikums:

Aufgaben¶

Aufgabe 1¶

Implementieren Sie die LU-Zerlegung: die Funktion LU nimmt als Input eine quadratische Matrix \(A\) und gibt Dreiecksmatrizen \(L, R\) zurück, so dass \(L R = A\) gilt. Im ersten Schritt lassen wir die Pivotisierung weg. Der Indexvektor idx bleibt dann unverändert.

[ ]:

# LU-Zerlegung der quadratischen Matrix A
# in: quadratische Matrix A
#out:
# - A wird überschrieben, unteres Dreieck = L (ohne Diagonale), oberes Dreieck = R
# - idx: Indexvektor der Zeilenvertauschungen
def LU(A):
    m = A.shape[0]
    idx = np.array(range(m))
    # code
    return A, idx

Testen Sie LU mit zufällig erzeugten Matrizen.

[ ]:

n = 7
#test LU
for k in range(1000):
    A = np.array( np.random.rand(n,n) )   # zufällige Matrix A erzeugen
    LR, idx = LU(A.copy())                  # LU-Zerlegung von A
    L,R = np.tril(LR,-1)+np.identity(n), np.triu(LR)  # Matrizen L, R extrahieren
    assert( np.linalg.norm(L@R - A[idx]) < 1e-8)

Aufgabe 2¶

Erstellen Sie eine Funktion linsolve(A, b), die für eine Matrix \(A\) und einen Vektor \(b\) das LGS \(A x = b\) mit Hilfe der LR-Zerlegung und dem Vorwärts- und Rückwärtseinsetzen löst.

[ ]:

# lineares Gleichungssystem A*x = b lösen.
def linsolve(A, b):
    # code

Testen Sie linsolve mit zufällig erzeugten Matrizen

[ ]:

#test linsolve
for k in range(1000):
    A = np.random.rand(n,n)
    rhs = np.random.rand(n)
    x = linsolve(A.copy(), rhs)
    assert( lin.norm(rhs - A @ x) < 1e-10)

Mit der folgenden Funktion rndCond(n, cond) können Sie zufällige \(n \times n\) Matrizen mit vorgegebener Konditionszahl cond erzeugen.

[ ]:

import numpy.linalg as lin
import numpy.random as rnd

# random orthogonal matrix
def rndOrtho(n):
    S = rnd.rand(n,n)
    S = S - S.T
    O = lin.solve(S - np.identity(n), S + np.identity(n))
    return O

# random matrix with specified condition number
def rndCond(n, cond):
    d = np.logspace(-np.log10(cond)/2, np.log10(cond)/2,n);
    A = np.diag(d)
    U,V = rndOrtho(n), rndOrtho(n)
    return U@A@V.T

Verwenden Sie diese, um Ihre Implementation von linsolve zu testen. Berechnen Sie den relativen Fehler \(\frac{|\Delta_x|}{|x|}\) und vergleichen Sie mit dem Ergebnis von linsolve aus der Bibliothek numpy.linalg.

[ ]:

for k in range(N):
    A = rndCond(n, 1e14)
    # code

Aufgabe 3¶

Implementieren Sie die Spaltenpivotisierung und wiederholen Sie die Testläufe aus Aufgabe 2