We present a new method for computing the action of the matrix exponential on a vector, $e^{At}v$. The proposed approach efficiently evaluates the solution for all $t$ within a prescribed bounded interval by expanding it into an orthogonal polynomial series. This method is derived from a new representation of the matrix exponential in the so-called $\star$-algebra, an algebra of bivariate distributions. The resulting formulation leads to a linear system equivalent to a matrix equation of Stein type, which can be solved by either direct or Krylov subspace methods. Numerical experiments demonstrate the accuracy and efficiency of the proposed approach in comparison with state-of-the-art techniques.
This paper introduces a new numerical method for evaluating the action of the matrix exponential, e tA v, for every t in a given finite interval, where A is a matrix and v a vector. The approach exploits the fact that the solution of a linear autonomous system of ordinary differential equations (ODEs) can be expressed by a matrix exponential. In turn, the new numerical method is obtained by exploiting the new formula for the solution of such an ODE by the ⋆-product [5,17] (a generalization of the Volterra convolution) followed by an appropriate discretization of the resulting expression [14,15]. The computation of matrix exponentials is a fundamental task in many areas, including quantum dynamics, control theory, network science, and the time integration of large-scale linear ODEs arising from spatial discretizations of PDEs (see, e.g., [10,18,9,7]).
Among the most widely used approaches are the ones based on Krylov subspace methods (e.g., [11]). The expv routine in Expokit [19] employs the Arnoldi or Lanczos process to construct a Krylov basis and approximates e tA v via projection, requiring only sparse matrix-vector products and avoiding explicit construction of e tA . In [2], expmv_tspan returns either a single matrix e tA B or a sequence of matrices e t k A B evaluated on an equally spaced grid of points t k . In the latter case, it uses the scaling phase of the scaling-and-squaring method combined with a truncated Taylor series approximation. The scaling factor and Taylor series degree are determined using the sharp truncation error bounds from [1], which are expressed in terms of A k 1/k for selected values of k, with the norms estimated via a matrix norm estimator.
The approach proposed in this work builds on the recently developed theory of ⋆-algebra [12,16,17]. By discretizing the ⋆-product via Legendre polynomial expansions, it reformulates the problem of matrix approximation into a linear algebra problem, which can be equivalently expressed in terms of a Stein-type matrix equation. The ⋆-approach gives the solution in terms of the coefficient of the truncated Legendre polynomial expansion of e tA v in a given interval I. Therefore, the solution can be computed for every t ∈ I. This paper also provides a bound for the truncation error of the method. Combined with an Arnoldi-based Krylov projection of dimension k, the method yields an efficient algorithm for large sparse matrices, allowing the solution to be evaluated at arbitrary times without recomputation.
Numerical experiments show that the proposed method can achieve accuracy comparable to that of state-of-the-art approaches. In particular, the Arnoldi-accelerated variant delivers competitive runtimes across a broad class of test problems, making it a promising alternative for large-scale exponential integrators, especially when solutions are required at many or a priori unknown time points.
The remainder of the paper is organized as follows. Section 2 presents the mathematical formulation of the ⋆-method, deriving the truncated Legendre expansion and the corresponding Stein equation. In Section 3, we derive an upper bound for the truncation error. Section 4 describes the numerical implementation of the ⋆-method, including its integration with a Krylov subspace approach, and reports detailed numerical experiments comparing it with expv and expmv_tspan. Section 5 concludes with a discussion of potential future research.
Consider the system of linear autonomous ODEs of the form
where à ∈ C N×N is a matrix and Ũs (t) denotes the matrix-valued solution parameterized by the initial time s. Note that we set the interval I = [-1, 1] for later convenience, but this can easily be generalized to any other bounded interval. Introducing the notation Ũs (t) := U(t, s), for t ≥ s the ODE can be expressed equally in the ⋆-product framework [15,5,12] as
where the bivariate matrix operator A(t, s) is defined by
where Θ(ts) is known as the Heaviside step function. Trivially, for t ≥ s we get Θ(ts) = 1 and the solution of Eq. ( 2) is the matrix exponential.
Following [15], we introduce a spectral approximation framework to numerically solve Eq. ( 2). Let {p k (t)} M-1 k=0 denote the Legendre polynomials rescaled and shifted so that they are orthonormal on the interval of interest
with δ kℓ the Kronecker delta. Since these polynomials form a basis for analytic functions over I, for a given positive integer M, we can approximate each entry of the solution U(t, s) by a truncated double expansion.
where the Legendre-Fourier coefficients are given by
For every pair (i, j), these coefficients are collected into the M × M matrix (block)
Hence, defining the vector
, the Legendre expansion of the solution element U i, j (t, s) can be rewritten in the matrix form
As a consequence, we can approximate the solution of Eq. ( 2) using the Kronecker product ⊗ by the formula
where I N is the identity matrix of size N, and U M is named the coefficient matrix of U
This content is AI-processed based on open access ArXiv data.