Exponential time differencing for matrix-valued dynamical systems

TL;DR

This work develops matrix exponential time differencing (METD) methods to solve stiff matrix-valued evolution equations of the form $\dot Q = LQ + QR + N(Q,t)$. It introduces explicit METD1, METD2, and METD2RK schemes for the commuting case $[L,R]=0$, and a general METDp framework with a rigorous local and global error analysis, plus a Baker-Campbell-Hausdorff (BCH) extension for non-commuting operators. A key contribution is the explicit formulation of $I_p$ in terms of backward differences and nested commutators, enabling high-order fixed-step integration without vectorization. Numerical experiments on differential Lyapunov/Sylvester problems, the stiff Allen–Cahn PDE, large-scale jet turbulence, and continuous graph neural networks demonstrate superior stability, accuracy, and efficiency compared to competing fixed-step or adaptive methods. The results show METD's broad applicability to large-scale, high-rank, stiff matrix dynamics in physics, control, and machine learning.

Abstract

Matrix evolution equations occur in many applications, such as dynamical Lyapunov/Sylvester systems or Riccati equations in optimization and stochastic control, machine learning or data assimilation. In many such problems, the dominant stability restriction is imposed by a stiff linear term, making standard explicit integrators impractical. Exponential time differencing (ETD) is known to produce highly stable numerical schemes by treating the linear term in an exact fashion. In particular, for stiff problems, ETD methods are the methods of choice. We extend ETD to matrix-valued evolution equations of the form $\dot Q = LQ + QR + N(Q,t)$ by deriving explicit matrix-ETD (METD) schemes. When $L$ and $R$ commute, we construct an explicit $p$-th order METD$p$ family and prove order-$p$ global convergence under standard assumptions; for the non-commuting case, we develop a Baker-Campbell-Hausdorff (BCH)-based extension. This allows us to produce highly efficient and stable integration schemes. We demonstrate efficiency and applicability on stiff PDE-derived and large-scale matrix dynamics, including an Allen-Cahn system, turbulent jet fluctuation statistics, and continuous graph neural networks. We further show that the scheme is more accurate, stable, and efficient than competing schemes in large-scale high-rank stiff systems.

Figures (3)

• Figure 1: Error at $T=10$ for the Lyapunov equation \ref{['eq:lyap']} measured as $\|C_\infty - C^{(n)}\|$ versus $\Delta t$, where $C_\infty$ solves the stationary Lyapunov equation. The observed slopes confirm the expected global orders for METD$p$ (shown here for $p=1,2,3$).
• Figure 2: Allen--Cahn equation \ref{['eq:allen_cahn_matrix']} with $\varepsilon=0.1$ on a $256\times 256$ grid. (a) The observed slopes confirm the expected convergence orders for METD$p$ with $p=1,2,3,4$ and the Runge--Kutta variant METD2RK. For $\Delta t < 10^{-3}$ the METD4 curve levels off at the numerical error floor. (b) Each METD marker shows the error achieved across each step size; the corresponding RK45 marker shows the runtime required to reach the same accuracy. RK45 runtime is roughly constant because its step size is stability-limited on this stiff system.
• Figure 3: Perturbation on top of zonal jets decaying to a fixed point through the evolution of fluctuations via the differential Lyapunov equation \ref{['dynamicLyapjet']}. The plot on the left shows the stable zonal jet profile $U^*$ with the stable $4$ jet configuration. Integrated using the METD2 scheme.

Theorems & Definitions (12)

• Remark 1
• Remark 2
• Definition 1
• Theorem 1
• Theorem 2
• Theorem 3: Rossmann rossmann2006lie
• Remark 3
• Lemma 1
• Lemma 2
• Lemma 3