Lösungsalgorithmen für nichtlineare Gleichungssystem

[SimPy]

Kurzfassung:
===============
Ich versuche verschiedene Lösungsalgorithmen wie newton(), broyden1(), broyden2(), broyden3(), broyden_generalized(), anderson() und anderson2() auf ein Gleichungssystem anzuwenden, das in meiner Simulation auftritt. Mein Ziel ist es, die jeweils benötigten Dauern zur Lösung des Gleichungssystems miteinander zu vergleichen.

Allgemeine Problembeschreibung:
=========================
In Abhängigkeit von meiner Simulation ändert sich das Gleichungssystem. Das bedeutet, dass ich das Gleichungssystem, welches ich Lösen will, dynamisch (also während der Programmausführung) erzeugen muss. Ich habe ein gehöriges Problem, ein dynamisch erzeugtes Gleichungssystem mit den oben genannten Lösungsalgorithmen zu lösen.

Konkrete Problembeschreibung:
=========================
Was ich auf jeden Fall kann, um folgendes, einfaches Gleichungssystem
x_0 * cos (x_1) = 4
x_0 * x_1 − x_1 = 5
zu lösen, ist Folgendes:

Code: Alles auswählen

def func1(x):  
  out1 = [x[0]*cos(x[1]) - 4]  
  out1.append(x[1]*x[0] - x[1] - 5)  
  return out1
x01 = fsolve(func2, [1, 1], StringCommand1)

Das funktioniert auch prima. Jetzt ist es allerdings so, dass ich mehrere tausende, wenn nicht sogar hundert tausende "append" Vorgänge habe durch mein Gleichungssystem. Daher möchte ich diesen "append" Vorgang nicht in der Funktion func1(x) haben, da die ja iterativ aufgerufen wird durch fsolve() (das sieht man recht gut wenn man in func1() eine "print" Ausgabe implementiert). Dadurch würde nämlich allein durch das Zusammensetzen der Gleichungen eine größere Dauer entstehen, was den Vergleich verfälschen würde. Denn ich will ja nicht die Dauer der "append" Vorgänge messen, sondern die der Lösung im Algorithmus. Also habe ich folgende Lösung gefunden:

Code: Alles auswählen

StringCommand1 = "f = [x[0]*cos(x[1]) - 4, x[1]*x[0] - x[1] - 5]"
def func2(x,StringCommand1):
  exec StringCommand1
  return f
x01 = fsolve(func2, [1, 1], StringCommand1)

"Works like a charm", wie man so schön sagt. Denn je nach Gleichungsystem kann ich nun "StringCommand1" dynamisch zur Laufzeit einmal erzeugen und dann jedes Mal ohne weiteren Zeitverlust an die Funktion func2() übergeben! So weit, so gut. Nun habe ich mich mal an den broyden1() Algorithmus gewagt und folgendes Problem festgestellt. Ich kann die Eingabevariable "StringCommand1" nicht mehr an die Funktion func2() übergeben! Ich verdeutliche das in dem Folgenden Code:

Code: Alles auswählen

StringCommand1 = "f = [x[0]*cos(x[1]) - 4, x[1]*x[0] - x[1] - 5]"
def func2(x,StringCommand1):
  exec StringCommand1
  return f
x02 = scipy.optimize.broyden1(func2, [1,1], StringCommand1, verbose=False)

Es kommt zu folgender Fehlermeldung: "TypeError: func2() takes exactly 2 arguments (1 given)". Ich habe also folgende Zeilen ausprobiert:

Code: Alles auswählen

x02 = scipy.optimize.broyden1(func2, ([1,1], StringCommand1), verbose=False)

ebenso

Code: Alles auswählen

x02 = scipy.optimize.broyden1(func2, [[1,1], StringCommand1], verbose=False)

leider funktioniert das alles nicht. Daher meine Frage an Euch:

WIE KRIEGE ICH MIT DER FUNKTION broyden1() MEINEN STRING "StringCommand1" AN DIE FUNKTION func2() übergeben???

