ParallelDo gives different solution to Eigensystem

I am trying to calculate the eigensystem of a large matrix (e.g. 256x256). I have found that when I do this within a ParallelDo (because I am actually calculating many of these eigensystems), the result is a different one than when I am calculating it on the main Kernel. Example:

SeedRandom[1234]
matrix = RandomReal[NormalDistribution[0, 1], {256, 256}];
(* Taking a Hermitian matrix for convenience*)
matrix = matrix + ConjugateTranspose@matrix;
mainkernel = Eigensystem[matrix];
parallelkernel = Table[0, {j, 4}];
SetSharedVariable[parallelkernel];
ParallelDo[parallelkernel[[j]] = Eigensystem[matrix], {j, 4}];

(*Check Eigenvalues*)
Max@Abs@(mainkernel[[1]]-parallelkernel[[1,1]])
(*!=0*)
Max@Abs@(parallelkernel[[1,1]]-parallelkernel[[2,1]])
(* =0*)

(* Make sure the largest entry in the eigenvector is positive to get the same direction *)
signofmain = Map[Sign@Total@MinMax[#] &, mainkernel[[2]]];
signofparallel = Map[Sign@Total@MinMax[#] &, parallelkernel[[1, 2]]];

(*Check eigenvectors*)
Max@Abs@(signofmain*mainkernel[[2]]-signofparallel*parallelkernel[[1,2]])
(*!=0*)
Max@Abs@(parallelkernel[[1,2]]-parallelkernel[[2,2]])
(*=0*)


Obviously, all the calculations done within ParallelDo show the exactly same result, but a different one than on the main kernel.

I am aware of the fact that the differences are extremely small here. However, a subsequent division by the difference of two eigenvalues (in my case eigenenergies) can lead to an error amplification up to 10^-2 in some cases, which is of course non-negligible.

Where does this difference come from and how can I avoid it?

A simple example ...

... to show that eigenvalue computations give slightly different results on the main kernel and parallel subkernels.

Create a symmetric matrix:

matrix = RandomVariate[GaussianOrthogonalMatrixDistribution[256]];


Compute eigenvalues on main kernel:

res = Eigenvalues[matrix];


Compute eigenvalues on subkernels:

parRes = ParallelEvaluate[Eigenvalues[matrix]];


The eigenvalues computed on the subkernels are not exactly the same as those computed on the main kernel:

Max@Abs[res - First[parRes]]
(* 7.4607*10^-14 *)


However, all subkernels give exactly the same result:

Outer[Max@*Abs@*Subtract, parRes, parRes, 1]
(* {{0., 0., 0., 0.},
{0., 0., 0., 0.},
{0., 0., 0., 0.},
{0., 0., 0., 0.}} *)


Why is there a (small) difference between the eigenvalues computed on the main and subkernels?

• I find your example a bit more complicated than necessary, and thus requiring more effort to read through than a better-phrased question. No doubt this is why the question received no attention in 16 hours, not even an upvote, despite the strangeness of the phenomenon you observed. I added a simple and commented example. – Szabolcs Nov 20 '18 at 13:24

I believe this is because eigenvalue computation is internally parallelized (in a way that has nothing to do with the Parallel Tools), thus the precise order of arithmetic operations depends on the number of parallel threads.

The main kernel will use 4 threads by default on my 4-core machine. The subkernels only use 1 by default.

You can check this as follows:

SystemOptions["ParallelOptions" -> "MKLThreadNumber"]
(* {"ParallelOptions" -> {"MKLThreadNumber" -> 4}} *)

ParallelEvaluate@

SetSystemOptions["ParallelOptions" -> "MKLThreadNumber" -> 1]