Compiling function solving for the smallest eigenvalue of a matrix

Question

I am trying to speed up my code, which many times finds a lowest eigenvalue of a Hermitian (symmetric, real) matrix. Using a small 2x2 matrix as an example, this works:

Energy = Compile[{{matrix, _Real, 2}}, Min[Re[Eigenvalues[matrix]]], {{Eigenvalues[_], _Complex, 1}}];
Energy[{{1, -2}, {-2, 3}}]
-0.236068

I have to use {{Eigenvalues[_], _Complex, 1}} here, otherwise Mathematica skips using the compiled function because of a type mismatch. Function Eigenvalues by default returns a vector of complex eigenvalues.

It should be in principle possible to compile a faster version, since I am only interested in the smallest eigenvalue, using Eigenvalues[matrix, -1], which only solves for one eigenvalue. However, I can't seem to make it to work:

Energy = Compile[{{matrix, _Real, 2}},Eigenvalues[matrix, -1][[1]], {{Eigenvalues[_], _Complex, 1}}];
Energy[{{1, -2}, {-2, 3}}]
During evaluation of CompiledFunction::cfex: Could not complete external evaluation at instruction 1; proceeding with uncompiled evaluation. >>
2 - Sqrt[5]

Where am I making a mistake? Also, is there a room for possible speedup if I somehow convince Mathematica, that I am only using Hermitian matrices and therefore expect only real eigenvalues?

Answer 1

Look at the output of CompilePrint[Energy]:

Needs["CompiledFunctionTools`"];
CompilePrint[Energy]


        1 argument
        1 Real register
        3 Tensor registers
        Underflow checking off
        Overflow checking off
        Integer overflow checking on
        RuntimeAttributes -> {}

        T(R2)0 = A1
        Result = R0

1   T(C1)1 = MainEvaluate[ Hold[Eigenvalues][ T(R2)0]]
2   T(R1)2 = Re[ T(C1)1]
3   R0 = MinRT[ T(R1)2]]
4   Return

The line MainEvaluate[ Hold[Eigenvalues][ T(R2)0]] means that the compiled program still needs to ask the main kernel to compute the Eigenvalues. This means that you gain nothing from compiling here at all. Looking at the list of compilable functions we confirm that Eigenvalues is not there.

Solution

The speediest way I know of finding eigenvalues in Mathematica is simply to give Eigenvalues a machine-number matrix:

Eigenvalues[N @ {{1, -2}, {-2, 3}}, {-1}]

This takes no measurable time on my laptop, while

AbsoluteTiming[Eigenvalues[{{1, -2}, {-2, 3}}, {-1}]]

shows 0.022 seconds. The reason is that Eigenvalues works differently for exact and numerical matrices. For a matrix with exact entries (like {{1, -2}, {-2, 3}}) Eigenvalues will interpolate the characteristic polynomial of the matrix and look for roots of that polynomial. When given a matrix with machine numbers on the other hand, Mathematica uses the very fast ARPACK library, which uses Arnoldi methods to find eigenvalues. I'm not too familiar with ARPACK, but I'm guessing it also specializes to the Lanczos algorithm when the input is Hermitian (see here for more info, e.g. Ctrl+F Eigenvalues)

Give this a go, and if you still need more speed, there are probably other places in the code that is slowing it down. Then we would need more details in order to help you ;)