[BlackJack]

@SimPy: Also mal davon abgesehen, dass ich das ``exec`` gruselig finde, sollte das mit `functools.partial()` möglich sein. Damit kannst Du aus Deiner Funktion mit den zwei Argumenten ganz einfach eine erstellen, die nur ein Argument haben möchte, in dem das andere auf einen festen Wert gesetzt wird.

Wobei ich schon das erste Codebeispiel ziemlich eigenartig finde. Die Funktion kann man völlig ohne `append()` schreiben, denn die ist äquivalent zu dem hier:

Code: Alles auswählen

def func1(x):  
    return [x[0] * cos(x[1]) - 4, x[1] * x[0] - x[1] - 5]

Ich weiss nicht warum Du das so kompliziert geschrieben hast. Das sollte auch schneller sein als erst eine Zeichenkette mit Python-Code mit ``exec`` zu kompilieren und auszuführen. Was sehr wahrscheinlich *langsamer* als das mit dem `append()` sein dürfte.

[SimPy]

@ BlackJack: erst mal Danke für Deine Antwort!

Die Funktion kann man völlig ohne `append()` schreiben, denn die ist äquivalent zu dem hier

Das ist schon richtig. Aber erst während der Simulation wird das Gleichungssystem "dynamisch erstellt" und dessen Größe ist zuvor nicht bekannt. Daher kann ich nicht eine "feste" Funktion vorgeben, wie Du das vorgeschlagen hast. Stattdessen brauche ich eine Implementierung, welche die Funktion dynamisch während der Laufzeit erstellt. Deshalb hatte ich das "append" als eine mögliche Lösung angeführt. Weiterhin werden Bestandteile der Gleichungen aus (sehr großen Arrays) extrahiert: nämlich die Vorfaktoren. Bspw. steht in einer Gleichung: a*x[o] während a=Array, und die Variable i wird dann wiederum für tausende von Gleichungen durchiteriert. Wenn man das mehrere hundert bis tausende Male in mehreren hundert bis tausenden Gleichungen macht, dann wird das schon sehr langsam. Die benötigten Rechenoperation nehmen ja quadratisch zu mit n^2. Das heißt in n Zeilen meines quadratischen Gleichungssystems muss ich jeweils n Variabeln aus dem entsprechenden Array extrahieren. Nehmen wir an n=1000. Dann habe ich in z.B. 100 Iterationsschritten (wieviel es genau sind, hängt vom Lösungsalgorithmus ab) 100*1000^2 Arrayzugriffe. Das ist eine Menge. Mache ich das einmal vorher und schreibe das Gleichungssystem in einen String, den ich über "exec" ausführen kann, muss ich das nie wieder machen.

Das sollte auch schneller sein als erst eine Zeichenkette mit Python-Code mit ``exec`` zu kompilieren und auszuführen.

@ BlackJack: Da gebe ich Dir Recht! Ich habe versucht direkt eine Funktion zu übergeben anstatt des entsprechenden Strings, leider habe ich das nicht geschafft. Mir scheint, Python lässt das nicht zu. Weißt Du, wie ich eine Funktion func(), die ich während der Programmlaufzeit zuvor erstellt habe, dem Befehl "scipy.optimize.broyden1" übergeben kann? Sozusagen als "Argument" (was natürlich ein bisschen im Kreis gedacht ist, das weiß ich schon Aber dadurch würde ich mir die "append" Aufrufe innerhalb der Iterationen von "scipy.optimize.broyden1" ersparen).

sollte das mit `functools.partial()` möglich sein

@ BlackJack: Ich habe die Funktion nachgelesen, jedoch steige ich gerade nicht durch, wie ich es auf mein Problem anwenden kann?

Ich habe weiterhin folgende Möglichkeit gefunden:

Code: Alles auswählen

StringCommand2 = ["f = [x[0]*cos(x[1]) - 4, x[1]*x[0] - x[1] - 5]"]
solve = scipy.optimize.root(func2, [1,1], args=StringCommand2, method="broyden1")

