Nonparametric Modeling of Continuous-Time Markov Chains

Inferring the infinitesimal rates of continuous-time Markov chains (CTMCs) is a central challenge in many scientific domains. This task is hindered by three factors: quadratic growth in the number of rates as the CTMC state space expands, strong dependencies among rates, and incomplete information for many transitions. We introduce a new Bayesian framework that flexibly models the CTMC rates by incorporating covariates through Gaussian processes (GPs). This approach improves inference by integrating new information and contributes to the understanding of the CTMC stochastic behavior by shedding light on potential external drivers. Unlike previous approaches limited to linear covariate effects, our method captures complex non-linear relationships, enabling fuller use of covariate information and more accurate characterization of their influence. To perform efficient inference, we employ a scalable Hamiltonian Monte Carlo (HMC) sampler. We address the prohibitive cost of computing the exact likelihood gradient by integrating the HMC trajectories with a scalable gradient approximation, reducing the computational complexity from $O(K^5)$ to $O(K^2)$, where $K$ is the number of CTMC states. Finally, we demonstrate our method on Bayesian phylogeography inference – a domain where CTMCs are central – showing effectiveness on both synthetic and real datasets.

💡 Research Summary

Continuous‑time Markov chains (CTMCs) are a cornerstone for modeling stochastic dynamics in fields ranging from phylogenetics to epidemiology and chemistry. Estimating the infinitesimal rate matrix Q becomes increasingly difficult as the number of discrete states K grows, because the number of off‑diagonal rates scales as K(K‑1) and many transitions are rarely or never observed. Traditional solutions either impose structural constraints on Q or adopt a log‑linear regression of rates on covariates. The former relies heavily on domain knowledge and can be overly restrictive; the latter assumes a simple linear relationship in the log‑space, which fails to capture complex, nonlinear drivers of transition dynamics.

The paper introduces a novel Bayesian framework that overcomes these limitations by (1) placing a Gaussian process (GP) prior on the relationship between covariates and the log‑rates, and (2) employing a scalable Hamiltonian Monte Carlo (HMC) sampler that uses an approximate gradient of the matrix exponential, reducing computational complexity from O(K⁵) to O(K²).

GP‑based non‑parametric modeling
Each off‑diagonal rate qᵢⱼ is expressed as qᵢⱼ = exp(θᵢⱼ). The latent parameters θᵢⱼ are modeled as a function f of covariates xᵢⱼ: θᵢⱼ = f(xᵢⱼ). Rather than assuming a linear form f(x)=β·x, the authors assign a zero‑mean GP prior to f with kernel k(·,·|α). This allows the model to learn arbitrary smooth functions, automatically adapting to the amount of information available for each state pair. For multiple covariates, an additive GP construction is used, where independent GPs are placed on each covariate dimension, enabling separate control of smoothness and relevance. The GP framework also yields full posterior uncertainty for the inferred rate‑covariate relationship, which is valuable for scientific interpretation.

Efficient HMC with gradient surrogate
The log‑likelihood of partially observed CTMC data involves the finite‑time transition matrix P(t)=exp(tQ) and its derivative with respect to θ. Exact evaluation requires O(K³) for each matrix exponential and O(K²) for each gradient entry, leading to O(K⁵) overall cost per HMC iteration. The authors adopt a recent approximation technique for the matrix exponential gradient (Didier et al., 2024) that replaces the exact derivative with a low‑cost surrogate while preserving the essential geometry needed by HMC. According to Li et al. (2019), using such surrogate trajectories does not compromise the correctness of the Markov chain, provided the approximation error is controlled. Consequently, each HMC iteration scales as O(K²), making inference feasible for state spaces of several thousand states.

Handling partially observed data and identifiability
When observation times are unknown, the authors normalize the rate matrix by the expected number of jumps under a reference stationary distribution, introducing a scaling constant ψ. This yields a normalized matrix Λ = Q/ψ whose time unit corresponds to an expected jump, eliminating the confounding between absolute rates and the time scale. The likelihood is then expressed as an integral over unknown observation times, which can be evaluated analytically or via numerical quadrature within the HMC loop.

Application to Bayesian phylogeography
The framework is applied to phylogeographic inference, where viral lineages migrate among geographic locations. Covariates such as Euclidean distance, population size, and air‑traffic volume are incorporated. The GP prior captures nonlinear effects such as saturation of migration with distance or thresholds due to travel restrictions. Synthetic experiments demonstrate that the method accurately recovers known nonlinear functions, while real‑world analyses of influenza and SARS‑CoV‑2 datasets show improved marginal likelihoods and more plausible covariate effect curves compared with traditional log‑linear models. Posterior predictive checks indicate better fit to observed location patterns, and the inferred rate surfaces provide actionable insights for public‑health interventions.

Broader implications
Beyond phylogeography, the authors argue that the combination of GP‑based rate modeling and O(K²) HMC is applicable to any domain where CTMCs with large state spaces arise, such as chemical reaction networks, queueing systems, and ecological metapopulation models. The paper provides a clear algorithmic description, discusses hyper‑parameter learning for the GP kernels, and supplies open‑source code for reproducibility.

In summary, the paper makes three major contributions: (1) a flexible, non‑parametric Bayesian prior that links covariates to CTCT rates via Gaussian processes, (2) a scalable HMC inference scheme that leverages an approximate matrix‑exponential gradient to achieve quadratic rather than quintic computational scaling, and (3) a demonstration that these innovations substantially improve inference quality in realistic phylogeographic scenarios. The work represents a significant step forward in the statistical methodology for high‑dimensional continuous‑time Markov processes.