# Why am I losing precision in computation of eigenvectors of sparse matrices?

I compute the eigenvectors and eigenvalues of the same matrix in its dense and sparse version. The eigenvalues are the same to great precision, but the eigenvectors are not: when I compute the eigenvalue of the solution eigenvector, I am losing precision around the 10th decimal place. This has serious consequences for the next computation with that eigenvector (computation of entropy for those who know about quantum information physics), it gives a totally different result when computing with dense or sparse matrices.

dim = 10;

Do[
If[n == 1,
ox = {{0., 1}, {1, 0}};
oz = {{1, 0.}, {0, -1}};
id = IdentityMatrix[2],

ox = SparseArray[ox];
oz = SparseArray[oz];
id = SparseArray[id];];

\[Lambda] = 0.1;
hx = 0.; hz = 0.; h = 0.;
Do[
If[i == 1, a = ox; c = oz, a = id; c = id];
hxpart = a; hzpart = c;

Do[
If[j == i || j == i + 1, b = ox, b = id;];
hxpart = KroneckerProduct[hxpart, b];
If[j == i, d = oz, d = id;];
hzpart = KroneckerProduct[hzpart, d];
, {j, 2, dim}];
If[i == dim, hxpart = 0];
hx = hx + hxpart; hz = hz + hzpart;
, {i, 1, dim}];

h = hx + \[Lambda] hz;
nordre = Ordering[Eigensystem[h, 2][[1]], 1][[1]];
gsreala = Chop[Eigensystem[h, 2][[2, nordre]]];
Print[NumberForm[Eigenvalues[h, 4], 15]];
Print[NumberForm[Conjugate[gsreala].h.gsreala, 15]];

(*Entropy computation*)

gsm = Table[gsreala[[2^(dim/2)*i - (2^(dim/2) - 1) + j - 1]], {i, 1,
2^(dim/2)}, {j, 1, 2^(dim/2)}];

singvalues = SingularValueList[gsm];
dimsing = Length[singvalues];
eigenvalues = Table[singvalues[[i]]^2, {i, 1, dimsing}];
epsilon = 10^-10;
H[x_] := If[Norm[x] > epsilon, -x*Log[2, x], 0.];
entreal = N[Sum[H[eigenvalues[[i]]], {i, 1, dimsing}]];
Print[entreal];

, {n, 1, 2}];


These are the results I obtain:

Dense matrix:

Four smallest eigenvalues: {9.03002193787516,-9.03002193787516,-9.03002193767716,9.03002193767714}

Eigenvalue of the smallest eigenvector: -9.03002193787517

Entropy: 1.00001

Sparse matrix:

Four smallest eigenvalues:{-9.03002193787516,9.03002193787516,9.03002193767717,-9.03002193767716}

Eigenvalue of the smallest eigenvector: -9.03002193781702

Entropy: 0.26277

So I am having a loss of precision on the computation of the eigenvector, that ultimately produces a totally different result in the entropy.

• If you're going to be totally paranoid about this, I would suggest attempting to evaluate to arbitrary precision; that means avoiding machine precision expressions like 0. in your code, and using the second argument of N[], judiciously, with perhaps slightly more precision than what you need. – J. M. is away Jul 8 '16 at 14:04
• Thanks, this solved the problem! Putting a precision inside N[ ] of 3 or higher gives me the correct result. – user3154647 Jul 8 '16 at 15:03

Eigensystem[h, 2, Method -> "Direct"]
• If memory serves, Arnoldi is indeed performed by default when the matrix is sparse and only the first few or last few eigenpairs are asked for. I would not recommend the use of "Direct" in general. – J. M. is away Jul 8 '16 at 13:59
• @user, that's why I said it's not recommended in general; the dense case uses QR by default, and here the setting "Direct" forces the use of QR on your sparse matrix. – J. M. is away Jul 8 '16 at 14:10
• @user3154647, there is no general (iterative) solution for this. Try the methods in the documentation, FEAST in a band (if your matrix is Hermitian), set "Arnoldi" with a tolerance. Did you have a look at the documentation? – user21 Jul 8 '16 at 14:22