Efficient scaling and squaring method for the matrix exponential

TL;DR

This work tackles the efficient computation of the matrix exponential $\,e^A\,$ to a user-specified tolerance by an adaptive algorithm that combines scaling and squaring with a curated set of polynomial and rational approximants, including Taylor, Padé, and diagonal Padé forms. Backward-error analysis guides the selection of approximation and the scaling parameter, while a cost-oriented framework decomposes Padé forms into fractions to reduce multiplications and enable parallelization. The authors provide an adaptive LUT and algorithm to minimize total cost across tolerances, and demonstrate significant computational savings over state-of-the-art software, while preserving Lie-group structures through diagonal Padé methods when desired. The approach is validated on diverse numerical tests, including random diagonally dominant, symplectic, and unitary matrices, highlighting practical impact for exponential integrators and structure-preserving computations.

Abstract

This work presents a new algorithm to compute the matrix exponential within a given tolerance. Combined with the scaling and squaring procedure, the algorithm incorporates Taylor, partitioned and classical Padé methods shown to be superior in performance to the approximants used in state-of-the-art software. The algorithm computes matrix--matrix products and also matrix inverses, but it can be implemented to avoid the computation of inverses, making it convenient for some problems. If the matrix A belongs to a Lie algebra, then exp(A) belongs to its associated Lie group, being a property which is preserved by diagonal Padé approximants, and the algorithm has another option to use only these. Numerical experiments show the superior performance with respect to state-of-the-art implementations.

Figures (4)

• Figure 1: Comparison of tolerance (dashed) and practical average cost plotted vs. total computational cost of superdiagonal methods. The thin vertical dotted line represents the cost in the current software that implements higham09tsa. Note that the distance between the dotted vertical line and the solid one represents cost savings.
• Figure 2: Comparison of tolerance (dashed) and practical average cost plotted vs. total computational cost of diagonal methods. The thin vertical dotted line represents the cost in the current software that implements higham09tsa. Note that the distance between the dotted vertical line and the solid one represents cost savings.
• Figure 3: Comparison of symplectic property preservation by the algorithm using only diagonal Padé methods for matrices \ref{['eq:sympl_err_best_case']} (left) and \ref{['eq:sympl_err_rand']} (right). Reference (dashed) is the LinearAlgebra.exp function from Julia 1.10.2. Machine $\varepsilon$ is added to avoid singularities.
• Figure 4: Comparison of unitarity preservation (where $J=I$) by the algorithm using only diagonal Padé methods for matrices \ref{['eq:sympl_err_best_case']} (left) and a random complex matrix (right). Reference (dashed) is the LinearAlgebra.exp function from Julia 1.10.2. Machine $\varepsilon$ is added to avoid singularities.