Reusing eigendecomposition for matrix exponentials

Looking into the documentation, exp is using the eigen decomposition to calculate the exponetial exp(A) of a matrix A. Is there a function in Julia/BLAS to compute exp(c*A) for a constant, possible imaginary, c for a lot of different c but the same A, reusing the eigen decomposition of A?

Only for symmetric/Hermitian matrices. For general dense matrices, it is scaling and squaring.

1 Like

Ah yes, sure. I meant A to be symmetric/hermitian.

I think it would be quite easy to implement exp(::LinearAlgebra.Eigen), so you could just compute eigen(A) once, scale the eigenvalues and compute the exponential of those. That would probably make quite a straightforward PR, if you want to give it a go.

Yes, I was just a bit curious if the implementation does some smart tricks, but it is implemented for hermitian matrices as

F = eigen(A)
Symmetric((F.vectors * Diagonal(($func).(F.values))) * F.vectors')

Where func are the normal candidates for functions implemented in a standard linear algebra library. The most general use case would be achieved by having a function applying an (arbitrary) function to the matrix A. Would something like

applyfunction(func, eigendecomp) = Symmetric((eigendecomp.vectors * Diagonal((func).(eigendecomp.values))) * eigendecomp.vectors')

make sense? I do not know how sensible the compiler is to passing functions as arguments.

Should I open an issue on github and move the discussion there?

This should do it:

exp_scaled(c::Number, F::Eigen) = F.vectors * Diagonal(exp.(c .* F.values)) * F.vectors'

I’m not sure a library function is really needed for something so simple…

2 Likes

The compiler is very good about higher-order functions, and can even inline the function if you add a type parameter.

2 Likes