你可以在没有特殊数组的情况下在Julia中并行不精确的雅可比计算吗？

Question

在Julia中,你想要根据向量函数f(x)来计算一个不精确的雅可比行列式,这需要大量的计算来评估.雅各比派的评价显然在概念上非常平行.我的问题是,这可以在Julia中完成而不需要使用DistributedArray,SharedArray等吗？

例如,假设您有代码:

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

是否有可能按照与手册中的200000000随机硬币翻转总和相同的方式对此进行并行化？也就是说,等同于

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

我试过这个:

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

其中"testfunc.jl"是找到此代码的文件的名称,以及f本身的定义.当我尝试这个时,用f简单地评估x.^ 2 + cos(x),我能够得到一个正确的(对角线)矩阵,但是这些值与非并行代码给出的值不相符(我可以确认是正确的值).进一步调查表明,当使用julia -p 4时,得到的雅可比行列式的某些值乘以2或3.

我所描述的方法是否合理(并且只需要调整以防止重复评估)？如果没有,是否有另一种方法可以在不使用更复杂的特殊数组类型的情况下评估雅可比行列式？

似乎在并行for循环内添加"J = zeros(n,n)"作为第一个操作可以纠正这个重复问题.如果不诉诸J阵列的强力清除,可以做同样的事情吗？

Answer 1

我从上面的代码中了解到的是，当你写：

  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 将 的副本发送J到新进程，然后评估f(x)结果并将结果相加。我认为更好、更有效的方法是防止J线程之间发送，并执行以下操作：

  @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

通过上面的代码，每个线程都在新的上工作J，因此求和返回正确的答案。