Finding eigenvalues in Mathematica: why so slow?

I am trying to find the eigensystem of a large sparse real symmetric matrix, and I only need the lowest 40 or so eigenstates. The relevant code is as follows:

TheModel[V0_, V1_, V2_, index_, period_] := (
Fn = Fibonacci[index]; Fn1 = Fibonacci[index - 1];
size = 2.0*period*Fn + 1;
RV0 = Table[(V0/4), size - Fn];
RV1 = Table[(V1/4), size - Fn1];
central =
Table[((2 \[Pi]/Fn)*(ind - Fn*period - 1)/\[Pi])^2 + (V0 + V1)/
2, {ind, 1, size}];
t1 = AbsoluteTime[];
HH = SparseArray[DiagonalMatrix[RV0, Fn] + DiagonalMatrix[RV1, Fn1] +
DiagonalMatrix[central, 0] + DiagonalMatrix[RV0, -Fn] +
DiagonalMatrix[RV1, -Fn1]];
t2 = AbsoluteTime[];
Print["Time is ", t2 - t1];
t1 = AbsoluteTime[];
{\[Epsilon], \[Phi]} = Eigensystem[HH, -Fn];
t2 = AbsoluteTime[];
Print["Time is ", t2 - t1];
Return[{\[Epsilon], \[Phi]}];
);
{\[Epsilon], \[Phi]} = TheModel[30.0, 2, 0, 9, 20]; // Timing


The matrix has a dimension of about 1200. I am surprised that without defining HH as a sparse array, the diagonalization step alone takes 1.33 seconds; in contrast, if I keep the sparse array definition, the diagonalization only takes 0.054 seconds. In comparison, it takes Matlab 0.6 seconds to diagonalize the same matrix (no sparse array definition and returning all eigenstates).

I am very curious about why sparse array definition speeds up the diagonalization process so significantly. Can sometime give me some insights?

Also, in general why is it so slow in Mathematica to generate all eigenstates? I thought both Matlab and Mathematica use LAPACK to do the diagonalization. Why is there such a big difference?

• Could you try AbsoluteTiming instead of Timing? Please check this question first. – Pavlo Fesenko Jan 30 '16 at 14:25
• The sparse case is fast because it can use Krylov methods to find the smallest several eigenvalues and vectors. – Daniel Lichtblau Jan 30 '16 at 22:26
• @PavloFesenko Thanks for the reference! I have read that post, and found it very helpful. However, I am wondering if AbsoluteTiming is the complete answer, because I can still see that calculations in Matlab completes appreciably faster than MMA, especially when I ask for the full eigenstates from both. Is it possible that the matrix diagonalization alone takes roughly the same time, but the overhead in MMA takes much longer? – Xiao Jan 30 '16 at 22:51
• You can get a bit of a speedup if you use Method->"Arnoldi" for both the sparse and non sparse case. – user21 Feb 3 '16 at 19:19
• You can also speed up by using real numbers rather than rational. Change the def of RV0 to RV0 = Table[(V0/4.), size - Fn]; and the same for RV1. – bill s Feb 23 '16 at 3:38