Znaleźć korzenie dużej liczby funkcji z jedną zmienną

Question

Pracuję z Python/numpy/scipy, aby napisać mały znacznik ray. Powierzchnie są modelowane jako dwuwymiarowe funkcje dające wysokość powyżej normalnej płaszczyzny. Zmniejszyłem problem znalezienia punktu przecięcia promienia i powierzchni z odnalezieniem korzenia funkcji z jedną zmienną. Funkcje są ciągłe i ciągle różniczkowalne.

Czy istnieje sposób, aby zrobić to wydajniej niż po prostu zapętlić wszystkie funkcje, używając wyszukiwarki scipy (a może przy użyciu wielu procesów)?

Edytuj: Funkcje są różnicą między funkcją liniową reprezentującą promień i funkcję powierzchni, ograniczoną do płaszczyzny przecięcia.

Answer 1

Poniższy przykład pokazuje obliczanie korzeni dla 1 miliona kopii funkcji x ** (a + 1) - b (wszystkie z różnymi a i b) równolegle przy użyciu metody bisekcji. Trwa około ~ 12 sekund.

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 

Podstawowym założeniem jest to, aby po prostu działać wszystkie typowe etapy Finder głównego równolegle do wektora zmiennych za pomocą funkcji, która może być oceniona na wektorze zmiennych oraz równoważnych wektora (ów) parametrów które definiują poszczególne funkcje komponentu. Warianty zastępowane są kombinacją masek i numpy.where(). Można to kontynuować, dopóki wszystkie korzenie nie zostaną znalezione z wymaganą precyzją lub na przemian, dopóki nie zostanie znaleziona wystarczająca liczba korzeni, że warto je usunąć z problemu i kontynuować z mniejszym problemem, który wyklucza te korzenie.

Funkcje, które wybrałem do rozwiązania są arbitralne, ale pomagają, jeśli funkcje są dobrze zachowane; w tym przypadku wszystkie funkcje w rodzinie są monotoniczne i mają dokładnie jeden dodatni rdzeń. Dodatkowo, dla metody bisekcji, potrzebujemy odgadnąć dla zmiennej, która daje różne znaki funkcji, a te są również dość łatwo wymyślić tutaj również (początkowe wartości x0 i x1).

Powyższy kod wykorzystuje prawdopodobnie najprostszą wyszukiwarkę pierwiastków (bisection), ale tę samą technikę można z łatwością zastosować do Newtona-Raphsona, Riddera itp. Im mniej czynników warunkowych znajduje się w metodzie znalezienia korzenia, tym lepiej nadaje się ona do tego. Jednak będziesz musiał ponownie wdrożyć dowolny algorytm, który chcesz, nie ma możliwości bezpośredniego użycia istniejącej funkcji odnajdywania katalogu głównego biblioteki.

Powyższy fragment kodu jest napisany z myślą o jasności, a nie szybkości. Unikanie powtarzania pewnych obliczeń, w szczególności oceny funkcję tylko raz na iteracji zamiast 3 razy, to przyspiesza do 9 sekund, co następuje:

... 
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) 
... 

Dla porównania, przy użyciu scipy.bisect(), aby znaleźć jeden root w czasie ~ 94 sekund: