Quanten Monte Carlo

von Ändy

Da soll ja wohl keiner sagen in der Uni lernt man nichts Interessantes. Neulich nämlich war ich schier begeistert von der Magie der Zahlen, und das in einer eher provisorischen Vorlesung bei welcher ich die meiste Zeit damit beschäftigt bin gegen das Einschlafen anzukäpfen und einen interessierten Eindruck zu machen (wir sind da nur zu zweit).

Inhaltlich ging's um oben genanntes Thema, auf das ich aber nicht weiter eingehen will, sondern es nur nenne weil es wie Musik in meinen Ohren klingt. 'Quanten Monte Carlo.' Schön.

Jedenfalls wurde eingehend als kleines Beispiel eine Methode zur Ermittlung der Kreiszahl Pi erläutert, das ich hier kurz wiedergeben will.

Aufgabe: Bestimme den Zahlenwert von Pi.

Das Prinzip besteht darin die Zahl Pi mit einer Wahrscheinlichkeit zu verbinden:


Solch ein Zusammenhang zwischen Pi und P resultiert bspw. aus Zufallsexperimenten. wie etwa dem Schießen auf ein Quadrat mit eingezeichneten Viertelkreis.


Betrachtet wird die Wahrscheinlichkeit eines Treffers in den Viertelkreis. Diese erhält man über das Verhältnis der Flächen:


Bei Kantenlänge 1 für das Quadrat und Kreisradius 1 ergeben sich die Flächen nach:


Und umgestellt nach Pi:


Somit hätten wir einen Bezug für Pi wobei wir nur noch die Wahrscheinlichkeit ermitteln brauchen. Diese ergibt sich bei Durchführung der Experiments aus dem Quotienten der Treffer in den Viertelkreis zur Gesamtschussanzahl:


Die Voraussetzungen für ein kleines Progrämmchen das das Experiment durchführt, die Treffer Zählt und Pi berechnet, wären hiermit geschaffen. Dieses arbeitet mit Zufallszahlen zwischen 0 und 1, die die Koordinaten (x,y) des Schusses angeben. Es wird geprüft, ob der Abstand von der unteren linken Ecke kleiner als 1 ist (Treffer), wobei man den Abstand über die Kreisgleichung oder den Satz von Pythagoras erhält.



Der Output wäre dann dieser:

****     Programm zur Ermittlung der Kreiszahl Pi     ****


Pi_exakt :  3.1415926536


  Schritte |    Pi_exp      |   Differenz    | Fehler in %
-----------+----------------+----------------+------------
   10^ 1   |  3.2000000000  |  0.0584073464  |    1.86E+00
   10^ 2   |  3.1200000000  |  0.0215926536  |    6.87E-01
   10^ 3   |  3.1040000000  |  0.0375926536  |    1.20E+00
   10^ 4   |  3.1188000000  |  0.0227926536  |    7.26E-01
   10^ 5   |  3.1387600000  |  0.0028326536  |    9.02E-02
   10^ 6   |  3.1416360000  |  0.0000433464  |    1.38E-03
   10^ 7   |  3.1423956000  |  0.0008029464  |    2.56E-02
   10^ 8   |  3.1416300400  |  0.0000373864  |    1.19E-03
   10^ 9   |  3.1415979160  |  0.0000052624  |    1.68E-04
   10^10   |  3.1416006086  |  0.0000079551  |    2.53E-04

Je mehr Schüsse ausgewertet werden, umso sicherer ist die Aussage über Pi. Aufgrund von irgendwelchen numerischen Effekten oder Unzulänglichkeiten des Zufallsgenerators ist die Genauigkeit aber beschränkt, was man an dem größeren Fehler bei 10^10 gegenüber 10^9 sieht. Auch lässt sich die Schrittanzahl nicht über 10^10 hinaus steigern. Die Klärung dieses Sachverhalts steht noch aus.
Unabhängig von diesen störenden Einflüssen ist der Fehler allgemein proportional zur Wurzel aus der Schussanzahl, also zur Halbierung des Fehlers müssen 4mal soviel Schüsse abgefeuert werden. Die Konvergenz des Verfahrens ist also schlecht, weshalb heute andere Verfahren zur Ermittlung von Pi benutzt werden.
Das Vorgehen zur Bestimmung einer Unbekannten über ein Zufallsexperiment nennt man allgemein Monte Carlo Verfahren. Wie gesehen ist derartiges gut zur Orientierung, aber nicht zur Exakten Bestimmung.

Damit wäre ich am Ende mit dem kleinen Exkurs in das Reich der Zahlen und der Geometrie.

---

