Computing the matrix exponential and the Cholesky factor of a related finite horizon Gramian

TL;DR

This work develops a scaling-and-squaring–style algorithm to simultaneously compute the matrix exponential $\Phi(A)=e^{A}$ and a Cholesky factor of the finite-horizon Gramian $G(A,B)=\int_0^1 e^{A\tau}BB^{*}e^{A^{*}\tau}\,d\tau$, without forming $G$ explicitly. It introduces a doubling recursion that propagates $\Phi$ and $G$ together, using Legendre-based initial approximations tied to diagonal Padé expansions, and computes the Gramian factor via QR after each step. A rigorous backward-error analysis guarantees unit-round-off accuracy in double precision, and a forward-error bound is derived via the Gramian’s condition number, with a subspace-preservation result ensuring the approximated Gramian shares the exact controllable subspace under appropriate scaling. The method, implemented in the Julia package FiniteHorizonGramians.jl under MIT license, demonstrates robust numerical performance on standard test problems, including nilpotent, randomly generated, and Laguerre-network A matrices, illustrating practical impact for robust state-space analysis and square-root filtering on moderate-dimension dense systems.

Abstract

In this article, an efficient numerical method for computing both the matrix exponential and a finite horizon controllability Gramian in Cholesky-factored form is proposed. The method is applicable to general dense matrices of moderate size and produces a Cholesky factor of the Gramian without computing the full product. It is a generalization of the scaling-and-squaring approach for approximating the matrix exponential and exploits a similar doubling formula for the Gramian to compute both at once. The required computational effort is thereby kept modest. Most importantly, a rigorous backward error analysis is provided, which guarantees that the approximation is accurate to the round-off error level in double precision. This accuracy is illustrated in practice on a large number of standard test examples. The method has been implemented in the Julia package FiniteHorizonGramians.jl, which is available online under the MIT license. Code for reproducing the experimental results is included in this package, as well as code for determining the optimal method parameters. The analysis can thus easily be adapted to a different finite-precision arithmetic.

Figures (4)

• Figure 1: The series $\theta_q$ and $\eta_q$ for $q=1, \ldots, 21$ (left), and the difference in the scaling parameter $s$ when selected to satisfy $\norm{A_s} \leq \theta_q$ or $\norm{A_s} \leq \eta_q$, respectively (right).
• Figure 1: The results of experiment 0. The relative error in the approximated Gramians (left) and Cholesky factors (right). The dashed line is the approximate error estimate of the computed Gramian given by \ref{['eq:erel']}.
• Figure 2: The results of experiment 1. Every value on the horizontal axis corresponds to a specific matrix $A$ from the builtin data set in MatrixDepot.jl. For each $A$, there are 50 grey dots, each corresponding to the relative error in the approximated Gramian for a single randomly generated matrix $B$. The dashed line is the approximate error estimate of the computed Gramian given by \ref{['eq:erel']}.
• Figure 3: The results of experiment 2. The relative error in the approximated Gramians, plotted against the dimension of the Laguerre network \ref{['eq:laguerre_network']}. The dashed line is the approximate error estimate of the computed Gramian given by \ref{['eq:erel']}.

Theorems & Definitions (21)

• Lemma 2.1
• Proof 1
• Proposition 3.1
• Proposition 3.2
• Proof 2
• Proposition 3.3
• Proof 3
• Theorem 3.4
• Lemma 4.1
• Proof 4