parallel processing - transicion - ¿Se puede paralelizar el cálculo Jacobiano inexacto en Julia sin arreglos especiales?

matriz de transicion de estados en java

Question

En Julia, quería calcular un Jacobian inexacto basado en una función vectorial, f (x), que requiere una gran cantidad de cálculos para evaluar. La evaluación del jacobiano es obviamente bastante paralela en concepto. Mi pregunta es, ¿se puede hacer esto en Julia sin recurrir a DistributedArray, SharedArray, etc.?

Por ejemplo, supongamos que tienes el código:

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

¿Es posible paralelizar esto de la misma manera que puede hacer paralelo a la suma de 200000000 saltos de moneda al azar, según el manual? Es decir, algo equivalente a

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

He intentado esto:

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

donde "testfunc.jl" es el nombre del archivo donde se encuentra este código, y la definición de f mismo. Cuando probé esto, con f simplemente evaluando x. ^ 2 + cos (x), pude obtener una matriz (diagonal) adecuada, pero los valores no coincidían con los dados por el código no paralelo (que puede confirmar que son los valores correctos). La investigación adicional sugiere que el Jacobian resultante tiene algunos valores que se multiplican por 2 o 3, cuando se usa julia -p 4.

¿El enfoque que he descrito es plausible (y simplemente requiere ajustes para evitar la duplicación de evaluaciones)? Si no, ¿hay algún otro método por el cual pueda evaluar el jacobiano sin utilizar los tipos de matriz especiales más complicados?

Parece que agregar "J = ceros (n, n)" como la primera operación dentro del bucle for paralelo corrige este problema de duplicación. ¿Se puede hacer lo mismo sin recurrir a la eliminación de la fuerza bruta de la serie J?

Answer 1

Lo que entiendo del código anterior es eso, cuando escribes:

 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 envía una copia de J al nuevo proceso y luego evalúa f(x) y suma los resultados. Creo que la forma mejor y más eficiente es evitar enviar J entre hilos, y hacer lo siguiente:

 @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 el código anterior, cada hilo funciona en una J nueva y, por lo tanto, la suma regresa la respuesta correcta.