Cholesky Zerlegung

Lernziele

  • Sie können einen mathematischen Algorithmus in Python umsetzen.

  • Sie können mit Hilfe der Cholesky Zerlegung Lösungen linearer Gleichungssysteme berechnen.

  • Sie kennen die beiden verschiedenen Formen der Cholesky Zerlegung.

Theorie

Die LU-Zerlegung, welche im letzten Praktikum betrachtet wurde, kann prinzipiell auf beliebige Gleichungssysteme mit regulären Matrizen angewandt werden. In vielen Anwendungen treten jedoch Matrizen auf, die zusätzliche Struktureigenschaften haben. In vielen Fällen ist die Matrix {em positiv definit} und symmetrisch. In dem Fall kann der Rechenaufwand um einen Faktor 2 reduziert werden. Zu dem ist keine Pivotisierung notwendig. Die Cholesky-Zerlegung setzt diese Eigenschaftenvoraus und wird mit unter auch zum Testen ob eine Matrix positiv definit ist, benutzt.

Aufgabe

In einem ersten Schritt dieses Praktikums wird die Cholesky-Zerlgung implementiert und getestet.

  • Implementieren Sie die Cholesky-Zerlegung gemäss Algorithmus 2.5 im Skript. Benutzen Sie dazu die folgende Funktionsdefinition und Rückgabe:

    def mychol(A):
    # Vektor für die Diagonale
    d = np.zeros(A.shape[0])
    # Matrix für die untere Dreiecksmatrix
    L = np.zeros_like(A) # oder np.eye(A.shape[0])

    <snipp selber machen>

return d,L

    Vorsicht bei der Anwendung: es wird davon ausgegangen, dass \(L\) eine normierte untere Dreiecksmatrix ist. Daher wird die Diagonale nicht mit 1 gefüllt, was abhängig vom Algorithumus für das Vorwärts- und Rückwärtseinsetzen auch nicht notwendig ist.

  • Testen Sie den Code in dem Sie

    • die Zerlegung,

    • das Residuum und

    • die Lösung testen.

    Letzteres benötigt die exakte Lösung. Sie können gemäss dem Skript unten vorgehen.

     1#Test Cholesky-Zerlegung
 2n = 7 # Dimension der Koeffizientenmatrix
 3for k in range(1000):
 4    A = np.array(np.random.rand(n,n))
 5    # positiv definite Matrix
 6    ATA = A.T@A
 7    x_ref = np.array(np.random.rand(n))
 8    # rechte Seite für gegebenes x berechnen
 9    b = ATA@x_ref
10    # Lösen des Gleichungssystems
11    d,L = mychol(ATA)
12    z = <snipp, selber machen> # Vorwaertseinsetzen
13    x = <snipp, selber machen> # Rückwärtseinsetzen
14    # Test der Zerlegung
15    assert( np.linalg.norm(ATA-L@np.diag(d)@L.T) < 1e-8)
16    # Test des Residuums
17    assert(np.linalg.norm(ATA@x-b) < 1e-8)
18    try:
19        assert( np.linalg.norm(x-x_ref)/np.linalg.norm(x_ref) < 1e-8)
20    except:
21        print('Vorsicht: relativer Fehler',np.linalg.norm(x-x_ref)/np.linalg.norm(x_ref))