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?

  • 1
    $\begingroup$ 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. $\endgroup$ – Szabolcs Aug 16 '14 at 21:39
  • 1
    $\begingroup$ Make sure that your matrix contains only machine precision floating point numbers. $\endgroup$ – Szabolcs Aug 16 '14 at 21:43
  • $\begingroup$ @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. $\endgroup$ – lagoa Aug 17 '14 at 20:10
  • $\begingroup$ Yes, I have tried with Eigenvalues. Eigensystem would take a little longer, see the extended comment I wrote in an answer. $\endgroup$ – Szabolcs Aug 17 '14 at 21:10
  • $\begingroup$ Sorry, I was wrong. I'll get back to this later today. $\endgroup$ – 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.

  • $\begingroup$ 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. $\endgroup$ – Szabolcs Aug 18 '14 at 22:21
  • $\begingroup$ @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. $\endgroup$ – 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.

  • $\begingroup$ 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? $\endgroup$ – lagoa Aug 17 '14 at 20:43
  • $\begingroup$ 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. $\endgroup$ – Szabolcs Aug 17 '14 at 21:10
  • $\begingroup$ 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 ? $\endgroup$ – coincoin Aug 17 '14 at 21:15

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.