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)?


  • 6
    $\begingroup$ 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,...?) $\endgroup$ – Daniel Lichtblau Oct 15 '15 at 15:57
  • $\begingroup$ @DanielLichtblau If the OP has version 10.2 this could be related to an actual bug - see the answer by aeolus below. $\endgroup$ – Jens Oct 15 '15 at 19:06

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}

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", "ℰ"}]


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.

| improve this answer | |

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.

| improve this answer | |
  • $\begingroup$ Thanks for this warning - I just answered but have version 10.1. $\endgroup$ – Jens Oct 15 '15 at 19:00
  • $\begingroup$ 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! $\endgroup$ – aeolus Oct 15 '15 at 20:04
  • $\begingroup$ 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. $\endgroup$ – Jens Oct 15 '15 at 21:18
  • $\begingroup$ Is that bug fixed in the current version? $\endgroup$ – Coolwater May 4 '16 at 17:55
  • $\begingroup$ I'm not sure, I've kinda moved to Jupyter/Python since this. I'll try to check soon! $\endgroup$ – aeolus May 4 '16 at 22:21

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.

| improve this answer | |
  • 8
    $\begingroup$ ...if the matrix has numeric entries. Of note is that Mathematica actually uses LAPACK under the hood for most linear algebraic operations. $\endgroup$ – J. M.'s discontentment Oct 15 '15 at 16:09
  • $\begingroup$ Related mathematica.stackexchange.com/a/88349/193 $\endgroup$ – Dr. belisarius Oct 15 '15 at 16:10
  • 5
    $\begingroup$ I fail to see any relevance here. For matrices of machine doubles Mathematica uses a storage format pretty much as condensed as any other, and uses Lapack and MKL libraries. This is not exactly a secret. $\endgroup$ – Daniel Lichtblau Oct 15 '15 at 20:45

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.