Incomplete LU decomposition

Wir erinnern uns an die klassische LR Zerlegung:

def LR(A):
   n = A.shape[0]
   LU = A.copy()

Wir iterieren über alle Spalten und berechnen:

for i in range(0,n-1):
    for j in range(i+1,n):
        LU[j,i] = LU[j,i]/LU[i,i]
    for j in range(i+1,n):
        for k in range(i+1,n):
            LU[j,k] = LU[j,k] - LU[j, i] * LU[i,k]
return LU

Das können wir nun nutzen um die LR Zerlegung zu berechnen:

A = numpy.array(
   [[ 3.0, 0.0, 0.0, -2.0, -1.0],
    [ 0.0, 3.0, 0.0,  0.0,  0.0],
    [-2.0, 0.0, 4.0,  2.0,  0.0],
    [-2.0, 0.0, 3.0,  4.0,  0.0],
    [-1.0, 0.0, 0.0,  0.0,  3.0]])
print "Matrix A"
print A
print "Normale Version"
print LR(A)

Ungefähr so ließe sich eine LR-Zerlegung in fast jeder Programmiersprache schreiben. Man kann dies zwar auch kompakter schreiben, aber für die spätere Implementierung der ILU Faktorisierung, ist diese kompakte Schreibweise leider nur bedingt anwendbar. Trotzdem wollen wir sie hier der Vollständigkeit halber einmal beschreiben. Wir nutzen hierbei wieder spezielle Features von python, insbesondere Index-Ranges:

def LR_compact(A):
   n = A.shape[0]
   LU = A.copy()
   for i in range(0,n-1):

mit Hilfe von Index-Ranges sparen wir uns die Schleife über die Spalteneinträge und aktualisieren direkt die komplette Sub-Matrix:

       LU[i+1:, i] = LU[i+1:,i]/LU[i,i]
       LU[i+1:, i+1:] = LU[i+1:, i+1:] - numpy.outer(LU[i+1:, i], LU[i, i+1:])
   return LU

print "Kompakte Version"
print LR_compact(A)

Jetzt wollen wir die Implementierung so anpassen, dass keine neuen Nicht-Null-Einträge erzeugt werden. Dazu müssen wir vorm schreiben prüfen, dass der Eintrag auch wirklich nicht Null ist:

def ILU(A):
   n = A.shape[0]
   LU = A.copy()
   for i in range(0,n-1):
       for j in range(i+1,n):
           if (A[i,j] != 0):
               LU[j,i] = LU[j,i]/LU[i,i]
       for j in range(i+1,n):
           for k in range(i+1,n):
               if (A[j,k] != 0):
                   LU[j,k] = LU[j,k] - LU[j, i] * LU[i,k]
   return LU

Wir wir in unserem Test jetzt sehen, erzeugen wir wirklich keine neuen Einträge:

print "ILU(A)"
print ILU(A)

Nun wollen wir diese Implementierung nutzen, um wirklich die Struktur einer dünn-besetzten Matrix zu nutzen:

S = sparse.csr_matrix(A)
print "ILU(S)"
print ILU(S).todense()

und erhalten direkt von python die Warnung, dass wir in der CSR Matrix neue Eintrage erzeugen. Aber warum? Das Problem ist, dass wir in unserer Implementierung auf Matrix-Einträge lesend zugreifen, welche nicht in der Matrix gespeichert sind. In manchen Programmiersprachen führt dies unmittelbar zu einem Fehler, in python erhalten wir eine Warnung und der Eintrag wird angelegt und mit 0 initialisiert. D.h. wir müssen zwei Dinge ändern:

def ILU_sparse(A):
   n = A.shape[0]
   LU = A.copy()

erstens müssen wir zum Testen, ob ein Eintrag 0 ist, nicht den Wert prüfen, sondern die Struktur Informationen. Wir prüfen, ob einer der nonzero Indizes zur Zeile i den Wert j hat:

def is_nonzero(i,j):
    return (A[i].nonzero()[1] == j).any()

Anschliessend überprüfen wir anhand dieses Tests, ob wir schreiben dürfen und zusätzlich ob wir lesen dürfen:

   for i in range(0,n-1):
       for j in range(i+1,n):
           if is_nonzero(j,i):
               LU[j,i] = LU[j,i]/LU[i,i]
       for j in range(i+1,n):
           for k in range(i+1,n):
               if is_nonzero(j,k) and is_nonzero(j,i) and is_nonzero(i,k):
                   LU[j,k] = LU[j,k] - LU[j, i] * LU[i,k]
   return LU

print "sparse ILU"
print ILU_sparse(S).todense()

Am besten würde man hier zwischen verschiedenen Varianten für unterschiedliche Sparse Datenstrukturen unterscheiden. Z.B. ist die aktuelle Variante für CSR Matrizen nicht gut geeignet, da diese einen effizienten zeilenweisen Zugriff erlauben, wir aber spaltenweise zugreifen.