0
$\begingroup$

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?

Improve this question
$\endgroup$
5
  • 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, 2014 at 21:39
  • 1
    $\begingroup$ Make sure that your matrix contains only machine precision floating point numbers. $\endgroup$
    – Szabolcs
    Aug 16, 2014 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, 2014 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, 2014 at 21:10
  • $\begingroup$ Sorry, I was wrong. I'll get back to this later today. $\endgroup$
    – Szabolcs
    Aug 19, 2014 at 12:22

2 Answers 2

Reset to default
1
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 1
    $\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, 2014 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, 2014 at 22:45
0
$\begingroup$

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.

Improve this answer
$\endgroup$
3
  • $\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, 2014 at 20:43
  • 1
    $\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, 2014 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, 2014 at 21:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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