Getting different answers when using mldivide in parfor-loop compared to in "ordinary" script

I am getting different answers for x and xp in the following code:

format long;

n = 1000;
A = rand(n);
b = ones(n,1);
x = A\b;

xp = cell(4,1);

matlabpool open local 4
parfor i=1:4,
    xp{i} = A\b;
end 

h = @(y) norm(x-y,1); 
disp(cellfun(h,xp(:)));

matlabpool close

.

An example of output at my computer is

1.0e-09 *
  0.189364648202524
  0.189364648202524
  0.189364648202524
  0.189364648202524

.

I am running Matlab 2011b with PCT version 5.2 on a 64 bit installation of Ubuntu Oneric.

The differences disappear when I start Matlab with "-singleCompThread".

Does anyone have an explanation for this behaviour?

Answer 1

[Geoff]

I had a go on my own 64-bit 2011b running on Win7 platform. Same deal. I found the following in the release notes for MatLab 4.3 (2010a) Parallel Computing Toolbox:

Other functions enhanced to support single-precision distributed arrays are chol, lu, mldivide, and eig...

Link: http://www.mathworks.com.au/help/toolbox/distcomp/rn/bsa8em4-1.html

So there is obviously some difference in how mldivide is computed between parallel and non-parallel contexts. You may be experiencing precision errors that are amplified by one of the methods. I did notice the snippet above refers to 'single-precision'. Could it be that your doubles are internally treated as singles?

MatLab says that the parallel toolbox only makes use of the outermost parfor loop. I suspect this relates to internally parallelised functions too.

With no thread pool open, I tested the following code:

while true
  x = A\b;
end

It maxed my 4 cores to 100%. So it's obviously parallel - if in doubt, compare that against something non-parallel like "x = 1+1".

Then I did this:

matlabpool open local 2
while true
  parfor i=1:2
    x=A\b;
  end
end
matlabpool close

All four cores maxed out to 50%. This suggests to me that the parallel mldivide might be disabled when invoked from a parfor. This is presumably the same effect as using the "-singleCompThread" switch.

So I would first consider that this is numeric precision error between two different algorithms, although the differences do seem to be a couple of orders of magnitude higher than I would expect.

It may be worth getting Mathworks to cross-validate the two algorithms.

-g-

[Jiro Doke]

The release note above does not refer to the backslash (mldivide) inside PARFOR. The note is about "distributed arrays" which are not used with PARFOR.

The difference is most likely due to the difference in the mutithreaded and single-threaded algorithms. Regular MATLAB uses BLAS for linear algebra operations, and it is multithreaded. But each MATLAB worker runs single-threaded. That's why you get the same results when you run MATLAB single-threaded.

[Anders Hoff]

Thanks for your answers.

This seems to make sense to me. The only thing is that I feel like the difference should have been somewhat close to eps('double').

You can confirm that implicit parallelism is off in workers with 'maxNumCompThreads'. This function is deprecated now though.

[Anders Hoff]

accepted geoff's answer, but read jiro doke's comment as well.