Matrix exponential of a skew-Hermitian matrix with fortran 95 and LAPACK

Question

I'm just getting tucked into fortran 95 for some quantum mechanics simulations. Honestly, I've been spoiled by Octave so I've taken matrix exponentiation for granted. Given a (small, $n\leq 36$) skew-Hermitian matrix of size $n\times n$, what is the most efficient way of using LAPACK to solve this problem? I'm not using the LAPACK95 wrapper, just direct calls to LAPACK.

Answer 1

Matrix exponentials of skew-Hermitian matrices are cheap to compute:

Suppose $A$ is your skew-Hermitian matrix, then $iA$ is Hermitian, and via zheevd and friends you can get the decomposition

$$iA = U \Lambda U^H,$$

where $U$ is the unitary eigenvector matrix and $\Lambda$ is real and diagonal. Then, trivially,

$$A = U (-i \Lambda) U^H.$$

Once you have $U$ and $\Lambda$, it is easy to compute

$$\exp(A)=\exp(U (-i\Lambda) U^H)=U \exp(-i\Lambda) U^H$$

by first exponentiating the eigenvalues, setting $B := U$ via zcopy, performing $B := B \exp(-i \Lambda)$ by running zscal on each column with an exponentiated eigenvalue, and finally setting your result to

$$ \exp(A) := B U^H$$

via zgemm.

Answer 2

Since I'm on my phone, I can't link things easily, and will add links later. You'll probably want to look at the paper "19 Dubious Ways to Calculate the Matrix Exponential", the Fortran library EXPOKIT, Jitse Niesen's paper on Krylov methods for calculating the Matrix exponential, and some of Nick Higham's recent papers on matrix exponentials. It's more common to need the product of a matrix exponential and a vector than the matrix exponential alone, and here, Krylov methods can be quite helpful. For smaller, dense matrices like the ones you describe, Padé methods might be better, but I've had a lot of success with Krylov methods when used inside exponential methods for numerical integration of ODEs.

Answer 3

The complex eigensolution approach is mathematically correct, but it does more work than necessary. Unfortunately, the improved approach I'm about to describe cannot be implemented with LAPACK calls.

Look at R. C. Ward and L. J. Gray, ACM Trans. Math. Soft. 4, 278, (1978). This describes the software that is available in TOMS algorithm 530, and which you can download from netlib. This describes how to factor the skew symmetric matrix $X$ as

$$X = U D U^T$$

where $U$ is real orthogonal and $D$ is real skew-symmetric and block diagonal. The diagonal subblocks are either $2\times 2$ or $1\times 1$. Because it is block diagonal, you can exponentiate each subblock separately. The $1\times 1$ blocks are zero, and $\exp(0)=1$, so those are trivial. The $2\times 2$ subblocks are done with

$$\exp \begin{pmatrix} 0 & -t \\ t & 0 \end{pmatrix} = \begin{pmatrix} \cos t & -\sin t \\ \sin t & \cos t \end{pmatrix}$$

The exponential matrix that you want is then given by

$$\exp(X) = U \exp(D) U^T$$

I have used this approach in my quantum chemistry codes for several decades and I have never had any problems with any of the software involved.

Answer 4

If all you need is the matrix exponential multiplied by a vector, then this fortran subroutine may be of some use to you. It computes:

$(e^A)v$

where v is a vector, and A is a regular hermitian matrix. It is a subroutine from the EXPOKIT library

Otherwise, you may want to consider this subroutine, which works for any general complex matrix A.