Der optionale Parameter "args" erlaubt es mir einen Tupel zu übergeben. In diesem verpacke ich dann meinen String. Ich habe diesen Weg in der SciPy Referenz Version 0.11.0 nachgelesen. In der Version 0.11.0 ist dieser Befehlsaufruf so (mit der Option verschiedene Lösungsalgorithmen über den Parameter "method" anzugeben) zum ersten Mal aufgeführt. (Zitat: "Find a root of a vector function. New in version 0.11.0")
??? Weiß jemand von Euch, ob es einen Unterschied macht, ob ich den Lösungsalgorithmus "broyden1" nun direkt aufrufe in Form von "scipy.optimize.broyden1" oder über "scipy.optimize.root(method="broyden1")" ??? (In der SciPy Referenz Version 0.11.0 habe ich dazu leider keine Anmerkungen gefunden)

[BlackJack]

@SimPy: Die feste Funktion ohne `append()` die ich angegeben habe, hat *exakt den selben Effekt* wie die mit dem `append()`, die Du angegeben hast. Wenn meine Funktion keine Lösung ist, dann war es Deine auch nicht. Dann frage ich mich jetzt allerdings warum Du die überhaupt geschrieben und hier gezeigt hast!?

Natürlich kann man in Python Funktionen als Argumente übergeben. Das tust Du doch sogar schon mit `func1` und `func2`! Die `append()`-Aufrufe kannst Du Dir wie schon gezeigt ganz einfach sparen, denn das machte in Deiner Beispielfunktion überhaupt keinen Sinn.

Zu `partial()`: Ich dachte Du hast eine Funktion die zwei Argumente entgegennimmt und möchtest gerne eine die nur noch eines entgegennimmt und für das andere einen festgelegten Wert verwendet. Dafür ist `partial()` da.

Ich sehe aber weiterhin nicht wozu diese ganzen Umstände gut sein sollen. Und habe gerade wegen der ersten Funktion mit dem sinnlosen `append()` auch das Gefühl das beschriebenes Problem und gezeigter Quelltext nicht wirklich zusammen passen.

[SimPy]

@ BlackJack: OK, dann mache ich ein einfaches Beispiel: Du lässt ein Programm zwei Mal laufen. Beim ersten Mal erzeugt das Programm ein Gleichungssystem bestehend aus 5 Gleichungen. Beim zweiten Mal erzeugt das Programm ein Gleichungssystem bestehend aus 10 Gleichungen. Das Gleichungssystem soll gelöst werden. Du weißt aber nicht, ob es 5, oder 10 oder beliebig viele Gleichungen sind, denn das wird erst während der Codeausführung entschieden.

Daher musst Du die entstehende Funktion func(x) "bauen". Hat das Gleichungssystem 5 Gleichungen, dann würde man bspw. mit einer for-Schleife 5 einzelne Gleichungen zusammenstellen zu einem Gleichungssystem. Hierzu kann man bspw. das "append" benutzen. Sind es eben 10 Gleichungen, dann läuft die for-Schleife eben von 0 ... 9.

Ich kann im Quellcode nicht einfach eine einzelne, "starre" Zeile schreiben, die alle Gleichungen mit einem Komma trennt. Stattdessen muss ich meinen Quellcode so schreiben, dass er für eine BELIEBIGE Anzahl von Gleichungen funktioniert. Ich weiß nicht, wie ich das ohne eine Schleifenkonstruktion schaffen soll. Das "append" ist hierbei lediglich eine Möglichkeit, wie ich innerhalb der Schleife die einzelnen Gleichungen zu einem großen Gleichungssystem zusammensetzen kann. Es geht ebenso mit der Zuweisung der Gleichungen zu Positionen eines Arrays, der letzlich alle Gleichungen enthält. Wie auch immer man es macht, es wäre ideal das "Bauen" des Gleichungssystems nicht in jedem Iterrationsschritt von "scipy.optimize.broyden1()" zu wiederholen, sondern nur einmal außerhalb dieser Funktion zu machen. Das ist es, was ich beabsichtige. Hast Du einen Vorschlag, wie man das erreichen kann?

