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,通常使用吉莱斯皮算法模拟)的转移概率之间的数学类比。在此过程中,我们提出了一种基于MJP的新方法,该方法依赖于根据矩阵元大小限制“轨迹”跳跃次数,并具有有利的计算缩放特性。随后,我们讨论了该方法对蒙特卡洛后验采样器混合特性的下游影响。我们还对其他两种适用于任意矩阵(不仅限于速率矩阵,更广泛地包括正定矩阵)的矩阵指数计算方法进行了基准测试,这两种方法与求解微分方程相关:龙格-库塔积分法和Krylov子空间法。在最大矩阵元和非零矩阵元数量均与$N$呈线性扩展的条件下——这是常需进行指数运算的速率矩阵的合理条件——最具竞争力方法(Krylov法和其中一种基于MJP的方法)的计算时间缩减为$N^{2}$,总内存需求为$N$。