Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It only takes a minute to sign up.
Sign up to join this community
Let $A, G, Q$ be real, square, dense matrices. $G$ and $Q$ are symmetric. Let
$$H = \begin{bmatrix} A & -G \\ -Q &-A^T \end{bmatrix}$$
be a Hamiltonian matrix. I want to compute the matrix exponential of $H$. I need the full matrix exponential, $e^{tH}$, not only the matrix-vector product. Are there any specialized algorithms or libraries available to compute the exponential of a Hamiltonian matrix?
Very quick answer...
The exponential of a Hamiltonian matrix is symplectic, a property that you probably wish to preserve, otherwise you would simply use a non-structure-preserving method. Indeed, there is no real speed advantage in using structured method, just structure preservation.
A possible way to solve your problem is the following. First find a symplectic matrix such that $\hat{H}=M^{-1}HM=\begin{bmatrix} \hat{A} & -\hat{G}\\ 0 & -\hat{A}^T \end{bmatrix}$ is Hamiltonian and block upper triangular, and $\hat{A}$ has eigenvalues in the left half-plane. You get this matrix for instance by taking $\begin{bmatrix}I & 0\\ X & I\end{bmatrix}$, where $X$ solves the Riccati equation associated to $H$, or (more stable since it's orthogonal) by reordering the Schur decomposition of $H$ and applying the Laub trick (i.e., replacing the unitary Schur factor $\begin{bmatrix}U_{11} & U_{12} \\ U_{21} & U_{22}\end{bmatrix}$ with $\begin{bmatrix}U_{11} & -U_{12} \\ U_{12} & U_{11}\end{bmatrix}$). You may have trouble doing it if the Hamiltonian has eigenvalues on the imaginary axis, but that's a long story and for now I will suppose it doesn't happen in your problem.
Once you have $M$, you have $\exp(H)=M\exp(\hat{H})M^{-1}$, and you can compute $$ \exp(\hat{H}) = \begin{bmatrix} \exp(\hat{A}) & X\\ 0 & \exp(-\hat{A}^T) \end{bmatrix}, $$ where $X$ solves a certain Lyapunov equation, I believe something like $$ \hat{A} X + X \hat{A}^T = -\exp(\hat{A}) \hat{G} - \hat{G} \exp(-\hat{A}^T) $$ (signs may be wrong; impose $\exp(\hat{H})\hat{H}=\hat{H}\exp(\hat{H})$ and expand blocks to get the correct equation. Look up "Schur-Parlett method" for a reference to this trick).
Then the three factors are exactly symplectic. Just use them separately: do not compute the product or you will lose this property numerically.
You might have an option to use hierarchical matrices ($\mathcal H$-matrices) and the corresponding functionality of the libraries that support them.
Effectively, if every matrix $A$, $G$, and $Q$ are well & efficiently represented in $\mathcal H$-matrix format, then the block Hamiltonian matrix $H$ is effectively a block $\mathcal H$-matrix. The question about the representation of $A$, $G$, and $Q$ in the hierarchical form comes down to their origin: if one can find low-rank structures within them (possible row/col indices permutations applied), then this approach is viable. A plausible example would be if $A$, $G$, and $Q$ come from an integral equation which will also explain their dense structure and potential for compression (depending on the kernel).
The formal requirement for this method to work, will be if $(H-\lambda I)^{-1}$ is representable in $\mathcal H$-matrix format; however, I would start directly from constructing $\mathcal H$-matrix representation of $A$, $G$, and $Q$ and hope for the best.
Application of $\mathcal H$-matrices to matrix exponentiation is well discussed in:
There are several libraries that support $\mathcal H$-matrices. I know that the old $\mathcal H$Lib supported some form matrix exponentiation and had the necessary pieces to build one.
Downsides of this approach:
Positives:
