A New Algorithm for Computing the Exponential of a Block Triangular Matrix

TL;DR

This paper introduces a specialized algorithm for computing the exponential of block triangular matrices by exploiting their structure to simultaneously obtain $e^A$, $e^B$, and $\mathcal{D}_{\exp}(A,B,E)$ without forming the full block matrix. It defines a linear operator $\mathcal{D}_f(A,B,E)$ that generalizes the Fréchet derivative, proves its key algebraic properties, and shows how to extend efficient Padé-based scaling-and-squaring schemes to compute $e^A$, $e^B$, and $\mathcal{D}_{\exp}(A,B,E)$ with a backward-error–guided selection of the scaling parameter. The authors provide a rigorous backward error analysis, derive sharp bounds, and establish that the scaling parameter can be chosen independent of $E$, with Schur decomposition employed to stabilize the squaring phase when needed. Numerical experiments demonstrate that the proposed method outperforms existing approaches in accuracy and efficiency, and they discuss diverse applications in exponential integrators, Hamiltonian systems, and polynomial-diffusion pricing. The work also sketches future directions, including expressing bounds in terms of $\alpha_p(A)$, allowing independent scaling of $A$ and $B$, and extending the framework to compute general matrix functions via Schur-based techniques.

Abstract

The exponential of block triangular matrices arises in a wide range of scientific computing applications, including exponential integrators for solving systems of ordinary differential equations, Hamiltonian systems in control theory, sensitivity analysis, and option pricing in finance. We propose a novel algorithm exploiting the block triangular structure for simultaneously computing the exponentials of the diagonal blocks and the off-diagonal block of the matrix exponential without direct involvement of the full block matrix in the computations. This approach generalizes the work of Al-Mohy and Higham on the Fréchet derivative of the matrix exponential. The generalization is established through a linear operator framework, facilitating efficient evaluation schemes and rigorous backward error analysis. The algorithm employs the scaling and squaring method using diagonal Padé approximants with algorithmic parameters selected based on the backward error analysis. A key feature is that the selection of the scaling parameter relies solely on the maximal norm of the diagonal blocks with no dependence on the norm of the off-diagonal block. Numerical experiments confirm that the proposed algorithm consistently outperforms existing algorithms in both accuracy and efficiency, making it a preferred choice for computing the matrix exponential of block triangular matrices.

Figures (3)

• Figure 1: Relative forward errors for the computation of $\mathcal{D}_{\exp}(A,B,E)$ by Algorithm \ref{['alg.expmF']}, Kenney-Laub's algorithm (KL), and expm. The solid line represents the condition number, $\mathrm{cond}_{\mathcal{D}_{\exp}}(A,B,E)$, multiplied by $u=2^{-53}$.
• Figure 2: The data in Figure \ref{['fig.test1']} presented as a performance profile.
• Figure 3: Relative errors for $\widehat{D}(\alpha) \approx \alpha \widehat{D}(1)$ with $\alpha = 1/\|E\|_1$ using Algorithm \ref{['alg.expmF']} and the KL algorithm. Data sorted as in Figure \ref{['fig.test1']}.

Theorems & Definitions (6)

• Definition 1.1
• Lemma 1.2
• Proof 1
• Theorem 3.1
• Theorem 3.2
• Proof 2