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?

  • 2
    $\begingroup$ Do you want the matrix exponential itself, or do you really just want to solve the ODE $\dot z = Hz$? $\endgroup$ – Daniel Shapero Apr 17 '18 at 15:30
  • $\begingroup$ I Need The Matrix Exponential itself. But equivalently i can solve the ODE $\dot{Z} = H Z,\ Z(0)=I$. $\endgroup$ – Max Behr Apr 17 '18 at 15:43
  • 2
    $\begingroup$ Benner's structure preserving eigensolvers can deal with the similarity transformation to ease the matrix exponential computaiton. $\endgroup$ – percusse Apr 17 '18 at 21:00
  • 1
    $\begingroup$ @RichardZhang The brutal way is the QZ decomposition. Check for example starting from link.springer.com/article/10.1007/s002110050315 for more details. $\endgroup$ – percusse Apr 17 '18 at 23:14
  • 2
    $\begingroup$ The paper 19 Dubious Ways to Compute the Exponential of a Matrix, 25 Years Later covers many bad (and a few good) ways to compute the matrix exponential. It isn't specific to Hamiltonian problems but nonetheless is really valuable if you're working on these kinds of problems. $\endgroup$ – Daniel Shapero Apr 18 '18 at 15:43

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.

  • 1
    $\begingroup$ Currently i am doing it in a slightly different way. I use the stabilizing spd solution of the ARE to setup the similiarity transformation for $H$ and obtain $\tilde{H}= \begin{bmatrix} \hat{A} & -G \\ 0 & -\hat{A}^T \end{bmatrix} $ like in your suggestion. Then let $X_L$ the solution of the Lyap.eqn $\hat{A} X_L + X_L \hat{A}^T = -G$ and setup a second similarity transf. $M_2=\begin{bmatrix} I & X_L \\ 0 & I \end{bmatrix}$. Applying this to $\hat{H}$ and obtain $\hat{\hat{H}}$, which is block diagonal with $\hat{A}$ and $\hat{-A^T}$ as blocks $\endgroup$ – Max Behr May 18 '18 at 15:02

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:

  • relies on the efficient representation of $A$,$G$, and $Q$
  • does not take advantage of the Hamiltonian structure


  • the compressed representation of the matrix exponential, though it is still a matrix, not just a way to do an MVP
  • linear-logarithmic complexity (provided the low-rank assumption is there)
  • the library might take advantage of the transposition and symmetry in the blocks

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.