ILU-Zerlegung

Als ILU-Zerlegung (von incomplete LU-Decomposition) oder unvollständige LU-Zerlegung bezeichnet man in der numerischen Mathematik die fehlerbehaftete Zerlegung einer Matrix $A\in \mathbb {R} ^{n\times n}$ in das Produkt einer unteren Dreiecksmatrix L und einer oberen Dreiecksmatrix U

$LU\approx A$,

bei der von den Zerlegungsmatrizen L und U nur die Einträge einer vorgegebenen Besetzungsstruktur berechnet werden.

Bei der Berechnung einer normalen LU-Zerlegung einer dünnbesetzten Matrix kann man die Besetzungsstruktur in der Regel nicht ausnutzen. Es wird daher sehr viel mehr Speicherplatz benötigt als für die ursprüngliche Matrix und auch die Anzahl der notwendigen Rechenoperationen ist nicht geringer als die für eine vollbesetzte Matrix. Durch die Vorgabe einer maximalen Besetzungsstruktur wird dieses Problem unter Inkaufnahme einer fehlerbehafteten Zerlegung umgangen.

Die ILU-Zerlegung wird erfolgreich als Vorkonditionierer zur Beschleunigung der iterativen Lösung großer dünnbesetzter linearer Gleichungssysteme $Ax=b$ mittels Krylow-Unterraum-Verfahren eingesetzt. Es werden dabei keine Eigenschaften des eigentlichen Problems (meist die numerische Lösung einer partiellen Differentialgleichung) ausgenutzt. Damit ist sie nicht auf bestimmte Problemklassen beschränkt und hat Einzug in viele Bereiche der numerischen Simulation gefunden, beispielsweise in der numerischen Strömungsmechanik ist die Technik weit verbreitet.

Zuerst erwähnt wurde das Verfahren 1960 von Richard S. Varga und Nikolai Iwanowitsch Bulejew (N. I. Buleev). Eine genauere Analyse wurde 1977 von J. A. Meijerink und van der Vorst veröffentlicht. Diese untersuchten Vorkonditionierungstechniken für das CG-Verfahren und schlugen eine unvollständige Cholesky-Zerlegung für symmetrische Matrizen vor. Gleichzeitig erwähnten sie eine Erweiterung auf allgemeine Matrizen.

Anwendung

Für die Anwendung als Vorkonditionierer wird das Gleichungssystem $Ax=b$ formal mit $(LU)^{-1}$ multipliziert,

$(LU)^{-1}Ax=(LU)^{-1}b.$

Wendet man darauf Krylow-Unterraum-Verfahren an, so konvergieren diese besser, da die Matrix $(LU)^{-1}A$ näher an der Einheitsmatrix als A ist. Für diese Verfahren benötigt man in jedem Schritt Matrix-Vektor-Multiplikationen. Da A dünn besetzt ist, ist die Berechnung von $y=Ax$ mit geringem Rechenaufwand möglich. Für die Lösung von $c=(LU)^{-1}y$ kann man das äquivalente Gleichungssystem

$LUc=y$

effizient durch Vorwärts-Rückwärts-Einsetzen lösen. Dabei lässt sich die dünne Besetztheit der Matrizen L und U ausnutzen.

Die Berechnung einer ILU-Zerlegung ist, etwa im Vergleich zu Splitting-basierten Vorkonditionierern wie Gauß-Seidel, relativ aufwändig, wobei der Aufwand direkt mit der Größe der erlaubten Besetzungsstruktur zusammenhängt. Aufgewogen wird dies durch die Beschleunigung der Krylow-Unterraum-Verfahren, die in der Regel besser ist, je größer die erlaubte Besetzungsstruktur. Werden direkt hintereinander mehrere Systeme mit derselben Matrix, aber verschiedenen rechten Seiten gelöst, empfiehlt es sich somit, mehr in die Berechnung der ILU zu investieren. Bei der numerischen Lösung zeitabhängiger partieller Differentialgleichungen, bei denen häufig Sequenzen tausender ähnlicher linearer Gleichungssysteme nacheinander zu lösen sind, wird eine einmal berechnete ILU-Zerlegung in der Regel eingefroren und periodisch neuberechnet.

ILU-Zerlegungen oder Varianten sind Teil jeder größeren Programmbibliothek zur Lösung dünnbesetzter linearer Gleichungssysteme, etwa von PetSc oder von MATLAB.

Grundform

In der Grundform wird als Besetzungsstruktur P die von A vorgegeben. Die Zerlegung in die Matrizen L und U wird dann durch folgende Bedingungen definiert:

1. $l_{ij}=u_{ij}=0,{\mbox{falls}}(i,j)\notin P$
2. $(LU)_{ij}=a_{ij},{\mbox{falls}}(i,j)\in P.$

Zusätzlich gilt eine Normierung, d. h. Festlegung über die Hauptdiagonale einer der beiden unvollständig besetzten Dreiecksmatrizen. Dabei werden entweder die Diagonalelemente der unteren Dreiecksmatrix auf 1 normiert:

$l_{ii}=1,\quad i=1,\dotsc ,n$

oder die Diagonalelemente der oberen Dreiecksmatrix:

$u_{ii}=1,\quad i=1,\dotsc ,n$

Je nach Normierung unterscheiden sich die Zerlegungsalgorithmen, was je nach Implementierung auch Auswirkungen auf die Berechnungseffizienz hat.

Da für $(i,j)\notin P$ nicht $(LU)_{ij}=0$ gefordert wird, ist $LU\neq A$ möglich (dies motiviert noch einmal den Namen unvollständige LU-Zerlegung).

Gegeben ist die $n\times n$-Matrix $A=(a_{ij})$. Die unvollständigen Zerlegungsmatrizen L und U werden dann gemeinsam in einer neuen Matrix $M=(m_{ij})$ abgespeichert, wobei die bereits vorher bekannten Einsen auf der Diagonale von L bzw. U nicht gespeichert werden. Die Matrix M wird derart initialisiert, dass sie für Einträge aus P identisch zu A gesetzt wird, andernfalls zu Null.

Bei Wahl der Normierung $l_{ii}=1,\,i=1,\dotsc ,n$ erfolgt die Berechnung der Zerlegung dann mittels des folgenden Algorithmus:

For $k=1,\dotsc ,n-1$, do
   For $i=k+1,\dotsc ,n$ and if $(i,k)\in P$, do
      $m_{ik}:=m_{ik}/m_{kk}$
      For $j=k+1,\dotsc ,n$ and if $(i,j)\in P$, do
         $m_{ij}:=m_{ij}-m_{ik}m_{kj}$

Die Reihenfolge der Schleifen im obigen Algorithmus kann verändert werden, um je nach Datenstruktur die Effizienz zu verbessern. Wird die Matrix beispielsweise zeilenweise abgespeichert, geschehen die Speicherzugriffe in der letzten Schleife nicht auf benachbarte Speicherblöcke. In solchen Fällen ist dann eine Vertauschung von Schleifen sinnvoll.

Existenz

Existenzaussagen der Zerlegung gibt es für M-Matrizen und H-Matrizen. Für allgemeine Matrizen gibt es Gegenbeispiele, bei denen der Algorithmus vorzeitig terminiert, weil eine Null auf der Diagonalen auftaucht, was zu einer Division durch Null führt. Trotzdem ist in der Praxis ein Abbrechen der Berechnung der Zerlegung nicht zu beobachten.

Varianten

Es gibt viele Varianten der ursprünglichen ILU-Zerlegung. Diese versuchen, entweder die Approximationseigenschaften zu verbessern oder bei ähnlicher Approximationsgüte den Berechnungsaufwand zu verkleinern.

ILU(p)

Weit verbreitet sind die ILU(p)-Zerlegungen, die erstmals von Watts 1981 bei der Simulation eines Ölreservoirs betrachtet wurden. Hierbei bezeichnet $p$ den Level of Fill. Die Basisversion der ILU hat den Level 0. Der Level 1 wird dadurch definiert, dass die Besetzungsstruktur des Produkts der Matrizen L und U aus der ILU(0) betrachtet wird. Level 2 ergibt sich aus den Zerlegungsmatrizen von ILU(1) usw. Zur Bestimmung der Besetzungsstruktur einer ILU(p)-Zerlegung ist es nicht nötig, die Zerlegungen der unteren Level vorab zu berechnen. Dazu weist man den Nichtnulleinträgen der Matrix $A$ anfangs den Level 0 zu, den Nulleinträgen dagegen unendlich. Dann durchläuft man die Elemente der Matrix so, wie es im regulären Algorithmus passieren würde und jedes Mal, wenn das Element $a_{ij}$ in der innersten Schleifen modifiziert werden würde, wird der Level aufdatiert mittels

$lev_{ij}=\min(lev_{ij},lev_{ik}+lev_{kj}+1).$

Somit ist es möglich, den Speicher für die ILU-Zerlegung vor Start des Algorithmus bereitzustellen. Bei der Benutzung einer ILU(p) ist zu beachten, dass zum einen die Berechnung der Zerlegung aufwändiger ist als bei der Basisversion und ferner die Anwendung teurer, da der Vorkonditionierer mehr Nichtnulleinträge hat. Damit führen bei hohen Levels etwa ab 3 die Reduktionen der Iterationszahlen im Krylow-Unterraum-Verfahren nicht mehr notwendigerweise zu einer Verkürzung der CPU-Zeiten. Darüber hinaus kann es vor allem bei indefiniten Matrizen sogar zu einer Verschlechterung der Iterationszahlen im Vergleich zur Basisversion kommen.

