An error control framework for computing the exponential of matrices arising from the finite element discretization

1. Problem Statement

The paper addresses the challenge of computing the action of the matrix exponential, $e^A b$, where the matrix $A$ has the specific structure $A = \tau M^{-1}K$. This structure arises frequently in the finite element discretization of time-dependent partial differential equations (PDEs), such as advection-diffusion equations, where:

• $\tau > 0$ is a time step.
• $M$ is a well-conditioned, symmetric positive definite (SPD) mass matrix.
• $K$ is the stiffness/advection matrix.

The Core Difficulty:
Existing error control frameworks for matrix exponentials rely on the numerical range (field of values), $W(A) = \{v^*Av : \|v\|_2=1\}$, to bound the approximation error. Specifically, if a rational approximation $r(z) \approx e^z$ is constructed on a region containing $W(A)$, the error $\|r(A) - e^A\|_2$ can be bounded.

However, applying this to $A = \tau M^{-1}K$ faces two major hurdles:

1. Computational Intractability: Computing the bounding box of $W(A)$ requires the extreme eigenvalues of the symmetric and skew-symmetric parts of $A$. For $A = \tau M^{-1}K$, these parts are not symmetric in the standard Euclidean inner product, making eigenvalue computation difficult and expensive.
2. Numerical Range Size: $W(A)$ can be extremely large and may extend into the right half-plane (positive real part). Constructing a rational approximation accurate enough to handle large values (e.g., $e^{17}$) within a strict error tolerance is often practically impossible or requires prohibitively high-degree polynomials.

2. Methodology

The authors propose a novel Error Control Framework that bypasses the difficulties of $W(A)$ by utilizing a similarity-transformed matrix.

Key Transformation

Instead of analyzing $A$, the authors define a transformed matrix:
$$\hat{A} = M^{1/2} A M^{-1/2} = \tau M^{-1/2} K M^{-1/2}$$
Since $M$ is SPD, $M^{1/2}$ exists. Note that $A$ and $\hat{A}$ share the same eigenvalues, but their numerical ranges differ significantly in properties.

Theoretical Foundations

The framework relies on three key properties of $W(\hat{A})$:

1. Error Bound via Condition Number: The authors prove (Theorem 1) that the error for the original matrix $A$ can be bounded using the numerical range of $\hat{A}$ and the condition number of $M$:
   $$\|r(A) - e^A\|_2 \leq (1 + \sqrt{2}) \kappa(M)^{1/2} \sup_{z \in W(\hat{A})} |r(z) - e^z|$$
   This allows the construction of $r(z)$ on $W(\hat{A})$ while still guaranteeing accuracy for $A$.

2. Computable Bounding Box: Unlike $W(A)$, the bounding box of $W(\hat{A})$ can be computed efficiently by solving generalized eigenvalue problems.

   • Let $D = (K+K^T)/2$ and $C = (K-K^T)/(2i)$.
   • The real and imaginary bounds of $W(\hat{A})$ correspond to the extreme eigenvalues of the pencils $(\tau D, M)$ and $(\tau C, M)$. These are standard problems solvable with existing numerical linear algebra tools.

3. Left Half-Plane Property: If the numerical range of $K$ lies in the left half-plane (a common property for advection-diffusion operators), then $W(\hat{A})$ is guaranteed to lie entirely in the left half-plane. This prevents the "right half-plane explosion" issue seen in $W(A)$, making the construction of rational approximations numerically stable and efficient.

The Proposed Algorithm (Algorithm 2)

1. Compute the symmetric ($D$) and skew-symmetric ($C$) parts of $K$.
2. Solve generalized eigenvalue problems $(\tau D)x = \mu Mx$ and $(\tau C)x = \mu Mx$ to find the extreme eigenvalues.
3. Construct a rectangular region $R$ enclosing $W(\hat{A})$.
4. Compute the condition number $\kappa(M)$.
5. Construct a rational approximation $r(z) \approx e^z$ on $R$ such that the maximum error on $R$ satisfies:
   $$\max_{z \in R} |r(z) - e^z| \leq \frac{\epsilon}{(1 + \sqrt{2})\kappa(M)^{1/2}}$$
6. Compute the final result as $r(A)b$.

3. Key Contributions

• Theoretical Shift: Introduced the use of the numerical range of the similarity-transformed matrix $\hat{A}$ rather than $A$ for error control in $M^{-1}K$ systems.
• Practical Computability: Demonstrated that the bounding box for $W(\hat{A})$ is computable via generalized eigenvalue problems, whereas $W(A)$ is not easily computable for this matrix structure.
• Stability Guarantee: Proved that if $K$ is stable (left half-plane), $W(\hat{A})$ remains in the left half-plane, avoiding the need to approximate exponentials over large positive real values.
• General Framework: Provided a rigorous error bound that incorporates the condition number of the mass matrix $M$, allowing for prescribed accuracy control.

4. Numerical Results

The authors tested the framework on 2D advection-diffusion equations discretized on square and star-shaped domains using P1 and P2 finite elements.

• Test Cases: Matrix sizes ranged from $n \approx 2,000$ to $n \approx 10,000$. Time steps $\tau$ were set to the mesh size $\bar{h}$ and $10\bar{h}$.
• Methods Compared:
  • Subdiagonal Padé Approximation: Used with a heuristic scaling parameter selection.
  • Rational Interpolation (AAA algorithm): Constructed in partial fraction form.
• Findings:
  • Accuracy: In all test cases, the computed error $\|e^A b - r(A)b\|_2$ remained below the prescribed tolerance $\epsilon$ (ranging from $10^{-8}$ to $10^{-2}$).
  • Efficiency: Using $W(\hat{A})$ significantly reduced the required degree of the rational approximation compared to using $W(A)$. In difficult cases (e.g., $P2$ elements with small diffusion), constructing an approximation on $W(A)$ failed (degree exceeded 320), while $W(\hat{A})$ succeeded with much lower degrees (e.g., degree 25).
  • Large Matrices: The framework scaled effectively to matrices of size $n \approx 10,000$ using sparse solvers and Lanczos methods for eigenvalue estimation.

5. Significance

This work resolves a critical bottleneck in exponential integrators for finite element methods. By enabling reliable error control for matrices of the form $M^{-1}K$, the framework allows for:

1. Adaptive Time-Stepping: Users can dynamically adjust time steps based on guaranteed error bounds.
2. Robustness: It prevents the failure of rational approximations due to the large numerical range of non-symmetric matrices.
3. Scalability: It makes high-order exponential integrators feasible for large-scale PDE simulations where storing full Krylov subspaces is prohibitive, as it relies on rational approximations that can be parallelized and implemented in limited-memory environments.

The paper effectively bridges the gap between theoretical error bounds (based on numerical ranges) and practical implementation for a specific, highly relevant class of matrices arising in scientific computing.