Trovare le radici di un gran numero di funzioni con una variabile

Question

Sto lavorando con Python/numpy/scipy per scrivere un tracciante di piccolo raggio. Le superfici sono modellate come funzioni bidimensionali che danno un'altezza sopra un piano normale. Ho ridotto il problema di trovare il punto di intersezione tra raggio e superficie per trovare la radice di una funzione con una variabile. Le funzioni sono continue e continuamente differenziabili.

C'è un modo per farlo in modo più efficiente rispetto al semplice ciclo su tutte le funzioni, utilizzando cercatori di root scipy (e magari utilizzando più processi)?

Modifica: le funzioni sono la differenza tra una funzione lineare che rappresenta il raggio e la funzione di superficie, vincolata a un piano di intersezione.

Answer 1

L'esempio seguente mostra il calcolo delle radici per 1 milione di copie della funzione x ** (a + 1) - b (tutte con differenti a e b) in parallelo utilizzando il metodo di bisezione. Prende circa ~ 12 secondi qui.

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 

L'idea di base consiste nell'eseguire semplicemente tutte le fasi usuali di un cercatore radice in parallelo su un vettore di variabili, utilizzando una funzione che può essere valutata su un vettore di variabili e vettoriale (s) equivalente di parametri che definiscono le singole funzioni dei componenti. Le condizioni sono sostituite con una combinazione di maschere e numpy.where(). Questo può continuare fino a quando tutte le radici sono state trovate alla precisione richiesta, o in alternativa finché non sono state trovate abbastanza radici che vale la pena rimuoverle dal problema e continuare con un problema più piccolo che esclude quelle radici.

Le funzioni che ho scelto di risolvere sono arbitrarie, ma aiuta se le funzioni sono ben educate; in questo caso tutte le funzioni della famiglia sono monotone e hanno esattamente una radice positiva. Inoltre, per il metodo di bisezione abbiamo bisogno di supposizioni per la variabile che dà segni diversi della funzione, e anche in questo caso è abbastanza facile venire in mente (i valori iniziali di x0 e x1).

Il codice sopra riportato utilizza forse il root finder più semplice (bisezione), ma la stessa tecnica potrebbe essere facilmente applicata a Newton-Raphson, Ridder's, ecc. I meno condizionali ci sono in un metodo di ricerca radice, il più adatto è a questa. Tuttavia, sarà necessario reimplementare qualsiasi algoritmo desiderato, non è possibile utilizzare direttamente una funzione di ricerca radice di libreria esistente.

Il frammento di codice sopra riportato è scritto con chiarezza in mente, non velocità. Evitando la ripetizione di alcuni calcoli, in particolare valutare la funzione solo una volta per iterazione invece di 3 volte, accelera questo fino a 9 secondi, come segue:

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

Per confronto, utilizzando scipy.bisect() per trovare root alla volta richiede ~ 94 secondi:

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)