Avoiding matrix exponentials for large transition rate matrices

Pedro Pessoa

doi:10.1063/5.0190527

Return

Avoiding matrix exponentials for large transition rate matrices

2024-03-04

Pedro Pessoa

Max Schweiger

Steve Pressé *

DOI：10.1063/5.0190527

Original

Original request for help

Save

Abstract

En 中文

Exact methods for the exponentiation of matrices of dimension N can be computationally expensive in terms of execution time (N3) and memory requirements (N-2), not to mention numerical precision issues. A matrix often exponentiated in the natural sciences is the rate matrix. Here, we explore five methods to exponentiate rate matrices, some of which apply more broadly to other matrix types. Three of the methods leverage a mathematical analogy between computing matrix elements of a matrix exponential process and computing transition probabilities of a dynamical process (technically a Markov jump process, MJP, typically simulated using Gillespie). In doing so, we identify a novel MJP-based method relying on restricting the number of trajectory jumps that incurs improved computational scaling. We then discuss this method's downstream implications on mixing properties of Monte Carlo posterior samplers. We also benchmark two other methods of matrix exponentiation valid for any matrix (beyond rate matrices and, more generally, positive definite matrices) related to solving differential equations: Runge-Kutta integrators and Krylov subspace methods. Under conditions where both the largest matrix element and the number of non-vanishing elements scale linearly with N-reasonable conditions for rate matrices often exponentiated-computational time scaling with the most competitive methods (Krylov and one of the MJP-based methods) reduces to N-2 with total memory requirements of N.

Keywords:

FINITE-STATE PROJECTION

STOCHASTIC SIMULATION