Czy możesz paralelizować niedokładne obliczenia jakobińskie w Julii bez specjalnych tablic?

Question

W Julii, chciałeś obliczyć niedokładny Jacobian w oparciu o funkcję wektorową f (x), która wymaga wiele obliczeń do oceny. Ocena jakobianina jest oczywiście dość starannie równoległa w koncepcji. Moje pytanie brzmi: czy można to zrobić w Julii bez uciekania się do DistributedArray, SharedArray, itp.?

Na przykład załóżmy, że masz kod:

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 

Czy możliwe jest parallelise to w taki sam sposób, że może parallelise sumowanie 200000000 losowej wygina monety, według instrukcji? Oznacza to, że coś równoważne

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

Próbowałem to:

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 

gdzie „testfunc.jl” to nazwa pliku, w którym ten kod, a definicja f siebie, znajduje się . Kiedy próbowałem tego, po prostu f oceniając x.^2 + cos (x), byłem w stanie uzyskać właściwą (ukośną) matrycę, ale wartości nie zgadzały się z wartościami podanymi przez nierównoległy kod (który może potwierdzić, że są prawidłowe wartości). Dalsze badania sugerują, że wynikowe wartości jakobianowe mają pomnożone wartości o 2 lub 3, używając julia -p 4.

Czy podejście, które opisałem, jest wiarygodne (i po prostu wymaga poprawienia, aby zapobiec powielaniu ocen)? Jeśli nie, czy istnieje inna metoda, dzięki której mogę ocenić Jacobian bez korzystania z bardziej skomplikowanych typów specjalnych Array?

Wygląda na to, że dodanie "J = zer (n, n)" jako pierwszej operacji w pętli równoległej poprawia ten problem z duplikowaniem. Czy można zrobić to samo, bez uciekania się do tak brutalnej wymazywania tablicy J?

Answer 1

Co rozumiem z powyższego kodu jest to, że kiedy piszesz:

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 wysyła kopię J do nowego procesu następnie ocenia f(x) i suma wyników razem. myślę, że lepszy i bardziej efektywny sposób to zapobiec wysyłaniu J pomiędzy wątkami, i wykonaj następujące czynności:

@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 

z powyższym kodem każdy wątek działa na świeżym J a więc podsumowanie zwraca prawidłową odpowiedź.