Performance and accuracy issue with eig

[jebej]

I am seeing a strange divergence between MATLAB and Julia 0.4.2 whereas the eigenvectors returned by both programs are not the same, and Julia is substantially slower. I have attached both files in a Gist here.

I am generating a bunch of transverse field Ising model hamiltonians and varying the transverse field strength h. I then look at the first component of the first eigenvector for each h value. This component represents the amplitude of the wavefuntion in the all spin down configuration.

You would expect, for h=0, that the two lowest energy eigenvectors have all the spins down and up. This is what I observe. As we increase h, the amplitude of the wavefunction in the all spin down component should decrease.

If we plot this first component of the first eigenstate as a function of h, MATLAB and Julia diverge! This is strange as I assume that they should both be using LAPACK...

Also, and this might be related, Julia is much slower than MATLAB (4sec vs. 0.8sec), with most of the time spent in the call to LAPACK in syevr!.

[KristofferC]

Matlab uses MKL so it would be interesting to build Julia with MKL and do the same experiment.

[andreasnoack]

The eigenvalues of these matrices have algebraic multiplicity larger than one so the eigenvectors are not uniquely defined. Only the subspaces for each value are uniquely defined. Numerically, they might be slightly different and then the vector are just ill-conditioned.

LAPACK returns one of many possible bases and small perturbations could change the result. A different BLAS would, in general, give a different numerical result even though the subspaces spanned are the same. Also, changing the number of threads could give different results. E.g.

julia> A = TFIM_1D_PBC(N,J,hrange[15]);

julia> blas_set_num_threads(4);

julia> _, vecs1 = eig(A);

julia> blas_set_num_threads(1);

julia> _, vecs2 = eig(A);

julia> norm(vecs1 - vecs2, 1)
21.464937248160055

With respect to the speed then I think that MKL might be using threading at the LAPACK level whereas with OpenBLAS you'll only get threading on the BLAS level.

[andreasnoack]

Just did some tests and it appears that MKL is just faster. Hard to tell what they are doing.

[jebej]

I am still not sure I understand the accuracy issue though. I am looking at a specific eigenvector, and the multiplicity of its eigenvalue is 1 as soon as the transverse field is turned on.

Note that there are errors for which the two lowest eigenvectors get interchanged since their eigenvalue is numerically very similar for small h, but I correct those variations by looking at the vector components (that's the purpose of the if block in the loop).

[andreasnoack]

Numerically, the multiplicity is effectively higher than one or the conditioning is just really bad. I get

julia> blas_set_num_threads(1);

julia> E1, V1 = diagonalize();

julia> (E1[18,2] - E1[18,1])/eps(E1[18,1])
2.0

so the two eigenvalues are almost the same and, therefore, any rotation (not just a permutation) could be applied to the first two vectors without significantly distorting the factorization. Maybe you could try a more fancy check that applies a rotation to the two first vectors when the values are very close.