[Sirius3]

@SimPy: Bisher sehe ich noch nicht die Vorschrift, wie Du Dein Gleichungssystem „bauen“ willst.
In Wirklichkeit berechnest Du ja für jeden Iterationsschritt nur einen Vektor aus dem Eingangsvektor x.
Dabei macht es kaum einen Geschwindigkeitsunterschied, ob das in einer for-Schleife geschieht oder nicht.
Auf der anderen Seite macht es aber einen großen Verständnisunterschied. Programmiere am besten so, dass es am einfachsten zu verstehen ist und am einfachsten zu warten und am einfachsten zu testen.
Die meisten komplizierten Konstrukte sind nicht die Mühe wert.

PS: in Deinem jetzigen Beispiel wird in jedem Iterationsschritt der String neu geparst, was so ziemlich die langsamste Method sein dürfte, das Gleichungssystem zu berechnen.

[BlackJack]

@SimPy: Mein Vorschlag wäre es das Gleichungssystem das sich nicht verändert nur einmal zu bauen und nicht bei jedem Schritt. Lösungsansätze wäre eine Klasse, oder ein Closure, oder `functools.partial()` — was letztendlich alles mehr oder weniger äquivalente Ansätze sind. Aber wie Sirius3 schon sagte, wir wissen nicht wie Du das überhaupt machst, also kann man da wenig konkretes zu sagen. Du könntest vielleicht auch mal tatsächlich den Quelltext für eine *Funktion* erstellen und die übersetzen (`compile()`) und dann direkt übergeben.

[SimPy]

Bisher sehe ich noch nicht die Vorschrift, wie Du Dein Gleichungssystem „bauen“ willst.

@Sirius: Hier mal ein grobes Schema, das der Sache aber recht nahe kommt:

Code: Alles auswählen

# Anzahl der Gleichungen meines Gleichungssystems (beispielhaft)
n_equations = 10**20

def func1(x, ArgumentTuple):
  # func1() "baut" das Gleichungssystem variabler Größe nxn
  # Matrix_1, Matrix_2 und Matrix_3 sind bereits zuvor im Programm berechnet worden!
  function = []                         # function enthält mein Gleichungssystem
  Matrix_1 = ArgumentTuple[0]           # Extraktion der Matrix_1
  Matrix_2 = ArgumentTuple[1]           # Extraktion der Matrix_2
  Matrix_3 = ArgumentTuple[2]           # Extraktion der Matrix_3
  n_equations = ArgumentTuple[3]        # Extraktion der Anzahl der Gleichungen
  
  for i in range (n_equations):
    function.append(Matrix_3[i,i]*x[i] + Matrix_1[i,i+1]*x[i+1] + Matrix_2[i,i+2]*x[i+2] + ... + Matrix_3[i,n_equations-1]*x[n_equations-1]))
  return function
  
StartValueArray = numpy.ones(n_equations)                     # hier werden die Startwerte für den Lösungsalgorithmus geschrieben
ArgumentTuple = [Matrix_1, Matrix_2, Matrix_3, n_equations]   # Matrix_1, Matrix_2 und Matrix_3 sind bereits zuvor im Programm berechnet worden!

# Aufruf des Lösungsalgorithmus: func1(), Startwerte, benötigten Argumente und der Solver "broyden1" werden übergeben
solve = scipy.optimize.root(func1, StartValueArray, args=ArgumentTuple, method="broyden1")
print solve

Und wenn ich diesen Pseudo-Code so anschaue, weiß ich gerade selber nicht, wie ich innerhalb der for-Schleife den Teil mit "..." programmieren soll Schließlich hängt dieser Inhalt ja auch von der Variable "n_equations" ab

[pillmuncher]

