Suche nach den Wurzeln einer großen Anzahl von Funktionen mit einer Variablen

Question

Ich arbeite mit Python/numpy/scipy, um einen kleinen Ray Tracer zu schreiben. Oberflächen werden als zweidimensionale Funktionen modelliert, die eine Höhe über einer normalen Ebene ergeben. Ich reduzierte das Problem, den Schnittpunkt zwischen Strahl und Oberfläche zu finden, um die Wurzel einer Funktion mit einer Variablen zu finden. Die Funktionen sind kontinuierlich und kontinuierlich differenzierbar.

Gibt es eine Möglichkeit, dies effizienter zu tun, als einfach alle Funktionen mit scipy root finders zu durchlaufen (und vielleicht mehrere Prozesse zu verwenden)?

Bearbeiten: Die Funktionen sind der Unterschied zwischen einer linearen Funktion, die den Strahl darstellt, und der Oberflächenfunktion, die auf eine Schnittebene beschränkt ist.

Answer 1

Das folgende Beispiel zeigt die Berechnung der Wurzeln für 1 Million Kopien der Funktion x ** (a + 1) - b (alle mit unterschiedlichem a und b) parallel unter Verwendung der Bisektionsmethode. Dauert ungefähr ~ 12 Sekunden hier.

import numpy 

def F(x, a, b): 
    return numpy.power(x, a+1.0) - b 

N = 1000000 

a = numpy.random.rand(N) 
b = numpy.random.rand(N) 

x0 = numpy.zeros(N) 
x1 = numpy.ones(N) * 1000.0 

max_step = 100 
for step in range(max_step): 
    x_mid = (x0 + x1)/2.0 
    F0 = F(x0, a, b) 
    F1 = F(x1, a, b) 
    F_mid = F(x_mid, a, b) 
    x0 = numpy.where(numpy.sign(F_mid) == numpy.sign(F0), x_mid, x0) 
    x1 = numpy.where(numpy.sign(F_mid) == numpy.sign(F1), x_mid, x1) 
    error_max = numpy.amax(numpy.abs(x1 - x0)) 
    print "step=%d error max=%f" % (step, error_max) 
    if error_max < 1e-6: break 

Der Grundgedanke ist einfach auf einem Vektor von Variablen alle üblichen Schritte eines Wurzel finder parallel auszuführen, unter Verwendung einer Funktion, die auf einen Vektor von Variablen und äquivalenten Vektors (n) der Parameter ausgewertet werden können Diese definieren die einzelnen Komponentenfunktionen. Conditionals werden durch eine Kombination von Masken und numpy.where() ersetzt. Dies kann fortgesetzt werden, bis alle Wurzeln mit der erforderlichen Genauigkeit gefunden wurden, oder alternativ, bis genügend Wurzeln gefunden wurden, dass es sich lohnt, sie aus dem Problem zu entfernen und mit einem kleineren Problem fortzufahren, das diese Wurzeln ausschließt.

Die Funktionen, die ich lösen wollte, sind willkürlich, aber es hilft, wenn die Funktionen gut sind; In diesem Fall sind alle Funktionen in der Familie monoton und haben genau eine positive Wurzel. Zusätzlich benötigen wir für die Bisektionsmethode Annahmen für die Variable, die unterschiedliche Vorzeichen der Funktion ergeben, und diese sind auch hier leicht zu finden (die Anfangswerte von x0 und x1).

Der obige Code verwendet vielleicht den einfachsten Wurzelfinder (Bisektion), aber die gleiche Technik könnte leicht auf Newton-Raphson, Ridder usw. angewendet werden. Je weniger Bedingungen in einer Wurzelfindungsmethode vorhanden sind, desto besser ist sie geeignet dazu. Sie müssen jedoch einen beliebigen Algorithmus neu implementieren, da es keine Möglichkeit gibt, eine vorhandene Bibliothekswurzelfinderfunktion direkt zu verwenden.

Das obige Code-Snippet wurde mit Klarheit geschrieben, nicht mit Geschwindigkeit. Die Vermeidung der Wiederholung von einigen Berechnungen, insbesondere die Funktion Auswertung nur einmal pro Iteration anstelle von 3-mal, beschleunigt dies auf 9 Sekunden nach oben, wie folgt:

... 
F0 = F(x0, a, b) 
F1 = F(x1, a, b) 

max_step = 100 
for step in range(max_step): 
    x_mid = (x0 + x1)/2.0 
    F_mid = F(x_mid, a, b) 
    mask0 = numpy.sign(F_mid) == numpy.sign(F0) 
    mask1 = numpy.sign(F_mid) == numpy.sign(F1) 
    x0 = numpy.where(mask0, x_mid, x0) 
    x1 = numpy.where(mask1, x_mid, x1) 
    F0 = numpy.where(mask0, F_mid, F0) 
    F1 = numpy.where(mask1, F_mid, F1) 
... 

Zum Vergleich wurde unter Verwendung von scipy.bisect() zu finden, Wurzel zu einer Zeit dauert ~ 94 Sekunden:

for i in range(N): 
    x_root = scipy.optimize.bisect(lambda x: F(x, a[i], b[i]), x0[i], x1[i], xtol=1e-6)