Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Computation of the exponential function of matrices by a formula without oscillatory integrals on infinite intervals

Masato Suzuki, Ken'ichiro Tanaka

TL;DR

This work develops a quadrature-based method for computing the matrix exponential $\exp(A)$ that avoids oscillatory integrals on infinite intervals by using a non-oscillatory $I_{\alpha}(A)$ and an oscillatory $J_{\alpha}(A)$ on a finite interval. The scalar integrals are approximated with a double-exponential transform for $I_{\alpha}$ and Gauss–Legendre for $J_{\alpha}$, with rigorous, explicit error bounds carried over to the matrix setting. The authors provide theoretical guarantees and extensive numerical experiments showing favorable performance, especially when the eigenvalues have large imaginary parts, outperforming several established quadrature-based methods. The approach is well-suited for large-scale, structured matrices due to parallelizability and stable numerical behavior, offering a practical tool for applications in ODEs, Markov chains, and network analysis.

Abstract

We propose a quadrature-based formula for computing the exponential function of matrices with a non-oscillatory integral on an infinite interval and an oscillatory integral on a finite interval. In the literature, existing quadrature-based formulas are based on the inverse Laplace transform or the Fourier transform. We show these expressions are essentially equivalent in terms of complex integrals and choose the former as a starting point to reduce computational cost. By choosing a simple integral path, we derive an integral expression mentioned above. Then, we can easily apply the double-exponential formula and the Gauss-Legendre formula, which have rigorous error bounds. As numerical experiments show, the proposed formula outperforms the existing formulas when the imaginary parts of the eigenvalues of matrices have large absolute values.

Computation of the exponential function of matrices by a formula without oscillatory integrals on infinite intervals

TL;DR

This work develops a quadrature-based method for computing the matrix exponential that avoids oscillatory integrals on infinite intervals by using a non-oscillatory and an oscillatory on a finite interval. The scalar integrals are approximated with a double-exponential transform for and Gauss–Legendre for , with rigorous, explicit error bounds carried over to the matrix setting. The authors provide theoretical guarantees and extensive numerical experiments showing favorable performance, especially when the eigenvalues have large imaginary parts, outperforming several established quadrature-based methods. The approach is well-suited for large-scale, structured matrices due to parallelizability and stable numerical behavior, offering a practical tool for applications in ODEs, Markov chains, and network analysis.

Abstract

We propose a quadrature-based formula for computing the exponential function of matrices with a non-oscillatory integral on an infinite interval and an oscillatory integral on a finite interval. In the literature, existing quadrature-based formulas are based on the inverse Laplace transform or the Fourier transform. We show these expressions are essentially equivalent in terms of complex integrals and choose the former as a starting point to reduce computational cost. By choosing a simple integral path, we derive an integral expression mentioned above. Then, we can easily apply the double-exponential formula and the Gauss-Legendre formula, which have rigorous error bounds. As numerical experiments show, the proposed formula outperforms the existing formulas when the imaginary parts of the eigenvalues of matrices have large absolute values.

Paper Structure

This paper contains 20 sections, 13 theorems, 105 equations, 4 figures, 1 table.

Table of Contents

  1. Introduction
  2. Integral formula for the exponential function
  3. Formula with a complex contour integral
  4. Derivation of a formula by taking a limit
  5. Approximation formula for the exponential function by numerical integration
  6. Approximation formula for $I_{\alpha}(z)$
  7. Approximation formula for $J_{\alpha}(z)$
  8. Choice of $\alpha$
  9. Approximation formula for the exponential function of matrices
  10. Numerical experiments
  11. Convergence of the Proposed and Existing Methods
  12. Setting of the Parameters $\alpha$ and $k$
  13. Comparison of the DE and Gauss-Laguerre Quadrature
  14. Proofs
  15. Proof of Theorem \ref{['thm:total_error_of_I_alpha_n_h']}
  16. ...and 5 more sections

Key Result

Theorem 2.1

Let $z$ be a complex number with $\mathop{\mathrm{Re}} z < 0$. Then, for any real number $\alpha$ with $|\mathop{\mathrm{Im}} z| < \alpha$, the following equality holds: In the special case $z \in \mathbb{R}$, this formula is simplified as follows:

Figures (4)

  • Figure 1: The horizontal axis shows the total number of resolvent calculations, and the vertical axis shows the absolute error in the matrix 2-norm. Note that the scale of the horizontal axis is different for each figure.
  • Figure 2: The horizontal axis shows the total number of resolvent calculations (i.e., $4n+N+2=(4+k)n+2$), and the vertical axis shows the absolute error in the matrix 2-norm. The yellow line in the figure "$k=4$" is the parameter our method generates, and is used in Section \ref{['sec:numerical_experiment_1']}.
  • Figure 3: The horizontal axis shows the total number of resolvent calculations, and the vertical axis shows the absolute error of $I_{\alpha}(A_i)$ in the matrix 2-norm. The three figures on the top left, top right, and bottom left show the case of $A_1,\, A_2,$ and $A_3$ respectively, where the parameters were chosen by our proposed method. In the bottom right experiment, we used a larger $d$ for $A_3$.
  • Figure 4: Explanation of the fact that the minimum value in the right-hand side of \ref{['eq:dist_between_sp_and_H_d']} is given by \ref{['eq:formula_of_dist_between_sp_and_H_d']}. The minimum value is the minimum distance between the two singular points (i.e., $-u + \mathrm{i} (\alpha - |v|)$ and $-u - \mathrm{i} (\alpha + |v|)$) and the boundary of the region $\mathcal{H}_{d}$. Candidates of the distance are A to D. First, it clearly holds that A $\leq$ B, C $\leq$ D. Since $\mathcal{H}_d$ is symmetric w.r.t. the real axis, we have A $\leq$ C. Therefore the minimum distance is given by A.

Theorems & Definitions (33)

  • Theorem 2.1
  • Remark 2.2
  • Remark 3.1
  • Theorem 3.2
  • Theorem 3.3
  • Remark 3.4
  • Remark 3.5
  • Theorem 4.1
  • Theorem 4.2
  • Remark 4.3
  • ...and 23 more