Was dem Camper sein Schlauchboot, das dem Wissenschaftler sein Hochleistungsrechner. Zur Bestärkung meines Egos habe ich es mir deshalb zur Aufgabe gemacht das geschriebene Programm auf mehreren Rechnern gleichzeitig laufen zu lassen um so schneller noch mehr Schüsse abzufeuern was die Genauigkeit des ermittelten Pi erhöht. Dabei möchte ich das Programm aber nicht auf jedem Rechner separat starten, sondern von einem zentralen Ort aus steuern. Dies geht über eine Erweiterung meiner Programmiersprache Fortran mittels OpenMP. Dieses Interface verteilt bei entsprechender Anweisung meinen "Job" auf mehrere Rechner, und sammelt am Ende mein Ergebnis ein. Wenn ich also z.B. 6mal schießen möchte und 2 CPU's habe, so bekommt jeder Prozessor 3 Schüsse zur Aufgabe. Die Trefferanzahl jedes Individuums wird am Ende an den Führer-CPU geschickt, der mir dann mit der Gesamttrefferanzahl Pi berechnet. Dadurch bin ich doppelt so schnell.


Zunächst habe ich aber noch das Ausgangsprogramm perfektioniert, so dass der oben beschriebene Effekt des größeren Fehlers bei 10^10 Schüssen nicht mehr auftritt, denn ich möchte ja möglichst noch mehr Schüsse machen: mc2.f90.
Für 10^10 Schüsse brauch ich jetzt 11min13sec.

Um das Programm nun auf mehrere Rechner zu Verteilen muss ich mir noch einen neuen Zufallsgenerator zulegen, da der bisherige nicht skaliert, wie man so schön sagt, was soviel heißt, dass ich bei Aufteilung des Programms keine Zeit gut mache. Dies liegt an der Zentralen Speicherung der Zufallszahlen. Wenn ein Prozessor dann eine haben will, "so muss er zunächst zur Zentrale laufen sich da eine holen, ggf. mit Anstehen." Besser ist es, wenn jeder seine eigenen Zahlen hat. Deshalb setzt ich den Ziggurat-Algorithmus ein.
Ein Rechner benötigt für 10^10 jetzt 8min29sec (seq:CODE).
Die 16 Rechner unseres Rechenclusters in der Chemie benötigen für das Gleiche nur 32sec (omp:CODE). Geil!!!

Für 10^11 Schüsse benötige ich mit 16 Rechnern 5min19sec statt 84min54sec. Ein Output hier:

****     Programm zur Ermittlung der Kreiszahl Pi     ****


 Pi_exakt :  3.1415927410


   Schritte |    Pi_exp      |   Differenz    | Fehler in %
 -----------+----------------+----------------+------------
    10^ 1   |  3.5999999046  |  0.4584071636  |    1.46E+01
    10^ 2   |  3.4800000191  |  0.3384072781  |    1.08E+01
    10^ 3   |  3.1800000668  |  0.0384073257  |    1.22E+00
    10^ 4   |  3.1487998962  |  0.0072071552  |    2.29E-01
    10^ 5   |  3.1285600662  |  0.0130326748  |    4.15E-01
    10^ 6   |  3.1405038834  |  0.0010888577  |    3.47E-02
    10^ 7   |  3.1414597034  |  0.0001330376  |    4.23E-03
    10^ 8   |  3.1414175034  |  0.0001752377  |    5.58E-03
    10^ 9   |  3.1416726112  |  0.0000798702  |    2.54E-03
    10^10   |  3.1416304111  |  0.0000376701  |    1.20E-03
    10^11   |  3.1415936947  |  0.0000009537  |    3.04E-05

Auch habe ich die Fehler aller Programme nochmal verglichen um zu sehen ob beide Zufallszahlgeneratoren gleich gut arbeiten und ob 16 Rechner genauso gut arbeiten wie einer, um Programmierfehler meinerseits auszuschließen:


Ihr seht, dass alles funktioniert und sich der Fehler in der tat wie 1/N^1/2 verhält (c=1.25).

Erwähnenswert ist diese Parallelisierungsstrategie besonders weil das Arbeiten mit mehreren Prozessoren gleichzeitig in Zukunft nicht nur Exoten vorbehalten sein wird, sondern man bereits jetzt mit Dual- und Quad-Cores Rechner unter seinem Schreibtisch stehen hat, die das auch können und bei der Entwicklung hin zu Octal- und Hexdecal-Cores das Prinzip unverzichtbar für eine effektive Ausnutzung des Rechenknechts wird.

So, bestärkt in meiner Hochschätzung mir selbst gegenüber geht's jetzt aber nicht mehr um Zahlen, sondern gut gelaunt in die Nacht hinein.