CGSI 2026

Scalable inference forcontinuous-time Markov processes

Applications to phylogenetics

Filippo Monti

UCLA Biostatistics

Evolution Procesess

Interactive demonstration

Single-lineage nucleotide evolution

The live version opens a one-tip nucleotide-only tree view.

Interactive demonstration

Interactive tree app

The live version opens the app-only interactive tree view.

Continuous-Time Markov Chains (CTMCs)

Nucleotide transition graphA, T, C, and G are placed at the corners of a square. Every nucleotide pair is connected by two parallel arrows, one in each direction.ATCGAAdenineTThymineCCytosineGGuanine
Geographic CTMC state graphNorth America, South America, Europe, Africa, Middle East, East Asia, and Australia are placed around a circle. Thin grey arrows connect every pair of states.NAEUMEEAAUAFSANorthAmericaEuropeMiddleEastEastAsiaAustraliaAfricaSouthAmerica

Modeling CTMCs

Questions

  1. How to infer Q?sparsity: not all transitions are represented in the data
  2. Are there external predictors for the entries of Q?

Proposed solution

Gaussian process model:

log⁡(qij)=f(xij)+ϵ,f(⋅)∼GP(0,K(⋅,⋅)).\log(q_{ij}) = f(x_{ij}) + \epsilon, \qquad f(\cdot) \sim \mathrm{GP}(\mathbf{0}, K(\cdot, \cdot)).

Global Flu example

Virus:
Global Influenza H3N2
Host:
Humans
When:
2002-2007
Data:

1441 H3N2 hemagglutinin sequences + locations

Log rates as a function of origin population density
Global influenza phylogeny with geographic reward summaries

Monti, F., Ji, X., Shao, Y., and Suchard, M. A. (2026). Nonparametric Modeling of Continuous-Time Markov Chains with Scalable Exact and Approximate Gradients.

Computational contribution

Let NN be the number of observations and KK the number of states.

Challenges

  • QQ has K2K^2 highly dependent entries gradient-based inference.

  • QQ enters the likelihood through etQ:=∑n=0∞(tQ)nn!e^{tQ} := \sum_{n=0}^{\infty} \frac{(tQ)^n}{n!}.

  • The gradient of etQe^{tQ} is super expensive:

    O(O(NK3)N K^3)

Solution

  • New adjoint-based matrix-exponential gradient to get:

    O(O(NK2)N K^2)
Gradient benchmark comparing exact and approximate methods against autodiff and PhyloGrad

Monti, F., Ji, X., Shao, Y., and Suchard, M. A. (2026). Nonparametric Modeling of Continuous-Time Markov Chains with Scalable Exact and Approximate Gradients.

Trait evolution as continuous-time diffusion

dX(t)=SDriftpulls morphology toward equilibrium{X(t) -θEquilibrium mean}dt+Σ1/2Diffusioninjects stochastic variation along the branchdW(t)
start of branch
end of branch
t = 0.00

Evolving traits: body size, head, forelimb, hindlimb, lamellae, and tail.

Inference of the drift matrix S

Challenges

(1) S must be (Hurwitz) stable: (Re⁡(λ)≤0)(\operatorname{Re}(\lambda) \le 0)

(2) computational complexity of branch variances(lots of etSe^{tS});

∫0te(t−τ)SΣe(t−τ)STdτ\int_{0}^{t} e^{(t-\tau) S} \Sigma e^{(t-\tau) S^{T}} d\tau

Current scalable solutions: diagonal or symmetric parametrizations

But: they can only describe “monotone” evolution

Using Schur blocks:

(αβ−βα)\begin{pmatrix} \alpha & \beta \\ -\beta & \alpha \end{pmatrix}

But: need to pre-specify the number of complex conjugate pairs

Proposed parametrization

A matrix S∈RK×KS \in \mathbb{R}^{K \times K} is represented under the Hurwitz-smooth spectral block parametrization (H-SSBP) if it can be written as

S=RBR−1,S = RBR^{-1},

where B=bdiag⁡(B1,…,BnB)B = \operatorname{bdiag}(B_1,\ldots,B_{n_B}) is block diagonal, with

Bb=B(ρb,θb,ωb):=(ρbcos⁡θbρbsin⁡θb+ωbρbsin⁡θb−ωbρbcos⁡θb)B_b = B(\rho_b,\theta_b,\omega_b) := \begin{pmatrix} \rho_b \cos \theta_b & \rho_b \sin \theta_b + \omega_b \\ \rho_b \sin \theta_b - \omega_b & \rho_b \cos \theta_b \end{pmatrix}

If KK is odd, one block is scalar.
Smooth transition between coupled complex eigenvalues and decoupled real eigenvalues
Stability is guaranteed by ρb<0,θb∈[−π4,π4]\rho_b \lt 0, \theta_b \in [-\frac{\pi}{4}, \frac{\pi}{4}]
Matrix support: all matrices diagonalizable over C\mathbb{C} and some defective cases   ⟹  \implies full Lebesgue measure over the stable cone

Monti, F., Holbrook, A., Glatt-Holtz, N. E., and Suchard, M. A. (2026+). Stable Matrix Parametrizations and Structured Adjoints for Ornstein—Uhlenbeck Processes.

Computational contribution

State of the art:Gaussian message passing algorithm

O(O(NK3)N K^3)

Using H-SSBP:

O(O(NK3)N K^3)
Benchmark figure for matrix exponential, matrix exponential adjoint, Lyapunov equation, full edge message, and full edge message gradient computations

Monti, F., Holbrook, A., Glatt-Holtz, N. E., and Suchard, M. A. (2026+). Stable Matrix Parametrizations and Structured Adjoints for Ornstein—Uhlenbeck Processes.

Thank you

Advisor: Marc Suchard

Collaborators: Xiang Ji, Yucai Shao, Andrew Holbrook, and Nathan Edward Glatt-Holtz.

Funding:

National Institutes of Health grants U19 AI135995, R01 AI153044, and R01 AI162611.