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:

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}];
ParallelDo[parallelkernel[[j]] = Eigensystem[matrix], {j, 4}];

(*Check Eigenvalues*)    
(* =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*)

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?

  • $\begingroup$ 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. $\endgroup$ – 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}} *)

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

If we set the main kernel to use only one thread, like the subkernels do, then the difference in results vanishes.

Set it using

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

My conclusion is that the difference in results is due to the different order in which arithmetic operations are carried out when doing the computation serially vs in parallel. Floating-point arithmetic is not associative, so the order of operations affects the roundoff errors.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.