Performance of Hermitian eigensolvers

General Usage Performance
1

I usually work with full diagonalisation of Hermitian matrices. After some research I realised that there are two options to retrieve the eigenvalues and eigenvectors of such a matrix.
LinearAlgebra.LAPACK.syev!() and eigen!(). From the source code of eigen!() it seems that it is calling a general eigensolver routine called ggevx().

The LinearAlgebra.LAPACK.syev!() routine required some digging into the source code to figure out that it is also intended for Hermitian matrices (docs.julialang.org only talks about symmetric matrices). I realised that it is calling specialised Hermitian routines in LAPACK.

Now after building up a 2000x2000 Hermitian matrix I called both routines the following way and timed their peformance:

H = weyl(10, 10, 10)
Ht = deepcopy(H)

@time eigsys = LAPACK.syev!('V', 'U', H)
@time eigsyst = eigen!(Hermitian(Ht))

To my surprise the eigen! routine is quite a bit faster:

  8.998029 seconds (6 allocations: 1.068 MiB)
  6.128692 seconds (16 allocations: 123.658 MiB, 1.71% gc time)

To me this is quite counterintuitive behaviour. I would have guessed that the specialised routines should work faster.

2

If you look more closely at LinearAlgebra/src/symmetric.jl you find that eigen!(::Hermitian) calls the syevr routines from LAPACK (note the trailing r). For full diagonalization these use Dhillon and Parlett’s “Relatively Robust Representations” which are more efficient than the traditional syev scheme.

3 Likes