3

I want to compute the eigenvalues (and later the corresponding eigenvectors) of an $n\times n$ Hermitian matrix. For this I use {evs, vecs} = Eigensystem[matrix] or the command Eigenvalues to evaluate only the eigenvalues, but this takes a long time when I already do it for a $40\times 40$ matrix, and that 30 times (I am computing energy bands and thus need to iterate the operation quite some times.)

How can I boost it to bigger matrices (time wise)?

Thanks!

J. M.'s missing motivation
  • 124,525
  • 11
  • 401
  • 574
Fransen
  • 31
  • 1
  • 2
  • 6
    Impossible to say without a representative example. Can even make it smaller, so long as it captures the class of matrices under consideration (numeric, symbolic, some mix,...?) – Daniel Lichtblau Oct 15 '15 at 15:57
  • @DanielLichtblau If the OP has version 10.2 this could be related to an actual bug - see the answer by aeolus below. – Jens Oct 15 '15 at 19:06

3 Answers3

8

You can in fact use Mathematica for rather large eigenvalue problems if they have floating-point entries.

Here is a band-structure calculation in one dimension where the eigensystem of a $40\times 40$ matrix is computed 81 times (at different wave numbers):

h[u_, k_, dim_] := N[
  DiagonalMatrix[(k + Range[-dim, dim])^2]
   + Table[u[j - n], {j, -dim, dim}, {n, -dim, dim}]]

bands[u_, k_, dim_, n_] := 
 Module[{hamil, max, eigenvalues, eigenvectors},
  hamil = h[u, k, dim];
  max = Norm[Flatten[hamil]];
  {eigenvalues, eigenvectors} =
   Eigensystem[hamil - max IdentityMatrix[2 dim + 1], n];
  {eigenvalues + max, eigenvectors}
  ]

Clear[potential];
potential[n_] := 1/(n^2 + 1);

l = Table[First[bands[potential, k, 40, 4]], {k, -1, 1, .025}];

ListLinePlot[Transpose[l], DataRange -> {-1, 1}, Frame -> True, 
 FrameLabel -> {"k", "ℰ"}]

bands

It doesn't take more than about a second to do this calculation.

The important thing is to wrap the matrix in a command such as N so that the entries are machine-precision numbers that can be handled by the MKL numerical library. Depending on the origin of the matrix, it could also help to wrap the matrix in SparseArray before doing the eigenvalue computation.

To describe the specific example here: potential is the Fourier amplitude of a periodic potential, and h uses it as input to construct the hermitian matrix corresponding to the Hamiltonian in Fourier space. This is then diagonalized with the command bands for a specific value of the wave number k. The dimension of the matrix returned by h can be specified in the third argument (dim) of h.

Jens
  • 97,245
  • 7
  • 213
  • 499
4

I was actually trying to solve the same sort of problem with Mathematica. As of 10.2, the Eigensystem function has a bug that will revert to using dense matrices even when fed sparse matrices.

You may want to consider using MATLAB until this is fixed. (Hopefully in 10.3!I will update when I get confirmation.)

Source: email correspondence with Wolfram Support.

Karsten7
  • 27,448
  • 5
  • 73
  • 134
aeolus
  • 41
  • 2
  • Thanks for this warning - I just answered but have version 10.1. – Jens Oct 15 '15 at 19:00
  • I actually encountered the problem in 10.1 as well, but was wrapping/using only sparse matricies 100sX100s in dimension. I never tried wrapping in N[] though... maybe worth a shot. Your implementation seems to work well and fast regardless! Up! – aeolus Oct 15 '15 at 20:04
  • Oh, OK - good to know. I heard 10.3 is just around the corner, so let's hope for a happy end. I guess I'll remove the caveat from my answer since it should work. – Jens Oct 15 '15 at 21:18
  • Is that bug fixed in the current version? – Coolwater May 04 '16 at 17:55
  • I'm not sure, I've kinda moved to Jupyter/Python since this. I'll try to check soon! – aeolus May 04 '16 at 22:21
-4

Generally working with large matrix needs large amount of RAM and Mathematica is not a good choice for this kind of problems. Using linear algebra libraries like LAPACK and MKL can be very efficient and easy. There is a good example and clear explanation here.

Abolfazl
  • 105
  • 4