Can you parallelise inexact Jacobian calculation in Julia without special arrays?

Question

In Julia, you wanted to calculate an inexact Jacobian based on a vector function, f(x), that requires a lot of calculation to evaluate. Evaluation of the Jacobian is obviously quite neatly parallelisable in concept. My question is, can this be done in Julia without resorting to DistributedArray, SharedArray, etc?

For example, suppose you have the 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

Is it possible to parallelise this in the same way that you might parallelise the summation of 200000000 random coin flips, as per the manual? That is, something equivalent to

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

I have tried this:

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

where "testfunc.jl" is the name of the file where this code, and the definition of f itself, is found. When I tried this, with f simply evaluating to x.^2+cos(x), I was able to get a proper (diagonal) matrix out, but the values did not match those given by the non-parallel code (which I can confirm to be the correct values). Further investigation suggests that the resulting Jacobian has some values being multiplied by 2 or 3, when using julia -p 4.

Is the approach I have described a plausible one (and simply requires tweaking to prevent duplication of evaluations)? If not, is there another method by which I can evaluate the Jacobian without using the more complicated special Array types?

It appears that adding "J=zeros(n,n)" as the first operation inside the parallel for loop corrects this duplication issue. Can the same thing be done without resorting to such brute force clearing of the J array?

Answer 1

What I understand from above code is that, when you write:

  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 sends a copy of J to new process then evaluates f(x) and sum results together. I think the better and more efficient way is to prevent sending J between threads, and do the following:

  @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

With the above code each thread works on a fresh J and so summation returns the right answer.