Exact methods for exponentiation of matrices of dimension $N$ can be computationally expensive in terms of execution time ($N^{3}$) and memory requirements ($N^{2}$) not to mention numerical precision issues. A type of matrix often exponentiated in the sciences is the rate matrix. Here we explore five methods to exponentiate rate matrices some of which apply even more broadly to other matrix types. Three of the methods leverage a mathematical analogy between computing matrix elements of a matrix exponential and computing transition probabilities of a dynamical processes (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 based on the magnitude of the matrix elements with favorable 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$.
翻译:精确计算维度为$N$的矩阵指数方法在执行时间($N^{3}$)和内存需求($N^{2}$)方面计算成本高昂,更不必说数值精度问题。科学计算中常需进行指数运算的矩阵类型之一是率矩阵。本文探索了五种率矩阵指数化方法,其中部分方法更广泛适用于其他矩阵类型。其中三种方法利用了矩阵指数元素计算与动力过程(严格来说为马尔可夫跳变过程MJP,通常采用Gillespie算法模拟)转移概率之间的数学类比。据此,我们提出了一种基于MJP的新型方法,该方法通过根据矩阵元素量级限制"轨迹"跳跃次数,实现了有利计算扩展性。随后讨论了该方法对蒙特卡洛后验采样器混合特性的下游影响。我们还基准测试了两种适用于任意矩阵(超越率矩阵及更广义的正定矩阵)的矩阵指数方法——龙格-库塔积分法和Krylov子空间法,这两类方法均与微分方程求解相关。在最大矩阵元素及非零元素数量均与$N$呈线性关系的条件下(这是常需指数运算的率矩阵的合理条件),最具竞争力的方法(Krylov法和基于MJP的其中一种方法)的计算时间复杂度降至$N^2$,总内存需求为$N$。