@SimPy: Ich verstehe nicht recht, was du sagen willst, aber ich vermute, du willst eine List Comprehension:

Code: Alles auswählen

>>> [x ** 2 for x in range(10)]
[0, 1, 4, 9, 16, 25, 36, 49, 64, 81]
>>> [x ** 2 + y for x in range(3) for y in range(3)]
[0, 1, 2, 1, 2, 3, 4, 5, 6]
>>> [[x ** 2 + y for x in range(3)] for y in range(3)]
[[0, 1, 4], [1, 2, 5], [2, 3, 6]]

Und da du anscheinend Python 2.x verwendest, willst du sicherlich xrange() statt range() verwenden.

[BlackJack]

@SimPy: Zusätzlich zu den Anmerkungen von pillmuncher: Da es hier um `numpy` geht, sollte man als erstes überlegen ob man diese Schleife nicht in die `numpy`-Datentypen verlegen kann. Bevor man anfängt ”exotische” Lösungen zu suchen.

Und die Beschreibung war auch die ganze Zeit falsch. Das was Du da an `function` bindest und zurück gibst, ist ja gar keine Funktion, sondern das *Ergebnis* der Berechnung.

Das `ArgumentTuple` kann man kompakter auf Namen verteilen (mit Namen die sich mehr an die Konventionen halten):

Code: Alles auswählen

    matrix_1, matrix_2, matrix_3, equation_count = arguments

Wobei ich mich gerade frage ob die Anzahl der Gleichungen überhaupt übergeben werden muss, denn ich nehme mal an, dass die sich aus der Grösse der Matrizen ergibt!?

Konkrete Datentypen haben in Namen nichts zu suchen. Wenn man den Datentyp dann mal ändert, und das kommt durchaus vor, muss man entweder den Namen überall ändern, oder man hat falsche, irreführende Namen im Programm. Schönes Beispiel ist `ArgumentTuple` das bei Dir an eine *Liste* gebunden ist.

[SimPy]

@pillmuncher wegen "xrange": VIELEN DANK! Das habe ich einfach nicht gekannt. xrange geht tatsächlich schneller.
@pillmuncher wegen ListComprehension: Ich verstehe leider nicht, wie ich die List Comprehension auf mein Problem anwenden kann. Was ich brauche ist eine Implementierung, die eine Funktion f(x) ausgibt. Denn "scipy.optimize.root()" verlangt als erstes Argument eine Funktion f(x) mit einer Variablen x. In Deinem Beispiel jedoch erzeugst Du, wenn ich das richtig verstehe, nur numerische Werte.
@BlackJack: Deine Anmerkungen zum Stil sind richtig. Wie gesagt, es handelt sich um PseudoCode, der das Prinzip veranschaulichen soll.

Vielleicht stehe ich auf dem Schlauch!? Dann zeigt mir einfach mal, wie Ihr Folgendes Implementiert: Nehmen wir an, n_equations=3. Dann brauche ich Folgenden Quellcode:

Code: Alles auswählen

def f(x)
  func = x[0] + x[1] + x[2] + exp(x[1] - x[2])
  return func

Nehmen wir an, n_equations=4. Dann brauche ich Folgenden Quellcode:

Code: Alles auswählen

def f(x)
  func = x[0] + x[1] + x[2] + x[3] + exp(x[2] - x[3])
  return func

Nehmen wir an, n_equations=5. Dann brauche ich Folgenden Quellcode:

Code: Alles auswählen

def f(x)
  func = x[0] + x[1] + x[2] + x[3] + x[4] + exp(x[3] - x[4])
  return func

Der Quellcode muss variabel sein, so dass n_equations einen beliebigen Wert annehmen kann. Wie implementiert Ihr so etwas?

[BlackJack]

@SimPy: Wie schon mal gesagt der Name `func` ist falsch und irreführend, denn das ist keine Funktion sondern eine Zahl die Du an diesen Namen bindest. Zumal der Name hier unnnötig ist, denn man könnte den Ausdruck auch gleich hinter das ``return`` schreiben.

