BEAGLE in Action: Phylogenetic Likelihoods for Very Large State Spaces

Abstract

Computing phylogenetic likelihoods becomes challenging when the continuous-time Markov chain has a large state space. At each branch, pruning algorithms require terms of the form $e^{Q t}\cdot B$, where $Q$ is an $S \times S$ rate matrix (e.g., for $S$ states), $t$ is a branch length, and the columns of matrix $B$ (e.g., $B \in \mathbb{R}_{+}^{S \times L}$) are partial-likelihood vectors for $L$ site patterns. Many implementations calculate the dense transition matrix $P(t)=e^{Q t}$ via eigendecomposition and multiply it by $B$. This strategy is efficient for small state spaces, but it can be inefficient for large sparse rate matrices because forming $P(t)$ or its eigen vector matrices destroys sparsity and its cost may not be well amortized when $L$ is small. We implement an alternative likelihood kernel based on the Al-Mohy and Higham (2011) algorithm for computing the action of the matrix exponential directly, avoiding explicit construction of $P(t)$ or its eigen vector matrices. The resulting method takes full advantage of the sparsity of $Q$. We implement this method as a CPU-based extension to the BEAGLE library and evaluate it on two large-state molecular models: a Markov-modulated codon model of heterotachy with $M$ selective regimes and $61\cdot M$ possible states, and a codon-pair model of interlocus gene conversion between two paralogs with $61^2=3\,721$ possible states. Our results characterize when action-based likelihood computation improves performance over calculating $e^{Q t}$ directly. We show that this method enables fast, scalable inference for complex models, removing a key computational barrier to the study of evolution with large state space models.