# Most efficient library to diagonalize exactly large hermitian or unitary matrices

I am working on a physics problem which requires obtaining the exact eigenvalues and eigenvectors of Hermitian and Unitary matrices numerically. Naturally I would like to ask the experts what are the most efficient libraries to diagonalize the referred matrices. In particular, I would like to diagonalize efficiently matrices of size 12870*12870. I have been using Mathematica for convenience, but I am open to use any language you might suggest. I also have access to a HPC cluster in which the diagonalization of a unitary matrix of size 12870*12870 takes around 7000 seconds in Mathematica using 12 cores. Can this time be improved? Note: The matrices are numeric (not symbolic).

To be more precise, I diagonalized a random Hermitian matrix and a random Unitary matrix of size 12870 on the cluster using Mathematica. Here is the code I used:

With[{mat = RandomReal[1, 12870 {1, 1}]}, AbsoluteTiming[{vals, vecs} = Eigensystem[mat+Transpose[mat]];]]

With[{matu = Transpose[vecs].(Exp[I*vals]*Conjugate[vecs])}, AbsoluteTiming[Eigensystem[matu];]]


The results were:

In[1]:= Out[1]= {265.369827, Null}
In[2]:= Out[2]= {7042.662152, Null}


So the unitary matrix diagonalization takes a lot more time. Can this be improved using other languages?

• Something's wrong with the timings you cite. If you have a numerical matrix of that size, finding all eigenvalues with Mathematica should not take more than 5-10 minutes on a laptop. Post the specific Mathematica code you used. – Szabolcs Aug 16 '14 at 21:39
• Make sure that your matrix contains only machine precision floating point numbers. – Szabolcs Aug 16 '14 at 21:43
• @Szabolcs I acted with N function on the unitary matrix before applying Eigensystem to it. Are you sure it should take 5-10 minutes? Have you tried an example? Thanks for the suggestion. – lagoa Aug 17 '14 at 20:10
• Yes, I have tried with Eigenvalues. Eigensystem would take a little longer, see the extended comment I wrote in an answer. – Szabolcs Aug 17 '14 at 21:10
• Sorry, I was wrong. I'll get back to this later today. – Szabolcs Aug 19 '14 at 12:22

If you want to diagonalize significantly larger matrices on large HPC clusters, I would look at Elemental, which is a state-of-the-art distributed dense linear algebra library. It's a better alternative than other libraries out there (e.g., ScaLAPACK, PLAPACK), and builds upon BLAS and LAPACK implementations.

• This library seems to be aimed at distributed memory computation. If the whole matrix fits in memory, how would its performance compare to e.g. the MKL running on a shared memory system? I'm asking because the OP's matrix likely fits in memory (about 1.3 GB for one matrix), and I always thought that switching to C from MATLAB/Mathematica will not help at all when the whole computation is a single function call to a linear algebra operation with a large matrix. The bottleneck is the MKL or other LAPACK implementation anyway, not the language it's being calling from. – Szabolcs Aug 18 '14 at 22:21
• @Szabolcs: It is aimed at distributed memory computation, and was mainly a response to the answer of coincoin. If the matrix fits in memory, and you run it in serial, I would think that yes, LAPACK would be the bottleneck, not the language used to call it. You could still run this computation in parallel on an architecture like BG/P and get a speedup relative to a serial implementation; see the case studies in dl.acm.org/citation.cfm?doid=2427023.2427030, which include an eigendecomposition of a 10,000 by 10,000 matrix on 8192 cores. – Geoff Oxberry Aug 18 '14 at 22:45

If you have access to a HPC cluster and you want to work on matrix of larger size (12k^2 order is not much), I would suggest to use Lapack routines. This is a base/reference in numerical linear algebra available in Fortran or C.

• Could you write a code example for reference, where you generate a random unitary matrix U( U=exp(i*H), where H is a random Hermitian matrix) and diagonalize it using Lapack in C++? How much time only in diagonalization does it take for a matrix of size 12870? – lagoa Aug 17 '14 at 20:43
• Mathematica, like all other similar software (MATLAB, etc.) also uses LAPACK (actually it uses the MKL). For these types of problems (large matrices, a single call to eigenvalue computation) there is typically no advantage to using C or Fortran. The limiting factor is going to be the LAPACK implementation anyway. – Szabolcs Aug 17 '14 at 21:10
• I thought there was still an overhead when using these softwares. And do these softwares implement parallel algorithms for a multi core environment since the original lapack does not provide that. And if he uses a HPC cluster maybe he may needsto get a code for production or Mathematica provides tools when using parallel architectures ? – coincoin Aug 17 '14 at 21:15