Wenn ich das Muster hier richtig erkannt habe, dann ist das ganz einfach das hier:

Code: Alles auswählen

def f(x)
    return sum(x[:-2]) + exp(x[-2] - x[-1])

Edit: Falls sicher ist, das `x` ein `numpy`-Array ist, würde ich die Berechnung natürlich so schreiben:

Code: Alles auswählen

    return x[:-2].sum() + exp(x[-2] - x[-1])

[SimPy]

@BlackJack: ICH VERNEIGE MICH VOR DIR! Das war ein riesen Schritt nach vorne! Ich habe nun folgenden Pseudo-Code:

Code: Alles auswählen

def fun(x,M):                                   
  n = len(M)            # Meine Matrix ist immer quadratisch --> Anz. der Gleichungen
  Result = []           # Liste, die am Ende ausgegeben wird
  Result = [sum( M[index,:n]*x[:n] + exp(x[-2] - x[-1]) ) for index in xrange (n) ]   # List Comprehension
  return Result

Ich habe die List comprehension verwendet (wie von pillmuncher vorgeschlagen), damit ich ich keine for-Schleife benutzen muss. Oder geht es doch schneller, wenn ich eine for-Schleife benutze und dann die Werte in einen Array "Result" schreibe?

[derdon]

In einer for-Schleife kann vieles passieren; innerhalb einer List-Comprehension ist die erlaubte Syntax schon eingeschränkter. Allein deswegen kann man sich schon vorstellen, dass bei einer LC mehr optimiert werden kann, weil eben mehr Annahmen über den Code gemacht werden können.

[BlackJack]

@SimPy: Die Zuweisung der leeren Liste an `Result` kann weg. Du bindest den Namen doch gleich in der nächsten Zeile an eine andere Liste ohne die leere Liste jemals benutzt zu haben.

Ob das Argument von `sum()` so funktioniert, hängt sehr stark von den Werten und deren Typen ab. Ich vermute mal das wird so nicht funktionieren beziehungsweise nicht das gewünschte Ergebnis liefern. Das „slicen” mit ``:n`` macht keinen Sinn, denn das ist per Definition immer die gesamte Sequenz, d.h. ``M[index, :len(M)]`` ist genau das gleiche wie ``M[index]`` — nur unnötig kompliziert ausgedrückt. Bei `x` ist die Frage welche Form/Länge das hat.

Wenn Du Geschwindigkeit willst, dann lass Python-Schleifen möglichst ganz weg, denn letztendlich ist auch eine „list comprehension” eine Schleife die in Python ausgeführt wird, und verwende `numpy`-Arrays. Also nicht nur als Datencontainer sondern tatsächlich auch als Datentypen mit den Operationen die sie zur Verfügung stellen. Rechne mit ganzen Arrays statt mit Einzelwerten in Schleifen.

[SimPy]

Wenn Du Geschwindigkeit willst, dann lass Python-Schleifen möglichst ganz weg ... und verwende `numpy`-Arrays. Also nicht nur als Datencontainer sondern tatsächlich auch als Datentypen mit den Operationen die sie zur Verfügung stellen.

@BlackJack: Genau das hat alles deutlich erleichtert und den Quellcode wesentlich verkürzt. Ich habe das nun so gemacht, wie Du gesagt hast. Es gibt jetzt keine Listen mehr. Ich habe stattdessen eine Vektormultiplikation eingeführt:

Code: Alles auswählen

dot(M,x)     # ganz schön kurz auf einmal, nicht wahr? :o)

Meine Exponentialfunktionen füge ich erst nach der Vektormultiplikation ein. Mache ich alles nur mit den Array-Operationen, die ich im SciPy-Tutorial gefunden habe. Ich werde erst in den kommenden Tagen sehen, ob ich alles korrekt implementiert habe. Für's erste danke ich Euch aber ganz arg für Eure SEHR HILFREICHEN Kommentare!

Euer SimPy