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