Computing the action of a matrix exponential on an interval via the $\star$-product approach

1. Problem Statement

The paper addresses the computational challenge of evaluating the action of a matrix exponential on a vector, $e^{At}v$, for all time points $t$ within a prescribed bounded interval $I = [t_0, t_{end}]$.

• Context: This problem arises frequently in quantum dynamics, control theory, network science, and the time integration of large-scale linear Ordinary Differential Equations (ODEs) derived from spatial discretizations of Partial Differential Equations (PDEs).
• Limitations of Existing Methods:
  • Krylov Subspace Methods (e.g., expv): Efficient for computing $e^{At}v$ at a single time point or a few points but require recomputation or projection for every new time instance.
  • Scaling-and-Squaring with Taylor Series (e.g., expmv_tspan): Can compute solutions on a grid of points but typically requires explicit evaluation at specific, pre-defined time steps.
  • Goal: Develop a method that provides a continuous representation of the solution over the entire interval, allowing evaluation at arbitrary time points without re-computation, while maintaining high accuracy and efficiency for large sparse matrices.

2. Methodology: The ⋆-Product Approach

The proposed method, termed the ⋆-method, is based on a novel algebraic framework involving the ⋆-algebra of bivariate distributions.

A. Mathematical Formulation

1. Reformulation as an ODE: The problem is cast as solving a linear autonomous ODE system $\frac{d}{dt}\tilde{u}(t) = \frac{1}{2}A\tilde{u}(t)$ with initial condition $\tilde{u}(-1) = v$ on the interval $I=[-1, 1]$.
2. The ⋆-Algebra: The solution is expressed using the ⋆-product (a generalization of the Volterra convolution) in the space of bivariate distributions. The matrix exponential is represented as a Neumann series or Taylor series within this algebra:
   $$\exp(A(t-s)) = (\delta I - A\Theta)^{-\star} \star \Theta$$
   where $\Theta$ is the Heaviside step function and $\delta$ is the Dirac delta.
3. Spectral Discretization: The bivariate functions are approximated using Legendre polynomials. The solution $U(t, s)$ is expanded into a truncated double series of Legendre polynomials.
4. Reduction to Linear Algebra:
   • The infinite-dimensional ⋆-algebra problem is discretized into a finite-dimensional linear algebra problem.
   • The approximation leads to a Stein-type matrix equation (also known as a discrete Lyapunov equation):
     $$X - \frac{1}{2} T_M X \tilde{A}^T = \phi_M(-1)v^T$$
     where $T_M$ is the coefficient matrix of the Heaviside function expanded in Legendre polynomials, and $\phi_M$ is the vector of basis polynomials.
   • Solving this equation yields a coefficient matrix $X$. The solution at any time $t$ is then recovered via:
     $$e^{At}v \approx X^T T_M^T \phi_M(t)$$

B. Implementation Strategies

The paper proposes two implementation variants:

1. Direct ⋆-Method: Solves the Stein equation directly using standard decomposition methods (e.g., Schur decomposition via MATLAB's dlyap). Suitable for moderate matrix sizes.
2. Arnoldi-Accelerated ⋆-Method: For large sparse matrices, the method projects the problem onto a Krylov subspace using the Arnoldi algorithm.
   • It constructs an orthonormal basis $V_k$ for the Krylov subspace $\mathcal{K}_k(A, v)$.
   • The Stein equation is projected onto this subspace, reducing the problem size significantly.
   • This allows the method to scale to large sparse systems while retaining the ability to evaluate the solution at arbitrary times.

3. Key Contributions

1. Continuous Interval Solution: Unlike standard methods that compute discrete points, this approach provides a global polynomial approximation valid for the entire interval, enabling evaluation at any $t \in I$ without re-solving.
2. Theoretical Error Bounds: The authors derive a rigorous upper bound for the truncation error (Theorem 3.5). They prove that the error decays super-exponentially with respect to the truncation parameter $M$ (the number of Legendre polynomials), provided certain spectral conditions are met.
3. Krylov Integration: The successful integration of the ⋆-method with Krylov subspace techniques makes it applicable to large-scale sparse matrices, a common requirement in PDE applications.
4. Stein Equation Formulation: The transformation of the matrix exponential action into a Stein matrix equation offers a new algebraic perspective that facilitates the use of established linear solvers and preconditioning techniques.

4. Numerical Results

The authors tested the method against state-of-the-art routines (expv and expmv_tspan) on a diverse set of test matrices, including:

• 2D Poisson matrices (finite difference discretization).
• Complex tridiagonal and Toeplitz matrices.
• Matrices with decaying eigenvalues and Chebyshev nodes.
• Finite Element Method (FEM) semidiscretizations of the Schrödinger equation.

Key Findings:

• Accuracy: The ⋆-method (both direct and Arnoldi-accelerated) achieves relative errors comparable to expv and expmv_tspan (typically $10^{-10}$ to $10^{-15}$).
• Efficiency:
  • For small to medium matrices, the direct method is accurate but computationally expensive due to the $O(N^3)$ cost of solving the Stein equation.
  • The Arnoldi-accelerated variant is highly competitive. It often matches or exceeds the runtime of expmv_tspan while providing a continuous solution.
  • For large matrices (e.g., $N \approx 3400$), the Arnoldi-⋆ method is significantly faster than computing solutions at many discrete points using traditional methods.
• Flexibility: The method allows for post-hoc evaluation at any time point, which is a distinct advantage over grid-based methods.

5. Significance and Conclusion

This paper introduces a paradigm shift in computing matrix exponentials by moving from point-wise evaluation to interval-based spectral approximation.

• Practical Impact: The method is particularly valuable for applications requiring solutions at many time points, or where the time points are unknown a priori (e.g., adaptive time-stepping or real-time control).
• Theoretical Foundation: It establishes a strong link between distribution theory (⋆-algebra), orthogonal polynomials, and numerical linear algebra (Stein equations).
• Future Work: The authors note that while the current error bounds are theoretically sound, they are not yet tight enough for all practical cases. Future research will focus on refining these bounds and developing specialized preconditioners for the resulting linear systems to further enhance performance.

In summary, the ⋆-method offers a robust, accurate, and flexible alternative for computing $e^{At}v$ over intervals, effectively bridging the gap between high-accuracy spectral methods and the scalability of Krylov subspace techniques.