Forschungsbericht 2008 - Max-Planck-Institut für Mathematik in den Naturwissenschaften

Schnelle Löser und der Fluch der Dimension

Autoren
Grasedyck, Lars
Abteilungen
Zusammenfassung
In Anwendungen treten oft große (lineare) Gleichungssysteme auf. Diese schnell zu lösen ist das Ziel der „Schnellen Löser“. Anfänglich wurden nur schwach besetzte Systeme betrachtet, aber Dank der am MPI-MIS entwickelten „Hierarchischen Matrizen“ ist dies für eine wesentlich umfangreichere Klasse von Matrizen möglich. Diese Technik stößt an ihre Grenzen, wenn statt der üblichen dreidimensionalen Probleme höhere Dimensionen auftreten – eine Herausforderung an moderne mathematische Methoden, der sich die Gruppe „Wissenschaftliches Rechnen“ am MPI-MIS stellt.

Entstehung von großen Gleichungssystemen

Die Behandlung von physikalischen, chemischen, biologischen oder wirtschaftlichen Prozessen führt häufig zu Gleichungssystemen, die sehr groß werden. Diese können dann nach einer Linearisierung im Prinzip gelöst werden:

Einfache Auflösung per Hand

Löst man ein Gleichungssystem mit 2 Unbekannten und 2 Gleichungen, so kommt man z.B. durch Einsetzen schnell zum Ergebnis. Bei 3 Gleichungen und 3 Unbekannten wird dies etwas umständlicher. Bei 10 Gleichungen und 10 Unbekannten braucht man schon eine ausgefeilte Technik, etwa die Gauß-Elimination, um ein solches System per Hand zu lösen – und auch dann kann es durchaus eine Stunde dauern. Bei mehreren hundert Unbekannten und Gleichungen wird man den Prozess automatisieren wollen und dem Computer die lästige Arbeit überlassen.

Automatisierte Auflösung per Computer

Doch selbst moderne Hochleistungsrechner stoßen an ihre Grenzen, wenn die Anzahl der Gleichungen in die Millionen geht. Der Grund dafür ist, dass der Aufwand der Gauß-Elimination kubisch in der Anzahl N der Unbekannten und Gleichungen ist, d.h. bei einem zehnmal größeren Gleichungssystem braucht man 1000- mal länger. Wenn also zum Beispiel zur Simulation elektromagnetischer Signale im Gehirn eines Menschen bei einer Genauigkeit von einem Zentimeter innerhalb von

1 Sekunde

die auftretenden Systeme gelöst werden können und nun zur sicheren Planung einer Operation eine Genauigkeit von einem Millimeter gewünscht wird, so ergibt sich ein 1000-mal so großes System dessen Auflösung dann eine Milliarde mal länger dauert. Das wären

32 Jahre.

Auch die geduldigsten Patienten empfinden dies als eine zu lange Wartezeit. Es lohnt sich daher, erst einmal einige Jahre zu überlegen, ob dies nicht mit einem geringeren Aufwand geht. Mit dieser Frage beschäftigt sich der Themenkreis „Schnelle Löser".

Schnelle Löser

Moderne und effiziente Verfahren zur Behandlung großer linearer Gleichungssysteme zeichnen sich dadurch aus, dass ihr Aufwand linear in der Anzahl der Unbekannten skaliert, also ein zehnmal so großes System in der zehnfachen Zeit aufgelöst werden kann. In unserem obigen Beispiel wären es dann eben nur

20 Minuten

zur Berechnung der gewünschten genauen Lösung. Mit dieser Wartezeit könnte ein Patient (im wahrsten Sinne des Wortes) gut leben.

Schnelle Löser für große schwachbesetzte Gleichungssysteme

Wie erreicht man die lineare Komplexität? Diese Frage ist nicht ganz so einfach zu beantworten. Auf jeden Fall muss die zugrunde liegende Matrix eine spezielle Struktur haben, denn ansonsten wäre schon das Auslesen der Matrixeinträge (dies sind bei N Unbekannten und N Gleichungen N ∙ N = N2) von quadratischer Komplexität. In unserem Beispiel von oben sind das immerhin noch gute

11 Tage.

Ein erster, sehr verbreiteter und sehr erfolgreicher Ansatz besteht darin anzunehmen, dass die Matrix nur eine Anzahl von Nichtnull-Einträgen proportional zur Anzahl der Unbekannten N besitzt.

Diese so genannten schwachbesetzten Matrizen (oben drei verschieden große Matrizen mit wenigen (roten) Nichtnull-Einträgen) können sehr effizient gespeichert und ausgewertet werden. Das alleine reicht zwar noch nicht, um das System auch mit linearer Komplexität zu lösen, aber es erscheint dann plausibel.

Schnelle Löser für große vollbesetzte Gleichungssysteme

