3. Schema gesucht: Sehr teure Gleichungen parallel lösen, ohne viel Rechenzeit zu verschwenden

Schema gesucht: Sehr teure Gleichungen parallel lösen, ohne viel Rechenzeit zu verschwenden

  • Mod

    Achtung, jetzt wird es erst einmal mathematisch (aber nicht schwer! Bloß abstrakt), ich suche aber einen Algorithmus, daher frage ich hier:

    Ich habe eine Funktion E(x)E(x). Die Funktionswerte sind sehr sehr teuer zu berechnen, sagen wir ein CPU-Tag, nicht parallelisierbar. Aber ich kann mehrere E(x)E(x) gleichzeitig parallel berechnen.

    Nun habe ich eine weitere Funktion
    F(E(x_1),E(x_2)) mit x_1<x_2F(E(x\_1),E(x\_2)) \text{ mit } x\_1 < x\_2,
    die auf diesen E(x)E(x) operiert. Diese Funktion ist billig. Es gilt (falls es hilft):
    F(E(x),E(x+δx_1))<F(E(x),E(x+δx_2)),dann und genau dann, wenn δx_1<δx_2F(E(x), E(x+\delta x\_1)) < F(E(x), E(x+\delta x\_2)), \text{dann und genau dann, wenn } \delta x\_1 < \delta x\_2
    In Worten: F ist monoton, wenn man eines der Argumente festhält. (Da gibt's bestimmt auch noch Fachsprech für, das mir gerade nicht einfällt)
    Außerdem ist F stetig, (numerisch) differenzierbar und alles was man sich sonst noch so wünschen kann.

    Problemstellung:
    Nun möchte ich eine Kette von x'en (x\_0 = x\_{start}, x\_1, x\_2, \dots , x_{end}) mit vorgegebenem xstartx_{start} und einem ungefähren Wunschwert xendx_{end} erzeugen. Als Nebenbedingung soll dabei gelten
    F(Ex_i,E_xi+1)=CF(E_{x\_i}, E\_{x_{i+1}}) = C
    wobei CC eine vorgegebene Konstante ist. Wegen der oben angegebenen Eigenschaft der Funktion F ist es immer möglich, zu einem gegebenen xix_i ein (eindeutiges) xi+1x_{i+1} zu finden, für das diese Bedingung exakt erfüllt ist. Aber eine Approximation reicht mir! Sie muss noch nicht einmal besonders gut sein. Wie genau formuliert werden soll, wann eine Approximation ausreichend ist, lasse ich erst einmal offen. Meinetwegen kann man es so formulieren, dass ich eine vorgegebene Abweichung δC\delta C toleriere.

    Bedingungen an den Algorithmus:
    Ich möchte nun möglichst effizient herausfinden, welchen Wert diese xix_i ungefähr haben müssen. Dabei möchte ich von dem Algorithmus zwei Dinge:
    1. Ich will nicht ewig auf das Ergebnis warten. Ich erwarte, dass die Kette 10-50 Werte enthalten wird.Die Angabe oben mit dem einen Tag pro E(xi)E(x_i) war nicht übertrieben, eher untertrieben. Ich sag mal, ein bis zwei Wochen kann ich mich gedulden.
    2. Ich habe nicht unendlich Ressourcen. Ja, ich kann Rechenzeit verschwenden, es muss nicht jeder berechnete Wert am Ende im Ergebnis auftauchen. Aber ich kann nicht einfach das Zielintervall in 10.000 Schritte aufteilen, zu jedem Schritt das EE ausrechnen und mir dann die passenden Werte rauspicken. Es muss schon irgendwie halbwegs gezielt passieren. Aber ein paar Dutzend Schüsse ins Blaue sind noch ok. Sagen wir ein paar CPU-Monate sollten insgesamt nicht überschritten werden.

    So, nun suche ich Ideen, wie man das möglichst effizient machen könnte. Bitte keine ganz naiven Vorschläge, die euch beim ersten Querlesen gekommen sind. Ich habe schon eine ganze Weile darüber nachgedacht und das Problem ist nicht einfach. Ich habe aber auch nie eine informatische Algorithmenvorlesung besucht, also wenn euch die Problemstellung bekannt vorkommt, dann nur raus damit.

    Was ich bisher gemacht habe:
    Es gibt hier zwei/drei Probleme:
    1. Eine iterative Lösung ist halbwegs schwer, da man dafür normalerweise die Ableitung von F berechnen müsste. Und dazu braucht man mehrere E. Und dann wird's teuer. Außerdem kann die Ableitung vom F durchaus wild sein, eine Approximation durch zwei weit auseinanderliegende Werte könnte schlecht sein. Aber wie schon gesagt, so schrecklich genau braucht es nicht zu sein.
    2. Die Kettenbildung scheint mir inhärent seriell. Ich muss bei xstartx_{start} anfangen und mich dann von x zu x weiterhangeln. Wenn ich zwischendurch E's berechne, dann kann ich mir diese zwar merken und später vielleicht nochmal nutzen. Aber insgesamt scheint mir dadurch doch die Zeitbeschränkung schwer. Die E's die ich berechne sollten also möglichst gut angelegt sein.
    (3. Metaproblem: Solche Arten von Problemen algorithmisch zu lösen habe ich nie gelernt. Hier ist man mal ganz praktisch durch Hardware und Zeit beschränkt. Ich brauche eine Lösung die in der Praxis gut funktioniert, nicht in meinem Theoretikerturm, indem ich mich bisher so wohl gefühlt habe)

    Naiver Algorithmus, den ich mir ausgedacht habe:
    Starte bei xix_i. Gucke schon berechnete E(x>xi)E(x > x_i) an und bestimme so ein Intevall, in dem xi+1x_{i+1} liegen muss. Dann Intervallschachtelung. Oder vielleicht besser ein Newtonverfahren? Ist halt nichtlinear 😞 . Zwischendurch berechnete E's natürlich aufheben.

    Nachteil: Ich nutze nicht die Möglichkeit, dass ich mehrere E's parallel berechnen kann, da dieser Algorithmus seriell ist.

    Mögliche Verbesserung: Es könnte was bringen, das Zielintervall (sowohl das konkrete Intervall zum Berechnen des nächsten xi+1x_{i+1}, als auch das noch übrige Gesamtintervall, in dem xi+2x_{i+2} und folgende liegen werden) im Vorhinein schon grob einzuteilen und einfach auf gut Glück ein paar E's darin verteilt zu berechnen. Das lässt aber zwei Fragen offen:
    1. Ist das gut angelegte Rechenzeit? Kann ich leider gar nicht einschätzen.
    2. Wie diese Werte verteilen? Ich weiß a priori gar nichts über das F, außer den oben angegebenen Bedingungen. Verteilt man die spekulativen Werte dann gleichmäßig auf dem Intervall? Oder gibt es da was besseres?

    0
  • C

    Ich kann dir wahrscheinlich nicht wirklich weiterhelfen, einfach weil mir der mathematische Hintergrund fehlt. Mich würde aber E(x) interessieren. Kannst (darfst) du die Funktion posten?

    0
  • Mod

    cooky451 schrieb:

    Ich kann dir wahrscheinlich nicht wirklich weiterhelfen, einfach weil mir der mathematische Hintergrund fehlt. Mich würde aber E(x) interessieren. Kannst (darfst) du die Funktion posten?

    Ich könnte jetzt zwar einen Klugscheißerausdruck hinschreiben, aber das was du wissen möchtest: Die Funktionswerte sind Ergebnis einer Simulation, das x ist ein Parameter dieser Funktion.

    0
  • S

    Kann man bzgl. E gar keine Annahmen treffen? Die Funktion springt wild umher, ist weder approximierbar noch sonstwie monoton, stetig, etc.?

    Wenn du nämlich "F(E(x), E(x+1)) = C +- deltaC" haben willst, und F monoton wachsend ist, dann kannst du ja ein ungefähres Wachstum von F angeben, oder? Bspw. ob F exponentiell steigt, etc.

    Wenn man nun über E besser Bescheid wüsste, könnte man ja sicher eine gute Heuristik zum Raten bauen, oder?

    MfG SideWinder

    0
  • Mod

    Das E ist stetig, es ist differenzierbar und sonstwie gutartig.

    Das E ist sicherlich auch nicht wild, sprich: Die erste und zweite Ableitung ist eher klein. Man kann erwarten, dass es im Gesamtintervall vielleicht ein bis zwei "Schwingungen" macht. (kein ganz zutreffendes Wort, da es keine reelle Zahl ist, aber das trifft's ganz gut)
    Wie könnte mir das helfen? Heuristik klingt gut. Man könnte die auch ein bisschen anfüttern, indem man vorher das Gesamtintervall grob absteckt.

    Was dem jedoch zuwieder läuft: In das F fließen immer zwei E's ein, es zählt eher so etwas wie die Differenz zweier Werte (ist leider etwas komplizierter). Daher nützt mir der absolute Wert des E nicht so viel, es zählt eher das Verhalten. Könnte jedoch gehen, indem man die erste Ableitung abschätzt.

    0
  • K

    F(Ex_i,E_xi+1)−C=0F(E_{x\_i}, E\_{x_{i+1}}) - C = 0 Damit reduziert es sich auf das Finden von Nullstellen, da xix_i bzw. xstartx_{start} vorgegen ist. Wenn F, wie du sagst, billig ist, kann der Wer von E(xi+1)E(x_{i+1}) berechnet werden. Dann bleibt 0=E(x)−E(xi+1)0 = E(x) - E(x_{i+1}) wieder als Nullstellenproblem uebrig. ... so als Anzatz.

    0
  • Mod

    knivil schrieb:

    F(Ex_i,E_xi+1)−C=0F(E_{x\_i}, E\_{x_{i+1}}) - C = 0 Damit reduziert es sich auf das Finden von Nullstellen, da xix_i bzw. xstartx_{start} vorgegen ist. Wenn F, wie du sagst, billig ist, kann der Wer von E(xi+1)E(x_{i+1}) berechnet werden. Dann bleibt 0=E(x)−E(xi+1)0 = E(x) - E(x_{i+1}) wieder als Nullstellenproblem uebrig. ... so als Anzatz.

    Kapier ich gerade irgendwie nicht.

    Wenn F, wie du sagst, billig ist, kann der Wer von E(xi+1)E(x_{i+1}) berechnet werden.

    Jain. Ein nötiger Wert kann berechnet werden. Es gibt aber mehrere Möglichkeiten dafür. Ich weiß aber nur, dass konstruktionsbedingt eine dieser Möglichkeiten auftreten wird. Aber nehmen wir mal an, ich kenne das nötige E(xi+1)E(x_{i+1}), dann will ich ja wissen, was xi+1x_{i+1} ist:

    Dann bleibt 0=E(x)−E(xi+1)0 = E(x) - E(x_{i+1}) wieder als Nullstellenproblem uebrig.

    Ja. Aber das zu lösen ist doch gerade mein Problem, bloß umformuliert.

    0
  • M

    Seien x_1<x_2<x3x\_1 < x\_2 < x_3. Wenn man F(E(x_1),E(x_2))F(E(x\_1), E(x\_2)) und F(E(x_1),E(x_3))F(E(x\_1), E(x\_3)) hat, welche Aussagen kann man dann für F(E(x_2),E(x_3))F(E(x\_2), E(x\_3)) machen?

    0
  • K

    SeppJ schrieb:

    Ja. Aber das zu lösen ist doch gerade mein Problem, bloß umformuliert.

    Ja, aber es steht jetzt explizit als Nullstellenproblem da. Nunn kannst du beispielsweise Intervallschachtelung anwenden und parallelisierst diese. Normalerweise berechnest du den naechsten Punkt E(x') und entscheidest dann, ob mit E(x+dx) oder E(x-dx) fortgefahren wird. E(x'), E(x+dx) und E(x-dx) kannst du parallel berechnen. PS: Das ist ein sehr naiver Ansatz.

    0
  • A

    http://ars.userfriendly.org/cartoons/?id=19980708

    0
  • R

    ohne deine funktion zu kennen ist es nicht ganz trivial etwas nicht-dummes vorzuschlagen, finde ich ;), da du dein problem ja schon sehr genau kennst und vermutlich alle ideen die einem beim durchlesen in den kopf kommen auch schon hattest.

    deine beschreibung ist auch ein wenig verwirrend. da du ueber stetigkeit/differenzierbarkeit sprichst, suchst du naehrungsverfahren die besser sind als z.B. newton rhapson? da du von parallelisierung der funktion sprichst, suchst du eine einfach loesung um sie mehrmals parralel fuer verschiedene aufgabenstellungen abzuarbeiten? sowas wie openCL oder http://ispc.github.com/ ?
    zudem sagst du, dass du das verfahren fuer einen wert nicht parallelisieren kannst, dann schreibst du dass du keine 10k schritte machen willst (und zudem ist die funktion stetig), somit koenntest du schon mehrere punkte gleichzeitig fuer ein problem berechnen und dadurch die naehrung verbessern? (z.b. newton-rhapson mit einer spline statt nur der funktions-tangente, wenn du z.b. immer 4 samples berechnen wuerdest).

    suchst du komplett neue loesungsansaetze oder optimierungen fuer was du schon hast?

    0

  • Hast du mal darueber nachgedacht E(x) auf einer GPU zu implementieren?

    0
  • B

    Als Nebenbedingung soll dabei gelten F(Exi,Exi+1)=CF(Exi,Exi+1)=C wobei C eine vorgegebene Konstante ist.

    Diese Bedingung ist irgentwie ein Knackpunkt.

    Wie wäre es wenn du E(x) durch eine einfache Funktion E'(x) approximierst und danach erst die Funktion F(x, y) mit E'(x) bestimmst ?

    Du musst du wohl erst einige E(xi)'s berechen und mit denen eine passende Approximation für E(x) suchen. Dabei würde ich versuchen den Grundtyp von E'(x) festzustellen. Nicht das E(x) vom Typ her eine Schwingung ist und du als Ansatz für E'(x) ein Polynom nimmst. Ist vermutlich eine Optimierungsaufgabe über Approximationen, abhängig von dem simulierten System.

    0
  • K

    Wie waere es mit einer Taylorentwicklung ...

    0
  • Mod

    Sehr viel dazugekommen. eine lange Antwort:

    rapso schrieb:

    deine beschreibung ist auch ein wenig verwirrend. da du ueber stetigkeit/differenzierbarkeit sprichst, suchst du naehrungsverfahren die besser sind als z.B. newton rhapson?

    Die Stetigkeit und Differenzierbarkeit ist nur Nebeninformation, falls es etwas nützt. Ich kann leider nicht sehr viel über die Funktion aussagen.

    da du von parallelisierung der funktion sprichst, suchst du eine einfach loesung um sie mehrmals parralel fuer verschiedene aufgabenstellungen abzuarbeiten?

    Das kann ich schon. Ich kann viele Funktionswerte gleichzeitig bestimmen, aber die Berechnung des jeweiligen Einzelwertes ist nicht weiter verbesserbar.

    zudem sagst du, dass du das verfahren fuer einen wert nicht parallelisieren kannst, dann schreibst du dass du keine 10k schritte machen willst (und zudem ist die funktion stetig), somit koenntest du schon mehrere punkte gleichzeitig fuer ein problem berechnen und dadurch die naehrung verbessern? (z.b. newton-rhapson mit einer spline statt nur der funktions-tangente, wenn du z.b. immer 4 samples berechnen wuerdest).

    Ja, das ist möglich und etwas, woran ich noch nicht gedacht hatte. (Erstaunlicherweise muss ich sagen. Jetzt wo du es sagst ist es dumm, dass ich daran nicht gedacht habe)

    suchst du komplett neue loesungsansaetze oder optimierungen fuer was du schon hast?

    Jain. Dieses Problem ist neu für mich. Ich habe noch nichts funktionierendes, da ich bei allen meinen bisherigen Ansätzen zu dem Schluss gekommen bin, dass sie zu lange dauern würden.

    TGGC schrieb:

    Hast du mal darueber nachgedacht E(x) auf einer GPU zu implementieren?

    Kommt drauf an: Zur Berechnung eines Wertes würde es nichts bringen, es würde sogar wesentlich langsamer, da der einzelne GPU-Kern nicht mit einer CPU mithalten kann. Es könnte was bringen, um parallel mehrere Werte zu berechnen. Aber der Portierungsaufwand wäre enorm und ich kann bereits dutzende Einzelwerte parallel berechnen (aber nicht hunderte, wie bei einer GPU). Und die vorhandene Hardware ist schon auf ein klassisches CPU-lastiges Programm ausgelegt. Aus diesen Gründen habe ich die Idee verworfen.

    knivil schrieb:

    Wie waere es mit einer Taylorentwicklung ...

    Taylorentwicklung ist etwas für theoretische Mathematik, nicht für angewandte Numerik. Oder möchtest du mir etwas konkretes mit dem Stichwort sagen, dass ich gerade nicht verstehe?

    Bitte ein Bit schrieb:

    Als Nebenbedingung soll dabei gelten F(Exi,Exi+1)=CF(Exi,Exi+1)=C wobei C eine vorgegebene Konstante ist.

    Diese Bedingung ist irgentwie ein Knackpunkt.

    Wie wäre es wenn du E(x) durch eine einfache Funktion E'(x) approximierst und danach erst die Funktion F(x, y) mit E'(x) bestimmst ?

    Du musst du wohl erst einige E(xi)'s berechen und mit denen eine passende Approximation für E(x) suchen. Dabei würde ich versuchen den Grundtyp von E'(x) festzustellen. Nicht das E(x) vom Typ her eine Schwingung ist und du als Ansatz für E'(x) ein Polynom nimmst. Ist vermutlich eine Optimierungsaufgabe über Approximationen, abhängig von dem simulierten System.

    Hmm, das ist eine Idee. Das Problem ist, dass ich derzeit nicht viel über die Funktion aussagen kann, erst recht nicht über einen möglichen Grundtyp. Und es ist auch nicht nur eine Funktion, sondern ich werde das Problem oftmals mit vielen ähnlichen Funktionen lösen müssen (darum frag ich ja auch nach einem möglichst allgemeinen Algorithmus, sonst könnte ich es auch einmal mit viel Geduld von Hand machen). Aber ich könnte mal über das Wochenende für ein paar der Funktionen mal ein paar Werte über den ganzen Bereich verteilt angucken. Vielleicht gibt es da ja schöne Muster.

    0
  • K

    Oder möchtest du mir etwas konkretes mit dem Stichwort sagen

    Nimm als Approximation die Taylorentwicklung.

    0
  • O

    mal ein anderer Ansatz, vielleicht hilft der dir:

    da du ja nicht nur einen Wert, sondern eine Kette (x_0,...,x_n) brauchst, kannst du eigentlich nach all diesen Punkten parallel suchen

    das heißt du wirfst einen Optimierungsalgorithmus auf der Funktion

    ∑i=1n(F(E(xi−1),E(xi))−C)2\sum\limits_{i=1}^n (F(E(x_{i-1}),E(x_i))-C)^2
    an. An der Stelle könntest du dann klassische multidimensionale Optimierungsverfahren wählen. Eventuell brauchst du noch einen Regularisierungsterm der dir versichert, dass die x_i monoton steigend sind oder musst das Problem passend umformulieren (zum Beispiel als x_i=x_i−1+y_i,y_i>0x\_i=x\_{i-1}+y\_i, y\_i > 0).

    Wenn du nach vielen x_i suchst, könntest du dann einfach einen effizienten Suchalgorithmus wie die CMA anwerfen. die CMA hätte den Vorteil, dass du alle Evaluationen eines Optimierungsschritts parallelisieren kannst. Andererseits brauchst du nur ca n\sqrt{n} parallele Evaluationen.

    0
  • Mod

    @otze: Das klingt schon einmal cool, daran hatte ich noch nicht gedacht. Wird eine Weile dauern bis ich mir die Stichworte angelesen habe. Aber wenn deine Beschreibung stimmt, dann geht dein Vorschlag mein Hauptproblem an, was natürlich recht günstig ist.

    0
  • C

    SeppJ schrieb:

    aber das was du wissen möchtest:

    Mich würde tatsächlich auch die Funktion interessieren. 🙂

    0
  • Mod

    ⚠ Achtung: Offtopic! ⚠

    cooky451 schrieb:

    SeppJ schrieb:

    aber das was du wissen möchtest:

    Mich würde tatsächlich auch die Funktion interessieren. 🙂

    Naja, das ist dann so etwas wie:
    E(x) = (<H(x)>, σH(x)\sigma_H(x) )
    Ich kann mit Latex leider keine spitzen Klammern schreiben. Das soll heißen: E(x) ist das Paar aus Erwartungswert von H(x) und Varianz von H(x), wobei der Erwartungswert von H(x) definiert ist als:
    <H> = −∂∂βlnZ(β)-\frac{\partial}{\partial \beta} \ln Z(\beta) bei vorgegebener Konstante β\beta (Einheit: Inverse Energie). Das Z ist wiederum
    Z(\beta) = \frac{1}{\text{Vorfaktor}} \int e^{-\beta E(\vec r\_1, \vec r\_2, \dots, \vec r\_N, V, x)}\text{d}^3 \vec r\_1 \text{d}^3 \vec r\_2 \dots \text{d}^3 \vec r\_N \text{d} V
    Das E ist jetzt ein anderes E als oben. Ich will mich aber an die gängige Konvention halten. Daher heißt das hier jetzt auch E, denn es ist eine Energie. Das E ist:
    E(\vec r\_1, \vec r\_2, \dots, \vec r\_N, V, x) = \sum\_{i,j} \text{Potential1_{i,j}}(\vec r\_i, \vec r\_j) + \sum_{i,j} \text{Potential2_{i,j}}(\vec r\_i, \vec r\_j) + \sum_{i,j,k}\text{Potential3_{i,j,k}}(\vec r\_i, \vec r\_j, \vec{i,j,k}) + E_{\text{V}}(V) + E_{\text{x}}(x, \vec r\_1, \vec r\_2, \dots, \vec{r_N})
    Das Z hat nun aber 3N+1 Freiheitsgrade, das Integral kann man niemals ausrechnen. Daher macht man einen Trick: Die meisten Zahlen in dem Exponenten sind sehr groß, tragen also kaum zum Integral bei. wenn man nun herausfinden könnte, welche E die kleinsten sind, dann kann man mit diesen das Z ziemlich gut abschätzen. Außerdem ist hier mein <H> (bis auf Vorfaktor, der dann passend gewählt ist) identisch mit dem Erwartungswert der E (das zweite E), jeweils gewichtet mit dem e−βEe^{-\beta E}. wenn ich nun also Konfigurationen der r⃗i\vec r_i und VV erzeugen kann, die mir eine Reihe von E's (das zweite E) liefern, die schon nach diesem Faktor gewichtet sind, dann kann ich einfach den Mittelwert und Varianz über diese E's (das zweite E) nehmen um mein erstes E zu erhalten. Das kann man aber recht einfach erreichen, wenn man eine gegebene Konfiguration (\vec r\_1, \dots ,\vec r\_N, V) leicht verändert, jeweils das E (das zweite) berechnet (nennen wir die Werte E1E_1 und E2E_2) und die neue Konfiguration mit einer Wahrscheinlichkeit von P=max(1,e−β(E_2−E_1))P=\text{max}(1, e^{-\beta(E\_2 -E\_1)}) akzeptiert. Wenn man so eine Kette von E_i erzeugt, dann kann man beweisen, dass deren Verteilung (bei unendlich vielen Schritten) genau den gewünschten Kriterien entspricht. Aber mit ziemlich vielen Schritten erhält man auch schon ganz gute Werte :p .

    Und das war klassische statistische Physik im kanonischen Ensemble für Einsteiger mit einem kleinen Exkurs in Monte Carlo Simulation. Vermutlich hat niemand, der es nicht schon kennt, auch nur ein Wort verstanden 🙂 . Und die, die es schon kennen, kreiden mir jetzt die ganzen Vereinfachungen und Ungenauigkeiten an, die ich gemacht habe, um einigermaßen verständlich zu bleiben.
    **
    Ich will es daher nochmal in Worte fassen:** Die Funktion E, um die es in diesem Thread geht, ist der Erwartungswert und die Schwankung einer Messung der Energie eines klassischen Systems von sehr vielen Teilchen mit gewissen Wechselwirkungen. Das kann im Allgemeinen niemand exakt berechnen, da man hundertausende Freiheitsgerade hat, aber die Computerphysik kennt sehr gute Näherungslösungen, die aber auch eine ganze Weile zum Rechnen benötigen.

    0
Anmelden zum Antworten