5
$\begingroup$

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?

Improve this question
$\endgroup$
6
  • $\begingroup$ Could you try AbsoluteTiming instead of Timing? Please check this question first. $\endgroup$
    – Pavlo Fesenko
    Jan 30, 2016 at 14:25
  • 4
    $\begingroup$ The sparse case is fast because it can use Krylov methods to find the smallest several eigenvalues and vectors. $\endgroup$
    – Daniel Lichtblau
    Jan 30, 2016 at 22:26
  • $\begingroup$ @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? $\endgroup$
    – Xiao
    Jan 30, 2016 at 22:51
  • 2
    $\begingroup$ You can get a bit of a speedup if you use Method->"Arnoldi" for both the sparse and non sparse case. $\endgroup$
    – user21
    Feb 3, 2016 at 19:19
  • 1
    $\begingroup$ 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. $\endgroup$
    – bill s
    Feb 23, 2016 at 3:38

0

Reset to default

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Browse other questions tagged or ask your own question.