Pouvez-vous paralléliser le calcul Jacobien inexact dans Julia sans tableaux spéciaux?

Question

Dans Julia, vous vouliez calculer un jacobien inexact basé sur une fonction vectorielle, f (x), qui nécessite beaucoup de calculs à évaluer. L'évaluation du jacobien est évidemment parfaitement parallèle au concept. Ma question est, cela peut-il être fait dans Julia sans avoir recours à DistributedArray, SharedArray, etc?

Par exemple, supposons que vous avez le code:

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 

Est-il possible de paralléliser ce de la même manière que vous pourriez paralléliser la somme de 200.000.000 pièce aléatoire flips, selon le manuel? Autrement dit, quelque chose d'équivalent à

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

J'ai essayé:

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 

où « testfunc.jl » est le nom du fichier dans lequel ce code, et la définition de f lui-même, se trouve . Quand j'ai essayé cela, avec f évaluant simplement à x.^2 + cos (x), j'ai pu obtenir une matrice (diagonale) correcte, mais les valeurs ne correspondaient pas à celles données par le code non parallèle (que je peut confirmer être les valeurs correctes). Une étude plus poussée suggère que le Jacobien résultant a des valeurs multipliées par 2 ou 3, en utilisant julia -p 4.

L'approche que j'ai décrite est-elle plausible (et nécessite simplement un peaufinage pour éviter la duplication des évaluations)? Sinon, y a-t-il une autre méthode par laquelle je peux évaluer le jacobien sans utiliser les types de tableaux spéciaux plus compliqués?

Il semble que l'ajout de "J = zéros (n, n)" comme première opération à l'intérieur de la boucle parallèle corrige ce problème de duplication. Est-ce que la même chose peut être faite sans avoir recours à une telle effraction brute de la matrice J?

Answer 1

Ce que je comprends du code ci-dessus est que, lorsque vous écrivez:

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 envoie une copie de J à nouveau processus évalue ensuite f(x) et somme les résultats ensemble. Je pense que la meilleure et plus efficace est d'empêcher l'envoi d'J entre les threads, et procédez comme suit:

@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 

Avec le code ci-dessus chaque thread travaille sur un nouveau J et ainsi sommation retourne la bonne réponse.