ILUT

Die ILU(p) haben den Nachteil, dass die Nichtnulleinträge nicht aufgrund von Approximationseigenschaften gewählt werden und somit Rechenzeit für Fast-Nulleinträge vergeudet werden kann. Dies wird in der 1994 von Yousef Saad vorgeschlagenen ILU-Zerlegung mit Threshold, genannt ILUT, berücksichtigt. Hier werden neben dem Einsatz einer Besetzungsstruktur noch zusätzliche Bedingungen zugelassen, nach denen Einträge nicht berücksichtigt werden, falls sie unterhalb einer gewissen Toleranz sind. Etwa für bestimmte diagonaldominante M-Matrizen kann dann wieder die Existenz der Zerlegung bewiesen werden. Die Implementierung einer effizienten ILUT ist schwieriger als bei den anderen Varianten, dafür sind häufig höhere Levels of Fill möglich als bei einer reinen ILU(p).

Fixpunktverfahren

Im Jahr 2015 wurde ein Verfahren vorgeschlagen, welches die ILU-Zerlegung als eine Fixpunktgleichung auffasst. Diese alternative Herangehensweise zeichnet sich durch ihre hohe Parallelisierbarkeit und ihre Einfachheit aus.

Weitere Varianten

Die ILU ist ohne große Probleme auf Blockmatrizen erweiterbar, hierbei muss statt der Division durch das Diagonalelement $d_{ii}$ mit der Inversen des entsprechenden Diagonalblocks multipliziert werden.

Vergleich von ICCG mit CG anhand der 2D-Poisson-Gleichung

Ein Spezialfall ist dagegen die unvollständige Cholesky-Zerlegung (IC). Diese wendet das Konzept der ILU-Zerlegung auf symmetrische und positiv definite Matrizen an, analog zur Cholesky-Zerlegung. Dieses 1977 als erste ILU-Variante eingeführte Verfahren wird häufig als Vorkonditionierer für das CG-Verfahren eingesetzt. Die Kombination CG mit IC wird auch als ICCG bezeichnet.

Parallelisierung

Die klassische ILU-Zerlegung ist streng sequentiell und daher schwer parallelisierbar. Allerdings wurden Varianten entwickelt, die die zentralen Ideen nutzen, um Parallelisierung möglich zu machen. Hierzu gehören insbesondere Multilevel-Techniken wie ILUM. Dabei werden unabhängige Mengen genutzt, um einen Satz Unbekannte blockweise zu eliminieren. Der entstehende Fill-In wird durch einen Threshold begrenzt. Daraufhin wird in den verbleibenden Unbekannten eine neue unabhängige Menge gesucht und der Schritt wiederholt, bis der verbleibende Block klein genug für eine direkte Lösung geworden ist.

Die iterative ILU mittels Fixpunktiteration ist intrinsisch in hohem Maße parallelisierbar.

Einfluss der Nummerierung

Die Nummerierung der Unbekannten in A hat einen nicht zu unterschätzenden Einfluss auf die Effizienz des Vorkonditionierers. Dies liegt daran, dass der Fill-In in der exakten Zerlegung genau davon abhängt und damit die Nummerierung Einfluss darauf hat, wie gut die fehlerbehaftete ILU-Zerlegung A approximiert. Darüber hinaus beeinflusst die Nummerierung die Größe der Einträge auf der Diagonalen und damit die Stabilität.

Auch hier gibt es keine handfesten Aussagen, welche Nummerierung für welche Art von Problemen sinnvoll ist. Insbesondere bei Verwendung der Grundversion ILU(0) sind keine überzeugenden Heuristiken bekannt. Für die stärkeren Vorkonditionierer ILUT oder ILU(p) mit p>0 hat sich in vielen Fällen die Reverse-Cuthill-McKee-Nummerierung als günstig herausgestellt.

Literatur

• Andreas Meister: Numerik linearer Gleichungssysteme, 2. Auflage, Vieweg, Wiesbaden 2005, ISBN 3-528-13135-7
• Gerard Meurant: Computer Solution of Large Linear Systems, Elsevier, 1999, ISBN 0-444-50169-X
• Yousef Saad: Iterative Methods for Sparse Linear Systems, 2nd edition, SIAM Society for Industrial & Applied Mathematics 2003, ISBN 0-89871-534-2