È possibile parallelizzare il calcolo Jacobiano inesatto in Julia senza matrici speciali?

Question

In Julia, si voleva calcolare un Jacobiano inesatto basato su una funzione vettoriale, f (x), che richiede molti calcoli da valutare. La valutazione dello Jacobiano è ovviamente abbastanza ordinatamente parallelizzabile nel concetto. La mia domanda è, può essere fatto in Julia senza ricorrere a DistributedArray, SharedArray, ecc?

Ad esempio, si supponga di avere il codice:

function Jacob(f::Function,x) 
    eps=1e-7 
    delta=eps*eye(length(x)) 
    J=zeros(length(x),length(x)) 
    for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    end 
    J 
end 

E 'possibile parallelizzare questo nello stesso modo in cui si potrebbe parallelizzare la somma di 200 milioni di monete a caso lancia, come da manuale? Cioè, qualcosa di equivalente a

nheads = @parallel (+) for i=1:200000000 
    int(randbool()) 
end 

Ho provato questo:

function Jacob(f::Function,x) 
    require("testfunc.jl"); 
    eps=1e-7 
    delta=eps*eye(length(x)) 
    J=zeros(length(x),length(x)) 
    J=@parallel (+) for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 
    J 
end 

dove "testfunc.jl" è il nome del file in cui questo codice, e la definizione di f per sé, si trova . Quando ho provato questo, con f semplicemente valutando x.^2 + cos (x), ero in grado di ottenere una matrice (diagonale) corretta, ma i valori non corrispondevano a quelli dati dal codice non parallelo (che io può confermare di essere i valori corretti). Ulteriori indagini suggeriscono che il risultante Jacobiano ha alcuni valori moltiplicati per 2 o 3, quando si usa julia -p 4.

L'approccio che ho descritto è plausibile (e richiede semplicemente un ritocco per evitare la duplicazione delle valutazioni)? In caso contrario, esiste un altro metodo con cui posso valutare Jacobian senza utilizzare i tipi di Array speciali più complicati?

Sembra che aggiungere "J = zeri (n, n)" come prima operazione all'interno del ciclo parallelo per correggere questo problema di duplicazione. Si può fare la stessa cosa senza ricorrere a una tale eliminazione bruta della forza della matrice J?

Answer 1

Quello che ho capito da sopra il codice è che, quando si scrive:

J=zeros(length(x),length(x)) 
    J=@parallel (+) for i=1:length(x) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 

Julia invia una copia del J al nuovo processo quindi valuta f(x) e somma i risultati insieme. Penso che il modo migliore e più efficace è quello di evitare l'invio di J tra i thread, e procedere come segue:

@parallel (+) for i=1:length(x) 
    J=zeros(length(x),length(x)) 
    J[:,i]=(f(x+delta[:,i])-f(x-delta[:,i]))/2/eps 
    J 
    end 

Con il codice sopra ogni thread lavora su una nuova J e così sommatoria restituisce la risposta giusta.