Die Situation erscheint hingegen aussichtslos, wenn die Matrix keine oder nur unwesentlich wenige Null-Einträge besitzt. Wie kann es dann möglich sein, das System schnell zu lösen? Die Antwort wurde bereits gegeben: Das System, also die Matrix, muss strukturiert sein, das heißt sie besitzt intrinsisch wesentlich weniger als N2-Informationen. Solche Matrizen, die mithilfe von wenigen Daten beschrieben werden können, nennt man datenschwach. Diesen Effekt kennt man von Bildern, die unkomprimiert sehr viel Speicher verbrauchen, bei geeigneter Kompression hingegen sehr viel weniger Speicherplatz benötigen. (Die nachfolgende Serie von Bildern zeigt die Ausrichtung von Leitfähigkeitstensoren im menschlichen Gehirn, die Größe der Bilder von links nach rechts ist 73 KB, 23 KB, 12KB.)

Bei Bildern sucht man nach Bereichen, in denen sich wenig verändert. Einfache Kompressionen vergleichen dabei nur Farbwerte, moderne Kompressionstechniken hingegen setzen beispielsweise auf Wavelet-Approximationen. Hierbei wird das Bild nicht mehr exakt gespeichert, sondern nur eine sehr gute Approximation, die für das Auge nicht mehr vom Original zu unterscheiden ist. Diese Idee überträgt sich auch auf Matrizen: Wenn bei der Modellierung und Linearisierung eines realen Prozesses ohnehin ein gewisser (kleiner) Fehler im Zuge der Abbildung in einem Computer gemacht wurde, dann muss das System nicht mehr exakt dargestellt und gelöst werden, sondern nur bis auf einen (kleinen) Approximationsfehler. Dieses Schlupfloch kann man benutzen, um eine vollbesetzte Matrix effizient zu speichern. Die Technik der hierarchischen Matrizen (H-Matrizen) basiert darauf, die Matrix in Blöcke unterschiedlicher Größe zu unterteilen und jeden Matrixblock approximativ mit einem niedrigen Rang, d.h. einer mathematisch gut fassbaren Struktur, darzustellen. (Nachfolgend sind zwei Beispiele solcher Unterteilungsmuster abgebildet.)

Dadurch wird der Speicherbedarf – natürlich abhängig von der zugrunde liegenden praktischen Anwendung – gesenkt und es ist möglich, ein solches System auch noch mit linearem Aufwand zu lösen [1-4]. Dies ist ein ganz fundamentaler Unterschied zu dem Bild-Beispiel: Ein komprimiertes Bild wird zur Weiterverarbeitung erst wieder dekomprimiert. Das wäre für eine sehr große Matrix ausgesprochen unpraktisch. Mit einer Matrix im komprimierten H-Matrix-Format kann man hingegen wie mit einer normalen Matrix rechnen, ohne sie dekomprimieren zu müssen.

Ein Beispiel aus der Quellenlokalisation im menschlichen Gehirn

Simulation elektrischer Felder

Ein physikalisch und mathematisch sehr gut untersuchtes Problem ist die Simulation elektrischer und elektromagnetischer Felder, die zum Beispiel in sehr schwacher Form im Gehirn auftreten. In der klinischen Praxis ist die direkte Simulation der Felder nicht von primärem Interesse, sehr wohl aber die folgende Frage: Wenn elektrische und magnetische Felder außerhalb des Kopfes mit hoher Genauigkeit gemessen werden, kann man dann rückschließen auf die Region des Gehirns, in der die Quelle der Aktivität liegt? Dies ist deshalb so interessant, weil es einerseits dank moderner Technik möglich ist, sehr genau zu messen, und der Patient sich dafür auch nicht einer Operation unterziehen muss oder großen Strahlungen ausgesetzt wird. Andererseits lassen sich mit der gewonnenen Information Abläufe im Gehirn gut studieren. Mathematisch spricht man von einem „Inversen Problem", da von den gemessenen Daten zurück auf die Quellen geschlossen wird. Zur Lösung des Inversen Problems werden zahlreiche Simulationen von elektrischen Feldern benötigt, welche wiederum auf große Gleichungssysteme führen [5, 6].

Anwendung in der Behandlung von Epilepsie-Patienten

Die Lokalisation von Hirnaktivität spielt eine Rolle bei der so genannten prächirurgischen Diagnostik, also der Untersuchung eines Patienten zur Vorbereitung einer Operation. Hier ist ganz offensichtlich eine hohe Präzision gefragt, gleichzeitig ist aber auch die zur Verfügung stehende Zeit vor der Operation eingeschränkt, sodass schnelle Verfahren notwendig sind.

(Obige Grafik zeigt die Lokalisation der von Epilepsie betroffenen Region im Gehirn eines Patienten. Das Bild wurde großzügig von Dr. Carsten Wolters bereitgestellt.)

Für Epilepsie-Patienten geben die Daten während eines epileptischen Anfalls Aufschluss über Lage und Ausbreitung des verantwortlichen Bereichs im Gehirn. Dieser Bereich wird, sofern andere Behandlungen nicht erfolgreich oder nicht anwendbar sind, operativ entfernt. Man möchte ihn daher sehr genau bestimmen, um weitere Anfälle ausschließen zu können und um die Menge des zu entfernenden Gewebes zu minimieren.

Die nächste Dimension – der Fluch der Dimension

In den bisherigen Betrachtungen sind stets deterministische Probleme zu lösen gewesen, also Probleme, bei denen sämtliche notwendige Daten exakt und vollständig gegeben sind. Dies ist in der Praxis allerdings nicht immer der Fall. Stattdessen sind die Daten mit Unsicherheiten behaftet, ungenau oder in einigen Fällen auch schlicht nicht bekannt. Dies ist zum Beispiel bei lebendigen Patienten der Fall, deren Gehirn ununterbrochen aktiv ist und damit auch schwache elektrische Felder erzeugt, deren Quellen nicht bekannt sind.

In diesem Fall ist man an einer Wahrscheinlichkeitsverteilung der Lage der Quellen interessiert, oder der Varianz der Quellorte. Die mathematisch präzise Modellierung der Wahrscheinlichkeit erhöht die Dimension des Problems, d.h. anstelle einer dreidimensionalen Simulation (das Gehirn) führt man vier- fünf oder noch viel höherdimensionale Simulationen durch. Man kann sich die zusätzlichen Dimensionen wie Parallelwelten vorstellen, in denen jeweils ein Gehirn deterministisch gerechnet wird. Da wir aber vorher nicht wissen, welche der Rechnungen dem tatsächlichen Patienten entspricht, müssen wir alle Rechnungen durchführen.

Wird nun zur Verbesserung der Genauigkeit die Auflösung von einem Zentimeter auf einen Millimeter vergrößert, so vervielfacht sich die Problemgröße für jede auftretende Dimension um 10. Bei einem zehndimensionalen Problem (3 für die Raumkoordinaten und 7 für den Wahrscheinlichkeitsraum) würde sich die Problemgröße um

10000000000 = 10 Milliarden

erhöhen. Mit immer höherer Dimension vergrößert sich dieser Faktor und wird schnell größer als die Anzahl der Atome in dem uns bekannten Universum. Man spricht daher vom „Fluch der Dimension". Ein aktuelles Feld der Mathematik von großem Interesse ist die effiziente Behandlung von hochdimensionalen Problemen, also Strategien, welche den Fluch der Dimension brechen können.

Beteiligte Projekte

Weitere Informationen zu „Hierarchischen Matrizen" findet man auf der Seite

http://www.hmatrix.org.

Die Modellierung und Lösung des EEG/MEG Inversen Problems ist Teil des DFG-Projektes „Entwicklung und Validierung von Verfahren zur Lokalisation von Hirnaktivität mithilfe der Methode der Finiten Elemente“ (1.12.2006 bis 30.11.2009) mit den Antragstellern Dr. Carsten H. Wolters (Westfälische Wilhelms-Universität Münster), Prof. Dr. Lars Grasedyck (MPI-MIS und Berlin Mathematical School an der TU-Berlin), Prof. Dr. Jens Haueisen (Friedrich-Schiller-Universität Jena) und Dr. Thomas R. Knösche (Max-Planck-Institut für Kognitions- und Neurowissenschaften in Leipzig).

Die hochdimensionalen Probleme sind der Fokus des 2008 startenden DFG-Schwerpunktprogrammes 1324 „Extraktion quantifizierbarer Information aus komplexen Systemen", siehe

http://www.dfg-spp1324.de/

Originalveröffentlichungen

W. Hackbusch:
A sparse matrix arithmetic based on H-matrices. Part I: Introduction to H-matrices.
Computing 62, 89-108 (1999).
L. Grasedyck, W. Hackbusch:
Construction and arithmetics of H-matrices.
Computing 70, 295-334 (2003).
S. Börm, L. Grasedyck, W. Hackbusch:
Introduction to hierarchical matrices with applications.
Engineering Analysis with Boundary Elements 27, 405-422 (2003).
S. Börm:
H2-matrix arithmetics in linear complexity.
Computing 77(1), 1-28 (2006).
C. H. Wolters, L. Grasedyck, W. Hackbusch:
Efficient computation of lead field bases and influence matrix for the FEM-based EEG and MEG inverse problem.
Inverse Problems 20(4), 1099-1116 (2004).
C. H. Wolters, H. Köstler, C. Möller, J. Härdtlein, L. Grasedyck, W. Hackbusch:
Numerical mathematics of the subtraction method for the modeling of a current dipole in EEG source reconstruction using finite element head models.
SIAM Journal on Scientific Computing 30(1), 24-45 (2007).
Zur Redakteursansicht