Variational Phylogenetic Inference with Products over Bipartitions
Authors: Evan Sidrow
Alexandre Bouchard-Cote
Lloyd T Elliott
Venue: ICML 2025 poster (2025)
Keywords: Phylogenetic Inference
Variational Bayes
COVID-19 Genetics
Linkage Clustering
Reinforce Estimators
Primary area: probabilistic_methods->variational_inference
TL;DR: We develop a variational inference approach for ultrametric phylogenetic trees that is differentiable, doesn't restrict tree space, and doesn't rely on MCMC subroutines.

Abstract: Bayesian phylogenetics is vital for understanding evolutionary dynamics, and requires accurate and efficient approximation of posterior distributions over trees. In this work, we develop a variational Bayesian approach for ultrametric phylogenetic trees. We present a novel variational family based on coalescent times of a single-linkage clustering and derive a closed-form density for the resulting distribution over trees. Unlike existing methods for ultrametric trees, our method performs inference over all of tree space, it does not require any Markov chain Monte Carlo subroutines, and our variational family is differentiable. Through experiments on benchmark genomic datasets and an application to the viral RNA of SARS-CoV-2, we demonstrate that our method achieves competitive accuracy while requiring significantly fewer gradient evaluations than existing state-of-the-art techniques.

Variational Phylogenetic Inference with Products over Bipartitions
Evan Sidrow Alexandre Bouchard-Côté Lloyd T. Elliott
Abstract
One alternative to MCMC is variational inference (VI), in which the posterior distribution over phylogenetic tree struc- tures is approximated using a variational distribution that minimizes some distance metric to the true posterior distri- bution. While VI over combinatorial spaces is also known to be difficult due to the complexity of the constraints on the support (Bouchard- Côtéé and Jordan, 2010; Linderman et al., 2018), there have been several recent advances in VI over phylogenetic trees that are tractable. Zhang and Mat- sen IV (2018) 8)) represented phylogenetic trees as Bayesian networks using subsplit Bayesian networks (SBNs), and later used SBNs to perform variational Bayesian phyloge- netic inference on unrooted trees (Zhang and Matsen IV, 2019). This approach has spawned many methodological ad- vancements. To improve the distribution for branch lengths, Zhang (2020) used normalizing flows, Molén et al. (2024) used mixtures, and Xie et al. (2024) used semi-implicit branch length distributions. To improve the variational fam- ily over tree topologies, Zhang (2023) used graph neural networks to learn topological features.
Bayesian phylogenetics is vital for understanding evolutionary dynamics, and requires accurate and efficient approximation of posterior distributions over trees. In this work, we develop variational Bayesian approach for ultrametric phylogenetic trees. We present a novel variational family based on coalescent times of a single-linkage clustering and derive a closed-form density for the resulting distribution over trees. Unlike existing methods for ultrametric trees, method performs infer- ence over all of tree space, it does not require any Markov chain Monte Carlo subroutines, and our variational family is differentiable. Through experiments on benchmark genomic datasets and an application to the viral RNA of SARS-CoV-2 we demonstrate that our method achieves compet- itive accuracy while requiring significantly fewer gradient evaluations than existing state-of-the-an techniques.
The number of parameters within an SBN grows exponen- tially with the number of taxa, so Zhang and Matsen IV (2019) used MCMC to find sets of most likely tree struc- tures, and they use these sets to restrict the SBNs. Other recent VI approaches sample tree topologies without the use of SBNs. For example, ViaPhy (Koptagel et al., 2022) uses a gradient-fre variational inference approach and directly sample from the Jukes and Cantor (1969) model, GeoPhy (Mimori and Hamada, 2023) uses distance-based met- ric in hyperbolic space construct unrooted phylogenetic trees, and ARTree (Xie and Zhang, 2023) uses graph neu- ral networks to construct deep autogressive model for VI over phylogenetic tree structures. Zhou et al. (2024) also introduce PhyloGFN, a phylogenetic VI technique based on reinforcement learning and generative flow networks (Bengio et al., 2023).
1. Introduction
The goal of Bayesian phylogenetics is to infer the genealogy of a collection of taxa given genetic model and aligned sequence data. Phylogeneties is used in fields such as epi- demiology (Li et al., 2020), linguistics (List al., 2014), and ecology (Godoy et al., 2018). 8). Bayesian phylogenetic infer- ence quantifies uncertainty and integrates over phylogenetic tree structures within phylogenetic model (Zhang and Mat- sen IV, 2019). Most Bayesian phylogenetic inference is per- formed using Markov chain Monte Carlo (MCMC) methods with candidate trees iteratively proposed and either accepted or rejected based on their consistency with the observed data. However, MCMC methods can struggle because the number of possible trees grows eypereexponentially in the number of taxa and posteriors on trees are highly multi-modal,
Phylogenetic inference also be performed using non- Bayesian methods, including RAXML (Stamatakis, 2014), neighbour-joining (Saitou and Nei, 1987), and more recently Phyloformer (Nesterenko et al., 2025). Phyloformer uses deep learning to construct pairwise representations of evo- lutionary distances between taxa. Phyloformer then uses pairwise distances construct a tree using a neighbour- joining algorithm similar to the single-linkage clustering
'Department of Statistics and Actuarial Science, Simon Fraser University, Burnaby, Canada *Department of Statistics, University of British Columbia, Vancouver, Canada. Correspondence to: Evan Sidrow Scsidrow@sfu..ca.
Proceedings of the nnnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Variational Inference with Products over Bipartitions
algorithm described here. However, non-Bayesian methods do not provide estimates of marginal likelihood, which are useful for model selection.
the difference in coalescent times between the parent and the child nodes. We measure t in backwards time, so each i, is positive and represents a time before the present. The leaves of correspond to the genomes of each measured taxon x E X. Additionally, an internal node u of 7 represents the (unobserved) genome of the most recent common ancestor of all taxa that have u as a parent node, As a binary tree, r contains a total of N internal nodes (including the root node). We can thus represent the tree with collection of bipartitions: = {{W,. Zn, J In this representation, an internal node uy, is the most recent common ancestor for all taxa x E Wn, U Zn. and the n-th coalescent event is represented by the bipartition {W,, Z.
In this work we focus on rooted ultrametric phylogenetic trees, for which branch lengths correspond to the amount of time between evolutionary branching events. This for- mulation is useful when time important, for example in applications involving rapidly evolving pathogens (Sagu- lenko al., 20. None of the aforementioned approaches incorporate time constraints into the branch lengths of the phylogenetic trees. To this end, Zhang and Matsen IV (2024) generalized their SBN-based approach to ultrametric trees, but they still rely on MCMC to restrict to a subset of tree space upon which to perform inference over. Alternatively, Bouckaert (2024) provides cube a method related to the one described in our manuscript. However, because the matrix representation of tree space in Bouckaert (2024) is not dense, they can only express a limited number of trees. As such, their variational family is not supported on many tree topologies. Further, they do not perform optimization on the tree structure, and instead rely on empirical values derived from MCMC.
Denote the set of possible characters within a set of aligned genetic sequences by Q (e.g., a DNA sequence may cor- respond to Ω {A, T, G. and an RNA sequence to Ω {A, U, G, cj). Further, denote the set of ob- served genomes by Y(o)) = {Y ..... ...,,,)(o)) y/ob)) y=(ah) = (Ym.r....... Ym,zx) corresponds to the base pairs at site m for all observed taxa æ E X. In addition to the observed genomes ( denote the unodserved genomes of all internal nodes by ( (r,('n) y-(um) = (Ym.u....., Ym,.=,)) corresponds the base at site m for all unobserved internal nodes u1,.... -I- Let the index of the root node be N - 11 (so, uN-1 is the root node). We denote the combined observed and unobserved genomes by Y {Y Y }. For further background phylogenetics, we refer to Hein al. (2004).
where
y-("m) J. where
We introduce Variational phylogenetic Inference with PRoducts over bipartitions (VIPR), new variational family for ultrametric trees based on coalescent theory and single- linkage clustering (Kingman, 1982). VIPR naturally per- forms variational inference on ultrametric trees and thus di- rectly incorporates time into phylogenetic inference. VIPR also performs inference over the entirety of tree space and does not rely on MCMC subroutines. In particular, we pa- rameterize a variational distribution over a distance matrix and use it to derive differentiable variational density over trees that result from single-linkage clustering. Through a set of experiments on standard datasets and an application to COVID-19, we show that our simple variational formu- lation achieves comparable results to existing methods for ultrametric trees in fewer gradient evaluations.
2.1. Phylogenetic Likelihood
For simplicity, we focus on the Jukes and Cantor (1969) model of evolution. We denote the stationary distribution by т, and the transition matrix by P(b) (a 4 x 4 matrix such that the i, j-h entry is the probability of transitioning from base i to base j given branch length b under Jukes and Cantor 1969). With this notation, the likelihood of an observed set of genetic sequences (ob) at site mn is as follows:
2. Background
Consider a set of N taxa X = {01, F2....,..... JN}. A nonempty subset X of X' is referred to as a clade of X. clade represents a collection of taxa which share a com- mon ancestor at a particular time in the past. Further, we represent an evolutionary branching events using biparti- tion (W. Z) of the clade X (x= = W U Z, W n Ø)).
p( (( m
|T, [T,t)=
Σ TT (en) (u.0)
II (P(b..o(T,) Yu.."
Here the product is over all edges (u, v) in T. We assume independence between sites, and so the likelihood of the observed genomes is:
We focus on ultrametric trees, in which the leaves of the trees are all equidistant from the root. We denote an ultra- metric tree with a rooted, binary tree topology T and a set of coalescent times {t., W where there one t, for each internal node in r. For ultrametric trees in particular, the branch length between child and its parent is equal to
p (( [T,t)= = II
p( (o T, (T,t).
(1)
Equation (1) can be evaluated in O(NM) time using the pruning algorithm (Felsenstein 1981, also known as the sum-product algorithm; Koller and Friedman 2009).
2

---

Variational Inference with Products over Bipartitions
2.2. Prior Distribution over Trees
Algorithm Single-Linkage Clustering(1] Ro)
We use the Kingman coalescent (Kingman, 1982) the prior distribution on trees. This coalescent process pro- ceeds backward in time with independent and exponen- tially distributed inter-event intervals. Events occur at rate A = (5) /Ne. Here k is the number of extant taxa and No is the effective population size, a parameter which governs the rate at which taxa coalesce. We fix N₂ = 5 in our ex- periments. At each event, a pair of extant taxa are chosen to coalesce into single taxon uniformly at random over all pairs of extant taxa, yielding the prior:
1: Input: Distances T E R200 and taxa set No
((11), (12). {xx)).
2: for I......N-- - do
3: w²,=" + arg min (w(w))) ww, 2 not coalesced)
4: Set Wn, E No to be the set containing w'
5: Set Zn E 2o to be the set containing 2*
In t Hw 3
Remove Wn, Zn from Ro and add W,, U Zn to No
8: end for
2 p(T,t) NI(N - 1)! k=2
N II
9: (WW,, Zn 1
A exp (~>o(th - tx-1))
10: t + {In J
11: Return (T,t)
2.3. Variational Inference for Phylogenetic Trees
Our goal is to infer distribution over tree structures and coalescent times given observed genetic sequences:
2.4. Matrix Representation of Tree Space
One way to construct a phylogenetic tree is to use distance matrix T (a symmetric N x NN matrix with positive and finite off-diagonal entries) and the single-linkage clustering algorithm, as described in Algorithm 1. We denote the distance between taxa u and v by ((w,2) and formulate the phylogenetic tree
p(T, t Y((b)) = p(Y((b))
t)/p()) (2)
Here p( (a)) is an intractable normalizing constant. Vari- ational inference involves defining a tractable family of probability densities parameterized by some variational pa- rameters é. Then, the posterior density is approximated by a variational density (T, t) whose parameters & should minimize a divergence measure D between the posterior p(- (ob)) and 4s- Here we use the reverse KL divergence:
algorithm return a representation using bipartitions and coalescent times that is consistent with our notation. We consider the distance matrix as an
(?)
element of R 50 by identifying the off-diagonal elements.
Algorithm is a naive implementation of single-linkage clustering with time complexity (NN) and space complex- ity O(N2). Sibson (1973) introduce SLINK, an implemen- tation of Algorithm 1 with time complexity ((V2) and space complexity ((N), and prove that both are optimal.
2o (T, (T) t.) p(T, t ( [ )]
DKI(gs p) = B( log
(3)
Bouckaert (2024) introduce cube a method that uses clustering to perform variational inference over ultrametric trees. First, they note that if exactly N entries of T are finite, then single-linkage clustering im- plies bijection between T and (r, tr. Then, they specify exactly N - entries of T to be random and finite, setting all other entries to infinity. Next, they run an MCMC al- gorithm over trees to obtain an empirical distribution over coalescent times. Finally, they estimate the parameters asso- ciated with the N finite entries of using the MCMC- generated empirical distribution. This method has two main drawbacks-all entries of T must be specified to form a distribution supported on the entire tree space, and it is un- clear how to select the N - best entries of T to cover the most posterior probability. To address these drawbacks we present a gradient- based variational inference method for ultrametric trees based on single-linkagee clustering that specifies variational distribution over the entire tree space.
Evaluating the exact posterior p(r, t Y (ab) is difficult. Instead, we equivalently (up to normalizing constant) max- imize the evidence lower bound (ELBO):
Φ² = arg max L(Ø).
(4)
p(r,t, ((())) )(
L(ø)= = Egg
log
(5)
(o t)
ELBO is also known as the negative variational free energy in statistical physics and some areas of machine learning. The expectation over do consists of sum over tree struc- tures 7 and an integral over coalescent times t, forming the following objective function:
I(O) = Σ (T, (o(T,,))log tog (p(T, t, Y(ob))
dt.
(6)
(T) t)
3

---

Proposition If the random variables (la.) are mutually Variational Inference with Products over Bipartitions
83 = 4
W {A, B}
Z3 {C, D)
12 =
W2 (C)
2
8
4
2
5
6
h =
Zz = (D)
B
# *
4.5
7
*
*
4
7
Wi {A}
T
T
C
$
#
*
3
#
8
8
3
Z1 (B}
D
#
#
*
A
B
B C D
Figure 1. Schematic showing the sampling process for VIPR. This diagram shows two possible example matrices (on the left) that could be drawn using gfu.,) ~ Lo (u,")) and result in the same phylogenetic tree (T, t) ~ qo (on the right) after running single-linkaee
(2,0) is quite general and can be provided by clustering. Entries of T that trigger coalescence event shown in bold. The form of 93
the practitioner, while expression for qs depends upon 4 (u,)).
3. Methods
(r, t) formed by rumning single-linkage clustering on T. Al- gorithm 2 presents pseudocode to sample from de if d is a log-now log-normal distribution, while Figure 1 visualizes the process of drawing T using [(u,,} Lo single-linkage clustering to map T to (T,
[u,v]
We now present our method Variational phylogenetic Inference with PRoducts over bipartitions (VIPR). We be- gin by outlining a generative process for sampling from variational distribution qe. We then describe how to evaluate the density of our variational distribution. Finally, we describe the optimization procedure used to maximize our variational objective function.
(u.r} and then using
3.2. Density Evaluation
In this section we describe how to evaluate the density of trees generated from Algorithm 2. The primary challenge is that a given tree (r, t) may have been generated from multiple distance matrices (see Figure 1). Luckily, this sampling procedure still yields a density with a closed-form solution, as shown in Proposition below.
3.1. Generative Process for Phylogenies
We begin by describing a generative process for sampling from 4s. as our variational distribution is best understood through the algorithm for sampling from it. Our algorithm to sample an ultrametric tree with leaf nodes X' proceeds similarly Algorithm 1 of Bouckaert (2024). Namely, randomly draw each element of the distance matrix T (tfu,e) for all {u, a with v E X)) using a set of indepen- dent variational distributions with densities y (u.e). Then, we run single linkage clustering on to form (T, t).
independent, and all d (u.w) are continuous in & and ffor all {u, u) with v E X, and Q2 (u..) is the survival function of (luo), then qo (T, 1) has the following form:
Algorithm 2 Sample-q(u, xx))
N-1 II
g
fv, (bn) [.
(T
1: Input: Parameters μa R()) and o R and taxa
(tm)
W
WEW. 7
2E Zn
set X [[11), (x.)...,
2: Draw (0.1) ~ NN0. 1) for all {u, w) C X (N)
Define matrix T e R such that:
A derivation of Proposition 1 using induction is provided in Appendix A. Every taxa pair {u, d} appears in the sum and product terms of Equation (7) exactly once, as each taxa pair coalesces exactly once within a rooted phylogenetic tree. Thus, evaluating both go(r, t) and V log (T t) takes ( time.
log (f(M.)) / (m.0) + 2 (1)) (v.))
4:
Return Single-Linkage Clustering
4/ Note that (u.v) (1tu.e) and (T, t) are closely related:
[u,r] describes the distribution over entry f(un)) T, do
If % [u,²]
is continuously differentiable, then do also continuously differentiable. our VIPR implementation,
while go describes the distribution over phylogenetic trees

---

Variational Inference with Products over Bipartitions
[1.t) is log-normal, so we can compute gradients with re- Ls
INFORCE (LOOR) estimator (Mnih and Gregor, 2014; Shi et al., 2022). The gradient V log (o(r((), (k), (()) can be cal- culated using automatic differentiation software such as Autograd (Maclaurin et al., 2015) or PyTorch (Paszke et al., 2019).
spect ø. Note however that Proposition I holds for any continuous mutually independent (u.))
Algorithm 3 VIPR(X, K)
1: Input: Integer K indicating number of samples to
3.3.2. THE REPARAMETERIIATTIN TRICK
use in gradient approximation and taxa set X
{221),22)......
The push out estimator (Rubinstein, 1992) popular in machine learning literature under the name of the reparam- eterization trick (Kingma and Welling, 2014). Recall that in our experiments q ('".)) is log-normal distribution for all u, v X. In Algorithm (2), the candidate coalescent times (la.o) log-normal(ufuw0) ( (u,w))) < ((n.)
2: Initialize variational parameters Φ
3: while not converged do
4: for k 1k=1,... 1,..... do
fun) 5: Draw T(A) E R500 with f(uo) ~ oz
exp(u(u.) + o (1 with 2( (0. De-
6: (r(h, (())), +
noting the set {2(u,²) JumEx by Z, we reparameterize the expectation in Equation (5) as follows:
Linkage Clustering (T()),
7: end for
8: Estimate gradient V (()) using (r(*), t(()))
P 96
(2)) (2)
((Ø) = Ez log
(8)
for k =. K.
9o (0 (Z))
9: Update e using gradient estimates and a stochastic optimization algorithm (Adam, SGD, etc.)
Here gøl (Z)
Clustering + O
2). Denoting the density of (*) -dimensional standard
10: end while
normal distribution as O. O, we have:
11:
Return
(( (o (z))
3.3. Gradient Estimators for do
[( = OZ log
dz. (9)
We now have almost everything we need to perform phy- logenetic variational inference: an (unnormalized) phylo- genetic posterior density p(T, t, Y( ). a variational fam- ily with density qT(T, t). and an objective function L()) to maximize in order find variational posterior distri- bution. We collect these steps together in Algorithm 3. Note that in this algorithm, optimizing L(Ø) with stochas- tic gradient methods such as Adam (Robbins and Monro, 1951; Kingma and Ba, 2014) requires random estimates of the gradient L(). Thus, we consider three meth- ods for gradient estimation and use them in our exper- iments: leave-one-out REINFORCE (Mnih and Gregor, 2014; Shi et al., 2022). the reparameterization trick (Ru- binstein, 1992; Kingma and Welling, 2014), and VIMCO (Mnih and Rezende, 2016). An overview of these methods are given in the remainder of this subsection, with details in Appendix C. Our code implementing VIPR is available at https://github. htt.a/ab..oom/vvaaaarrowo//////e /VIP..
Using this formulation and some additional algebra, in Ap- pendix C.2 we derive an estimator for the full gradient V (() using random samples of Z. Unfortunately. this estimator is biased because the gradient and integral can- not be interchanged without introducing some error (see Appendix C.2). Therefore, this optimization procedure is not guaranteed to converge to a local optimum of the objec- tive function. Nonetheless, these gradient estimates tend to perform at least comparably to the LOOR estimator.
3.3.3. THE VIMCO ESTIMATOR
One drawback of the single-sample ELBO in Equation (5) is that variational distributions that target the ELBO tend to be mode-seeking (i.e., they can underestimate variance of the true posterior). As an alternative, Mnih and Rezende (2016) suggest a K K-sample ELBO (VIMCO: variational inference for Monte Carlo objectives) that encourages mode-covering behaviour in the posterior. We derive VIMCO Estimator for model in Appendix C.3.
3. THE REINFORCE ESTIMATOR
Define (T, log(p(r, t, Y((h))) log(gs(T, (T,t)., that [(Ø)= E.y., [T. t)]. We can interchange the gradient and the finite sum over 7 in Equation (6), and we assume that we can interchange the gradient and integral (see L' E'Ecuyer 1995 for technical conditions). After performing some al- gebra (see Appendix C. we obtain the leave-one-out RE-
4. Experiments
We compared the performance of our VIPR methods with that of Zhang and Matsen IV 2024 (denoted VBPI this section). We do not compare VIPR to cube VB (Bouck-
5

---

Variational Inference with Products over Bipartitions
aert, 2024) because Cu does not involve maximization and therefore is not directly comparable. However, we do investigate how many tree topologies in the posterior fall outside of the restricted cube space of CB in Appendix B.1. We studied eleven commonly used genetic datasets that are listed in Lakner et al. (2008) denoted DS1 through D (these are the names that are given to these datasets in Lakner et al. 2008).
likelihood values using the stepping-stone estimator (Xie et al., 2010). For each dataset, we ran 100 path steps of 500,000 MCMC iterations and repeated this process ten times to obtain 10 independent estimates of the MLL.
4.2. The VBPI Baseline
We compared VIPR to the VBPI algorithm as implemented by Zhang and Matsen IV (2024), which requires MCMC runs to determine likely subsplits (i.e evolutionary branch- ing events). To provide these runs, we used BEAST to obtain a rooted subsplit support. We ran 10 indepen- dent MCMC chains for 1,000,000 iterations, with the first 250,000 discarded as burn-in. We then thinned to every 1,000-th iteration, yielding 7,500 trees for the VBPI subsplit support.
We also studied dataset of 72 COVID-19 genomes ob- tained from GISAID (Global Initiative on Sharing All In- fluenza Data; Khare et al. 2021). In particular, we obtained COVID-19 19 RNA sequences that were collected in Canada on January 2, 2025; submitted to GISAID prior to Jamuary 20, 2025; contained at least 29,00 29,000 sequenced base pairs; and were of the strain JN The 72 COVID-19 genomes studied here are all of the COVID-19 genomes provided by GISAID that satisfied all of these criteria. After obtaining these genomes, we aligned them using multiple sequence alignment in MAFFT using the FFT-NS-1 algorithm (Katoh and Standley, 2013). Finally, we subset the genomes to M = 3,101 hon-looolog sites (i.e., we omitted all sites that were the same across all 72 taxa). The final datasets ranged from 27 to 72 total taxa N and 378 to 3,101 total sites M. See Appendix B.2 for the number of taxa and sites by dataset as well as other summary information.
To fit the VBPI baseline, we used the VIMCO gradient esti- mator with K K-sample ELBO for K = 10 and K = 20 (indi- cated by VBPI1 and VBPI20 our plots and tables below). Zhang and Matsen IV (2024) use an annealing schedule dur- ing optimization, but we omitted the annealing schedule to be consistent with our optimization for VIPR. We used the Adam optimization algorithm implemented in PyTorch with four random restarts and learning rates of 0.003, 0.001, 0.0003, and 0.0001 (Kingma and Ba, 2014; Paszke et al., 2019). We estimated the marginal log-likelihood (MIII) every 100 iterations (i.e., parameter updates) using 500 importance-samplee particles.
For all methods considered (BEAST, VBPI and our VIPR methods), we used a Kingman coalescent prior on the phy- logenies. We fixed the effective population size at N₂- (Kingman, 1982) and assumed the Jukes- Lukes-Cantor model for mutation (Jukes and Cantor, 1969). These assumptions are described above in Sections 2.2 and 2.3. We also measure the branch lengths in terms of expected mutations per site, which is in line with BEAST and VBPI.
Of the 16 runs for each VBPI batch size condition (4 learn- ing rates and 4 random restarts), we retained the run with the highest average MLL in the last 10 estimates of the run. This run (with highest average MLL) was included in our plots and figures. We used the primary subsplit pair (PSP) parameterization of VBPI. Code for these experi- ments was adapted from https:// http:/gttthbhhooooee com/zcrabbit/ vbpi- /teaee/eee//eeeeeeeeeeeee weaecc//tee/mminniooieee main/rooted. See Zhang and Matsen IV (2024), Section 6 for more implementation details.
Each run for the experiments on the DS1 to DS11 datasets was executed on supercomputer node. The runs were al- located 12 hours of wallclock time, 1 CPU, and 16GB of RAM. The supercomputer had heterogeneous infrastruc- ture involving in which each CPU make and model was Intel v4 Broadwell, Intel Caskade Lake or Skylake, or AMD EPYC 7302. Experiments on the COVID-19 dataset were run with identical conditions to those for DS1 to DS11, but without the 12 hour limit on wallelock time. Instead, they were run for 10,000 iterations (i.e. parameter updates) or 12 hours (whichever took longer).
4.3. The VIPR Methods
For our VIPR methods, we set the variational distribu- tions g fu,") to be log-normal, so the variational parameters & were the means and standard deviations corresponding to the logarithm of the entries log(t((.r) of the matrix T. After running BEAST, we plotted histograms of pair- wise coalescent times across sampled trees for all datasets. In most cases these histograms looked approximately log-normal, motivating our choice of log-normal distri- butions (see https://githut com/ heem//vaaSiidw///ie /Viro////////////// tree/ amainnssupmat hists for the histograms). We also simulated data and ran posterior predictive checks in Ap- pendix B.3 for model checking.
4.1. The BEAST Gold Standard
To approximate the true posterior distribution of each dataset we ran 10 independent MCMC chains using BEAST, each with 10,000,000 iterations. We discarded the first 250,000 iterations as burn-in and thinned to every 1,000-th iteration.
This yielded in a total of 97,500 trees that were used as a "gold standard." " We estimated ground- truth marginal log-
To initialize the parameters o, we computed the empirical
6

---

Variational Inference with Products over Bipartitions
distribution of coalescent times between taxa (u.v) from the short MCMC runs used to establish the support for the VBPI baseline. We then set the initial mean and standard deviation of (u,0) to be the mean and standard deviation of the empirical distribution.
half of the datasets). This is likely because VBPI targets multi-sample ELBO for optimization, which produces mode-covering behaviour. In contrast, our VIPR methods target the single-sample ELBO.
We experimented with three gradient estimation techniques. We estimated (o)( using (1) the LOOR estimator and (2) the reparameterization trick, both with batch sizes of 10 samples. We also estimated V (Ø)) using the VIMCO estimator with batch size of K = 10.
Table 1. Gap between gold standard and estimated marginal log- likelihoods for variational inference methods (in nats). Marginal log-likelihoods for VI methods were estimated using importance sampling with 1,000 random samples from each variational dis- tribution. Values indicate difference between gold standard MLLs and each method's s MLLs. Gold standard MLLs (indicated by the BEAST column) are derived from 10 independent chains of the stepping-stone algorithu in BEAST. Datasets (DATA column) DSI to DS11 are from Lakner el al. (2008). Dataset COV is the COVID- 19 dataset obtained from GISAID. VI methods are specified by columns: Variational Bayesian Phylogenetic Inference with K- sample BLBO, K = 10 (VBPT10; Zhang and Matsen IV 2024); VIPR using the leave-one out REINFORCE estimator (LOOR); VIPR using the reparameterization trick (REP); VIPR using the Variational Inference for Monte Carlo Objectives estimator with 10 (VIMCO). Results for VBP120 and standard errors are in Appendix B.4.
For each gradient estimation technique, we used the Adam optimizer in PyTorch with ten random restarts and learning rates of 0.001, 0.003, 0.01, and 0.03. We recorded the estimated MLL every 10 iterations (i.e., parameter updates) with 50 Monte Carlo samples. Of the 40 runs (4 learning rates and 10 random restarts), we retained the run with the highest average MILL in the last 10 estimates of the run. As for VBPI, the retained run (for each dataset and technique) is reported in the plots and figures below.
4.4. Simulated Datasets
To explore the runtime and asymptotic complexity of the methods as a function of the number of taxa, we created a
DATA
BEAST
VBPIIO
LOOR
REP
VIMCO
DS1
-7154.26
-0.53
2.29
-1.83
-0995
series of simulated datasets. Through this simulation, we
DS2
-26566.42
0.16
-0.76
-0.14
0..7
could keep the number of sites and the underlying tree dis-
tribution and substitution model all constant. We simulated seven datasets using the ms software (Hudson, 2002) with 1,000 sites, an infinite sites model (Kimura, 1969), a neutral model of evolution, and the Kingman coalescent. These assumptions imply the software command 'ma (N) -T -s 1000.' Here (N) indicates the number of taxa, which we varied in the set {8, 16, 32, 64, 128, 256, 512). We refer to these datasets as the MS datasets.
DS3
-33787.62
-0.44
-3.66
-1.91
-266
DS4
-33506.05
0.03
-2.48
-0.47
- -1.73
DS5
8271.26
-1.70
-0.29
-4.01
0.94
DS6
-6745.31
-0.6
33.96
33.26
-272
DS7
37323.88
0.27
2.73
2.82
10.42
DS8
8650.20
0.82
33.8
4.95
-2.88
DS9
4072.66
-5.32
3.12
-5.79
- -7.60
DS10
10102.65
0.88
5.38
33.8
6.82
5. Results
DS11
-6272.57
-18 -1.79
-6.79
-7.31
-9.62
COV
7881.61
-39.1
-611
-374
2214
Tables 1 and 2 show the estimated MLLs and ELBOs for our variational inference experiments after 12 hours of com- pute time. Results for VBPI with in Appendix B.4. The stepping-stone stone algorithm is not variational method, and so it has entry in Table 2.
Figure 2 shows the trace of estimated log-likelihne likelihood ver- sus iteration number for all VI methods on DS1 and on the COVID-19 19 dataset. See Appendix B.4 for results on DS2-11. 11. These figures also display empirical distributions of tree met- rics for each VI method's learned variational distribution in addition to the BEAST gold standard run (plotted with mat- plotlib's kde function with default parameters; Hunter 2007). We removed of the 1,000 trees sampled from VBP120 for the COVID-19 19 experiment because they had extremely low log- lngllikelssoods (< - 90,000), resulting in flat densities.
20 (VBPI20) are
The MLLs in Table are reported by the gap between the MLL of the gold standard (BEAST/stepping stone run) and the method's 's MLL (the difference between the MLLs). Methods with smaller gaps are therefore closer to the gold standard, and the method with the highest MIL is bolded. Note that some VI methods surpass the gold standard, likely due to Monte Carlo error.
VIPR tended to underestimate the variance of tree length compared to BEAST, while VBPI tended to overestimate. In addition, the reinforce and reparameterization gradient estimates result in variational distributions with higher tree log-likelihoods on average, while VBPI and our VIPR with
VBPI tends to slightly outperform our VIPR methods in terms of MLL, but all methods are comparable in terms ELBO (with our methods outperforming VBPI on exactly
7

---

Variational Inference with Products over Bipartitions
(a)
Tree Lengthe
(b)
Tree Lengths
40
VINCO
Density
80
VIMCO
30
60
20
40
10
20
0.46
0.48
Expected Mutations per Site
Expected Mutations per Site
(c)
Tree Log-likalihood
(d)
Tree Log-likelhhood
0.12
REAST
REAST
0.10
VGP110
VSP120
0.04
VBPR0
VBF
LOOR
LGOR
0.06
REP
VIMCO
Denaity
0.03
VIMGO
0.06
0.02
0.04
0.02
0.01
0.00
0.00
-7000
-6995
-6990
-6986
-6980
-6975
-6365
-6960
-8000
-7803
-7600
-7400
-7200
-7000
(e)
COVID-19 Dalaset
-7150
VBP11
Eslimated MLL (nats)
VEP00
-7175
LOOR
REP
-7200
-9000
VIMCO
Estimated
-7225
VSP11 VSP123
-250
LOOR
-7275
RE
VNCD
-11000
-7300
10'
10²
10
10"
lteration
Iteration
Figure 2. Variational inference results for DSI (left) and the COVID-19 dataset (right). (a Density estimation for tree lengths. (c- Density estimation for tree log-iikelsooods. Estimates are formed from 1,000 samples from variational posterior of each VI method and 97,500 samples from the BEAST gold standard. (e-f) Trace plots of estimated marginal log-likelihood vs. iteration number. The number of importance samples used to estimate the marginal log- likelihood was 500 for VBPI and 50 for VIPR.
from
VIMCO tended to produce trees with more variable log- likelihood values. Again, this is likely because the multi- sample ELBO results in mode-covering behaviour.
Macbook Pro with 16GB of RAM and a 2.6 GHz 6-core In- tel i7 CPU. The wallclock time in seconds for each method versus the number of taxa shown Figure 3, This also plotted in terms of log ratios in Appendix B.5 (with slope indicating complexity). Our method is approximately twice as slow as VBPI per iteration for 8 taxa, but it scales better and outperforms VBPI for 512 taxa. Even though evaluating the variational density of VIPR takes O(N2) time, VIPR has an empirical time complexity of roughly O(N), indicat- ing that the primary bottleneck is calculating the likelihood, which takes ((NM) time. We also demonstrated that the number of parameters for VBPI grows super-linearly with the number of taxa N (see Appendix B.5), so the asymptotic computational complexity of VBPI may also be super-linear.
VIPR converged quickly on DS1 because its parameters were initialized in region of high ELBO, while VIPR con- verged slower for the COVID-19 19 dataset since its parameters were initialized in a region of low ELBO. The optimization may have been caught relatively flat regions of the param- eter space, highlighting the need for intelligent parameter initializations or annealing schedules.
5.1. Computational Complexity
To compare the time complexity of our algorithm against VBPI, we considered the MS datasets (described in Sec- tion 4.4) with 1,000 sites and between 8 and 512 taxa. We ran each method for either 5 minutes or 1,000 iterations (whichever came first) and plotted the wall clock time per 1,000 iterations. These experiments were run on a 2019
For SBNs, as the number of taxa grows, the number of parameters grows with the number of trees in the SBN. There is closed form for this number-it depends on the MCMC and the posterior concentration. In Table 3 show
8

---

Variational Inference with Products over Bipartitions
Table Estimated evidence lower bounds for variational infer- ence methods (in nats). ELBOs were estimated using importance sampling on 1,000 random samples from each variational distri- bution. Our VIPR methods beat the VBPI baseline on half of the datasets. Dataset names, method acronyms, and conditions are the same those described in the caption for Table
Table 3. Number tree structure parameters versus number of taxa (NTAXA) on simulated data with 1,000 sites.
NTAXA
VBPI
VIPR
8
4
56
16
44
240
DATA
VBPI10
LOOR
REP
VIMCO
32
55
992
DS1
7157.99
7159.56
71599.5
7161.60
DS2
-26573.03
266569.56
26569.50
26570.77
64
3,826
4,032
DS3
33793.96
33794.96
733794.77
33796.53
128
29,939
16,256
DS4
-13541.39
13512.54
13512.60
13513.41
DS5
8281.03
8279.93
$280.35
8282.03
256
127,217
65,280
DS6
6751.77
6754.36
6755.29
6756.10
512
319,533
261,632
DS7
37331.12
37333.36
37332.04
37352.10
DS8
8657.78
8062.26
8661.8
8064.54
DS9
4088.64
4085.61
4087.25
4090.52
DS10
10111.81
10114.76
10113.16
10119.70
B030
DS11
6329.37
6289.60
6289.70
6294.31
1000 Itensiora
COV
8100.96
8489.82
8244.41
8087.43
8000
4000
the number of parameters in the SBNs for VBPI on each of the seven MS datasets. We compare VIPR, showing quadratic behaviour for VIPR and that VBPI uses more parameters when the number of taxa is greater than 128.
2000
6. Discussion
100
300
In this work, we introduce a new variational family over ultrametric, time-measured phylogenies that models the coa- lescent time between each pair of taxa. The family is formed by deriving closed-fore forrsdd-foooooo expression for the marginal distri- bution on phylogenies induced by single linkage clustering on a distance matrix. Methods using this variational family require only ((2)) parameters in total, and each param- eter has an intuitive interpretation as a description of the distribution on pairwise coalescents. For example, in this work place independent log-normal distributions the entries of the distance matrix, yielding 2(N) parameters (one mean and one standard deviation for each pair of taxa). VIPR is unique in that it does not require aspects of MCMC runs in its iterations in order to make inference computa- tionally tractable. In contrast, Zhang and Matsen IV (2024) used MCMC fix the support the trees described by their variational family. (Note that we also used short MCMC runs to initialize parameters.)
Figure 3. Seconds per 1,000 iterations vs. number of taxa. Each VI method was run for 1,000 iterations or 5 minutes (whichever took less) on simulated datasets.
by fixing the distribution on the time to coalescence of taxa u and v that are far away in genetic space at infinity, thus reducing the number parameters to learn. Expanding the variational family to include conditional parameters may also improve performance: if taxa and v coalesce first, we may define new parameter (((u,)),w) describing the coalesce between clade containing {u, and another taxon Further, fast and accurate parameter initializations and well-tuned annealing are essential for top performance in variational Bayesian phylogenetics. Our experiments may be improved by using an annealing schedule similar to Zhang and Matsen IV (2024) to prevent convergence to local maxima.
We have focused the difficult task of inferring tree topol- ogy and branch lengths. Inference for aspects such relaxed clock (see Douglas et al., 2021) and effective pop- ulations size can also be done starting from this new varia- tional family. Our method thus promising foundation on which more intricate variational families can be built.
Our methods may be further developed in many ways-for example, by moving from the normal distribution on pairwise coalescent times to mixture distributions Molén et al. (2024), or by using normalizing flows similar to Zhang (2020). We could also directly enforce sparsity in the prior

---

Variational Inference with Products over Bipartitions
Acknowledgements
Hunter, J. D. (2007). Matplotlib: A 2D graphics environ-
ment. Computing in Science & Engineering, 9(3).
This research was enabled in part by support provided by the Digital Research Alliance of Canada. We gratefully acknowledge all data contributors, i.e., the authors and their originating laboratories responsible for obtaining the speci- mens, and their submitting laboratories for generating the genetic sequence and metadata and sharing via the GISAID Initiative, on which this research is based. We also thank the anonymous reviewers for their helpful comments.
Jukes, T. and Cantor, C. R. (1969). Evolution of Protein
Molecules. Academic Press.
Katoh, K. and Standley, D. M. (2013). MAFFT multiple sequence alignment software version 7: Improvements in performance and usability. Molecular Biology and Evolution, 30(4).
Khare, S., Gurry, C., Freitas, L.,, Schultz, M. B., Bach, G., Diallo, A., Akite, N., Ho, J., Lee, R. T., Yeo, W., GISAID Core Curation Team, and Maurer-Stroh, S. (2021). GI- SAID's role in pandemic response. China CDC Weekly, 3(4)).
Impact Statement
This paper presents work that advances the field of machine learning. There are many potential societal consequences of machine learning (see for example Metz 2023). However, applications of our methods may also advance basic science in diverse fields that use phylogenetics such as epidemiology, linguistics and ecology.
Kimura, M. (1969). The number of heterozygous nucleotide sites maintained in a finite population due steady flux of mutations. Genetics, 61)).
Kingma, D. and Welling, M. (2014). Autoencoding vari- ational Bayes. In Proceedings of the 31 International Conference on Learning Representations.
References
Bengio, Y., Lahlou, S., Deleu, T., Hu, J., Tiwari, M., and Bengio, E. (2023). GFlowNet foundations. Journal of Machine Learning Research, 24(210).
Kingma, D. P. and Ba, J. (2014). Adam: A method for stochastic optimization. arxiv preprint 1412.6980.
Bouchard-Cóté, A. and Jordan, M. I. (2010). Variational inference over combinatorial spaces. In Proceedings of the 23rd Conference on Advances in Neural Information Processing Systems.
Kingman, J. (1982). The coalescent. Stochastic Processes
and their Applications, 13(3).
Koller, D. and Friedman, N. (2009). Probabilistic Graphical Models: Principles and Trchhiques-Adqpeee Computa- tion and Machine Learning. MIT Press.
Bouckaert, R. R. (2024). Variational Bayesian phyloge- nies through matrix representation of tree space. PeerJ, 12:e17276.
Koptagel, H., Kviman, O., Melin, H., Safinianaini, N., and
Lagergren, J. (2022). VaiPhy: A variational inference based algorithm for phylogeny. In Proceedings of the 36th Conference on Advances in Neural Information Pro- cessing Systems.
Douglas, J., Zhang, R., and Bouckaert, R. R. (2021). Adap- tive dating and fast proposals: Revisiting the phyloge- netic relaxed clock model. PLOS Computational Biology, 17(2).
Lakner, C., van der Mark, P., Huelsenbeck, J. P., Larget, B., and Ronquist, F. (2008). Efficiency of Markov chain Monte Carlo tree proposals in Bayesian phylogenetics. Systematic Biology, 57(1).
Felsenstein, J. (1981). Evolutionary trees from DNA se- quences: A maximum likelihood approach. Journal Molecular Evolution, 17(6).
L' L'Ecuyer, P. (1995). On the interchange of derivative and ex- pectation for likelihood ratio derivative estimators. Man- agement Science, 41(4).
Godoy, B. S., Camargos, L. M., and Lodi, S. (2018). When phylogeny and ecology meet: Modeling the occurrence of Trichoptera with environmental and phylogenetic data. Ecology and Evolution, 8(11).
Li, T., Liu, D., and Yang. Y. (2020). Phylogenetic supertree reveals detailed evolution of SARS-CoV-2. Scientific Reports, 10, 22366.
Hein, J.. Schierup, M., and Wiuf, C. (2004). Gene Genealo- gies, Variation and Evolution: A Primer in Coalescent Theory. Oxford University Press.
Linderman, S., Mena, G., Cooper, H., Paninski, I.., and Cunningham, J. (2018). Reparameterizing the Birkhoff polytope for variational permutation inference. In Pro- ceedings of the 2 International Conference on Artificial Intelligence and Statistics.
Hudson, R. R. (2002). Generating samples under a Fright-Fisher neutral model of genetic variation. Bioin- formatics, 18(2).
10

---

Variational Inference with Products over Bipartitions
List, J., Nelson-Sathi, S., Geisler, H., and Martin, W. (2014). Networks of lexical borrowing and lateral gene transfer in language and genome evolution. Bioessays, 36(2).
Shi, J., Zhou, Y., Hwang, J., Titsias, M., and Mackey, L. (2022). Gradient estimation with discrete Stein operators. In Proceedings of the 35th Conference on Advances in Neural Information Processing Systems.
Maclaurin, D., Duvenaud, D., and Adams, R. P. (2015). Autograd: Effortless gradients in numpy. In the ICML 2015 Workshop on AutoML.
Sibson, R. (1973). SLINK: An optimally efficient algorithm for the single-link cluster method. The Computer Journal, 16(1).
Metz, C. (2023). 'The Godfather of A.I.' leaves Google and warns of danger ahead. The New York Times.
Stamatakis, A. (2014). RAxML version 8: tool for phylo- genetic analysis and post-analysis of large phylogenies. Bioinformatics, 30(9).
Mimori, T. and Hamada, M. (2023). GeoPhy: Differentiable phylogenetic inference via geometric gradients of tree topologies. In Proceedings of the 37th Conference Advances in Neural Information Processing Systems.
Xie, T., Matsen IV, F. A., Suchard, M. A., and Zhang, C. (2024). Variational Bayesian phylogenetic inference with semi-implicit branch length distributions. arxiv preprint 2408.05058.
Mnih, A. and Gregor, K. (2014). Neural variational infer- ence and learning in belief networks. In Proceedings of the 31st International Conference Machine Learning.
Xie, T. and Zhang, C. (2023). ARTree: A deep autoregres- sive model for phylogenetic inference. In Proceedings of the 36th Conference on Advances in Neural Information Processing Systems.
Mnih, A. and Rezende, D. (2016). Variational inference for Monte Carlo objectives. In Proceedings of the 33rd International Conference on Machine Learning.
Xie, W., Lewis, P. O., Fan, Y., Kuo, L.,, and Chen, M.- H. (2010). Improving marginal likelihood estimation for Bayesian phylogenetic model selection. Systematic Biology, 60(2).
Molén, R., Kviman, O., and Lagergren, J. (2024). Im- proved variational Bayesian phylogenetic inference using mixtures. Transactions on Machine Learning Research, 3353.
Zhang, C. (2020). Improved variational Bayesian phyloge- netic inference with normalizing flows. In Proceedings of the 33rd Conference on Advances in Neural Information Processing Systems.
Nesterenko, L., Blassel, L., Veber, P., Boussau, B., and Ja- cob, L. (2025). Phyloformer: Fast, accurate, and versatile phylogenetic reconstruction with deep neural networks. Molecular Biology and Evolution, 42(4).
Zhang, C. (2023). Learnable topological features for phylo- genetic inference via graph neural networks. In Proceed- ings of the 11th International Conference on Learning Representations.
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang, E., Detito, Z., Rai- son, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, I.,, Bai, J., and Chintala, S. (2019). PyTorch: An im- perative style, high-ppeffrrance deep learning library. In Proceedings of the 32nd Conference on Advances in Neural Information Processing Systems.
Zhang, C. and Matsen IV, F. A. (2018). Generalizing tree probability estimation via Bayesian networks. In Pro- ceedings the 31 Conference on Advances in Neural Information Processing Systems.
Robbins, H. and Monro, S. (1951). A stochastic approxi- mation method. The Annals of Mathematical Statistics, 22(3).
Zhang, C. and Matsen IV, A. (2019). Variational Bayesian phylogenetic inference. In Proceedings of the 7th Inter- national Conference on Learning Representations
Rubinstein, R. Y. (1992). Sensitivity analysis of discrete event systems by the 'push out' method. Annals of Oper- ations Research, 39(1).
Zhang, C. and Matsen IV, F. A. (2024). A variational ap- proach to Bayesian phylogenetic inference. Journal of Machine Learning Research, 25(145).
Zhou, M. Y., Yan, Z., Layne, E.,, Malkin, N., Zhang, D., Jain, M., Blanchette, M., and Bengio, Y. (2024). PhyloGFN: Phylogenetic inference with generative flow networks. In Proceedings of the 12th International Conference on Learning Representations.
Sagulenko, P., Puller, V., and Neher, R. A. (2018). TieeTime: Mikeloodlikelihood phylodynamic analysis. Virus Evo- lution, 4(1).
Saitou, N. and Nei, M. (1987). The neighbor-joining method: a new method for reconstructing phylogenetic trees. Molecular Biology and Evolution, 4(4).

---

Variational Inference with Products over Bipartitions
A. Proof for Proposition 1
We provide a proof by induction for the derivation of Equation (7), Consider the coalescent events in Algorithm For 1 $ K S N I, let ti:k be the times of the first K coalescent events (tı:x = [t / Let TıKK be the bipartitions of the ). Let qe.k (e.k(TK, (x) be the probability density function of the marginal distribution the times and the bipartitions of the first K coalescent events. Let S, be the set {[w, =} w E W., & E Z (here (Wn. n, Zn is the n-th bipartition). Note that S, is the set of all unordered pairs of taxa that have not coalesced before bn and that coalesce f. Let S be U Sn (i.e., Si:K is the set of all unordered pairs of taxa that coalesce by time t.)) By the definition of Sn, sum En, ZweW..eEZ, is equal to the same sum indexed by {w, 2} Sn- (And the same is true of products.) Our induction hypothesis for I K S N 1 is as follows:
first K coalescent events (TI-K [[W,, Zn J m==))
K II ==1
(orse) (ln) go Qg fmm)) (tm)
l (T.KK,
II fw.s)ES.
Q(u=a) (tn)
II {xe==)$8..
Q Qlu (tx))
(tk).
(10)
fw.ztes
In Equation (10), the product over (w, 2} & SI.K is outside of the product over n = L,..., ......, K. (Throughout this derivation, a product has more than one factor in its operand, they are all enclosed by the pair of brackets appearing immediately after There exists an unordered pair of taxa {w", {=",2*} such that (W1; 3 ((w"), ". {==}). The probability density go.1 (TI:1; t1:1) ) is the density of the event that taxa w" and =* coalesce at time t (this density is "==""=") (t)) times the probability that all other taxa coalesce after time t (as {W1, Z1} is the first bipartition). Therefore, we have:
the product sign.) Consider the base case of the induction where
16.1 (T.a, t1a) 9
ww0 (00','}(11)
(11)
II (w,e)<(w";"))
Auw (h)
(11)
40
(w'
(t1)
stw"ts") (tt)
II (w.=)es
Awre (ti)
(12)
(h)
x
furzz) (ti) (11)
II (w.S)ES
Q Auuszt( (41)
II
(13)
Q (11)))(11))
(witEs
(wea)&Sia
Here in Equation (11) we use the mutual independence of (6->)) to split the joint probability of all taxa coalescing after 4 other than {w", }*}, into product. Thus, the base case (where K = 1) is established. Assume that the induction hypothesis in Equation (10) holds for a given - (here < K - 1< < N - 1). Consider the conditional probability density function K K ([WK, ZkY } tx| T1:K-1, t1: When we condition on T1:K t1:k-1, the K-th coalescent event with bipartition (WK, 3 occurs time 1 and only if the following hold:
1. There exists an unordered pair of taxa (w",2" 2²) E SK such that z(w" ,-) = tx. (We are conditioning on the event that taxa w* and == have not coalesced before time ( (x-).) The conditional probability density of the event that w" and =* coalesce at time tı is thus df
{w" =x") (fx)/Q( fw" ,-*) (tx-1))
(<r)'s are continuously 2. All other taxa pairs that have not coalesced by time f coalesce after time (x. (As the g differentiable, they are continuous and so {w*, 2²) is unique almost surely.) Note that we are conditioning on the event that all taxa pairs {w,, >} & S1.-11 have not coalesced before time f 1--1- The conditional probability density of the event that w" and z" have not coalesced by time ( is thus Q (w" ,2*) (lx) (w (IK-1).
The conditional probability density thus:
12

---

Variational Inference with Products over Bipartitions
(o.K(WWK. bx),
71K--1
(14)
(w' (''x='((tx)
Q
1 ((x)
Σ
II
Qa +
(15)
Quee" f/'" (1x-))
s (1
(w JESK
(w K-1
(m,=) (10
Awezt ((x))((k)) (Lx)
((x)
mm= (lx)
II
(16)
Qa fuuee) (tx)
Q (w,z) (1k-1)
(W.R)ESK
g guwre (tx)
Q (w)
(tk)
(wrz) (tx) o
(17)
2¹ fu0,0) (tk)
2
(tk-1)
tw,z) (tx-1)
(w,
(w.2)ESK
(w.a)&S..x
Note that we drop the stars the taxa and after Equation (15) because the indices no longer need to be distinguished once the sum is isolated. Also, in Equation (15) we use the mutual independence of g(x>) to form the product. Multiplying this conditional probability with the induction hypothesis Equation (10) for K - 1 yields the total probability density:
(o.x(T.x.t..k) = 96 ((Wx ZK},t tk TI:K-1, tixx-i)
4s.K- (
(18)
Lo
[1,22
(lx)
II
Q (1 (lx)
(wxz) (lx)
(19)
Q
(1m)) (tx)
(l (tx-1)
Qa (m,,2) (tx-1)
(w,
(w,,))
(w,=t&SS
K-1
1s (tn)
(bn)
II (w,s)sS.
Q ffu
II
Q Quw (1
(20)
Q (1) ((m)
(kn)
n=1
(w,
(w
We then move the first component of Line (19) into the first component of Line (20), and we move the numerator of the second component of Line (19) into the second component of Line (20). These rearrangements yield the following:
(ur,s) (1x)
(TiK(Ti tr:k)=
II (S.
II (w,,)&S
of (tx-1)
of (w,e) (tx-1))
K
<w) (tn)
II
2 ((u,))) (tm)
II (w,a)ES,
Q fm,=) (tn)
II
Q (w,=) (tx-1).
(21)
n=1
(w.aeES
(w,2)ES
For the final display, in Equation (22) split the numerator and denominator of of [w.=)) (lx) (/. into separate
products. And in Equation (23) we cancel the ((w=22 (tx-1) factors involving {w, 2) Sk. And in Equation (24)
cancel the (1122) (Ix-1) factors involving {we, =} & S1--1--
13

---

Variational Inference with Products over Bipartitions
(T.,K
II (1x
Quma)) (1x)
(mr))
Qluiaa (tx-1)
2 (urrtt (tx-
(w,=)&SLix
(m,,))&S1xx
(m)22) (tn)
II
II
Q ((x)))(tn))
II
Q
fw (1x-1)
(22)
Q
( (tn)
n=1
(w.s)ES.
(w.=)ES.
(w.e)EšL.x
Qlurae) (tk)
Q (f (tx-1)
(10,2) (1n)
II
II
Q (uua)t(t)))
II
ol (IK-1)
(23)
Q 1..2 (km)
n=l
(1,
(w,
[mm (tn)
Σ
II (u,
olussa)(1s) (lx)
of (lx)) {w,=)&S1:K
II
(24)
11 (kn)
(M,
Thus, the inductive step is established and Equation (10) holds for all < K < N - 1. To complete the derivation, note that N- = gø: and S = - S. = {(w, =} w, E } (all taxa coalesce after N coalescent events); and indices over w E Wn, 2 E Zn are equivalent to indices over {w, 2 E Sn.. Thus, for = N 1 the last term of Equation (24) an empty product yielding the desired result:
g0(T,t) 90 1
(25)
N-1
Lo (w,z) (Ir)
II
II (w.=tEs.
Qfuie) (tn)
(26)
Q (ux,))
(uui)) (tra)
(1)))
(w,=tES,
(w,=) (In)
II
Σ
II
(27)
( (fm)
W. n
WEW,, zEZn
2E
14

---

Variational Inference with Products over Bipartitions
B. Additional Results
B.1. Coverage of cubeVB
Our main comparison is between our method VIPR and the method of Zhang and Matsen IV (2024). Both these methods are based on optimization. However, another variational inference method for ultrametric phylogenetic trees described by Bouckaert (2024). We do compare this method in terms of likelihood estimation because is not optimization-basel Nonetheless, compare VIPR with Bouckaert (2024), constructed the maximum clade credibility (MCC) using gold standard BEAST run, selected ordering of taxa at random using the MCC tree, and calculated the percentage of tree topologies from the BEAST gold-standarr that are within the "cube space" implied by the ordering. This process estimates the percentage of the posterior that is impossible to reach using the restricted tree space from Bouckaert (2024):
Table 4. Percentage of trees sampled from BEAST that are outside of cube space defined by Bouckaert (2024). Percentages are listed for each dataset considered in our likelihood experiments. many cases, coverage of cube space is extremely small. Dataset sizes are provided below in Table 5.
DATA
OUTSIDE CUBE SPACE
DS1
29.2
DS2
15.2
DS3
76.8
DS4
79.7
DS5
98.0
DS6
94.7
DS7
69.9
DS8
42.7
DS9
99.9
DS10
84.6
DS11
99.9
COv
99.9
Some datasets correspond to posteriors where 99.9% of sampled trees lie outside of cube space. Although these results striking, the discussion of Bouckaert (2024) mentions that CB may struggle on high-entropy posteriors.
15

---

Variational Inference with Products over Bipartitions
B.2. Dataset Characteristics
To understand the characteristics of each dataset, we calculated pairwise Hamming distances between each taxa for each dataset (dropping sites with missingness). These dataset features may be cross referenced with the likelihood estimation in Section 5 to better understand relative model performance.
Table 5. Number of taxa N, number of sites M, and the average Hamming distance between sites for each dataset. Values in parentheses indicate standard deviations.
DATA
N
M
HAMMING DISTANCE / M
DSI
27
1949
.040(.01))
DS2
29
2520
.214(.057)
DS3
36
1812
.230(051)
DS4
41
1137
.138(.055)
DS5
50
378
.192(.04))
DS6
50
1133
.056(.029)
DS7
59
1824
200((069)
DS8
64
1008
.082(.01))
DS9
67
955
.025(.014)
DS10
67
1098
.070(.026)
DSI
71
1082
.082(.053)
COV
72
3101
.008(.003)
16

---

Variational Inference with Products over Bipartitions
B.3. Posterior Predictive Checks
In this Appendix, we run simple posterior predictive checks using 10 trees simulated from Kingman coalescent model with N 8 taxa and an effective population size of Ng 1.0. For each tree, simulated genomes using the lukes-Cantor model of evolution with M 1, 000 sites. We then ran VIPR using the leave-one-out reinforce estimator for 10, 000 iterations and the same optimization parameters as the primary experiments. Figure 4 displays histograms of sampling 10, 000 trees from the resulting variational distributions. Our method tends to underestimate tree length (the sum of all branch lengths) for trees of length or greater. This is likely because the posterior distribution for very long branch lengths approximately exponential, but we use log-normal distributions for pair-wise distances in VIPR. Besides that notable exception, the true tree length and ogg-liklihood comfortably within the posterior distribution estimated by VIPR.
Tee Lengths
Tree Log. licelihood
The Lengths
Tee Log llkelhhood
0.7
0.25
0.8
0.70
0.20
Doncity
0.2
0.15
0.15
0.10
0.10
0.05
0.00
-5725 -700 -88---88--888 -550
-7540
7830
Thee Longths
Tree Log- likelihood
Tree Lengths
Tree Log-ikkllhood
0.20
Denatt
1.0
0.15
0.3
0.15
0.10
0.10
0.1
0.05
3.0
3.5
4.0
4.5
77770
-7755
77760 7755 -7550
-5770
-5700
-5750
-574D
Tee Langths
Treel -oel--hllhoodd
Lengths
Log-likelihood
0.4
0.20
Denatt
0.3
0.30
0.1
0.00
0.5
0.05
-0740
-0730
-8720
2.0
3.0
-1230
-5240
-5230
55229
Lenitths
Log-licefitadd
Lengihs
Log-ikelinood
0.25
0.20
0.15
Denatt
0.3
0.10
0.05
0.05
0.0
0.00
0
0.00
Tie Lengths
Tiee Log-lioelitood
Tree Lengthis
Tree Log-litelitood
0.8
0.6
Dendit
0.15
6.15
6.0
0.00
0.0
0.00
-1620
-7000
-7040
-7080
-7000
-7640
Figure 4. Posterior predictive checks for VIPR. Checks include lengths (the of all branch lengths) and log-likelihoods. True trees were generated using N = 8 taxa and Kingman coalescent with N₂ = 1.0. Genomes were sampled using Juke-Cantor model of evolution and M 1, 000 sites. Results are shown for 10, 000 samples from the variational posterior (histograms) and the true tree used to generate the genomes (vertical black lines).
17

---

Variational Inference with Products over Bipartitions
B.4. Marginal Log Likelihood and ELBO Values
In this Appendix, include more complete versions of Tables and from the main text. In particular, Tables 6 and 7 below include the number of taxa and sites for each dataset, results for VBPI with batch size of K 20 (VBPI20), and estimated standard emrors using 100 bootstrapped samples. Further, Figure 5 below shows trace plots of marginal log-likell logllikelhhodd versus iteration number for all 12 datasets.
Table Gap between gold standard and estimated marginal log-likelihoods for variational inference methods (in nats). Marginal log- likelihoods for VI methods were estimated using importance sampling with 1,000 random samples from each variational distribution. Values indicate differences between gold standard MLLs and each method's MLLs. Gold standard MLLs (indicated by the BEAST column) are derived from independent chains of the stepping-stone algorithm in BEAST. Datasets (DATA column) DSI to DS1I are from Lakner al. (2008). Dataset COV the COVID-19 dataset obtoined from GISAID. methods are specified by columns: Variational Bayesian Phylogenetic Inference with K K-sample ELBO, K = 10 (VBPI10; Zhang and Matsen IV 2024); Variational Bayesian Phylogenetic Inference with K sample ELBO, K 20 (VBP120; Zhang and Matsen 2024); VIPR using leave-one-our REINFORCE estimator (LOOR); VIPR using the reparameterization trick (REP); VIPR using the Variational Inference for Monte Carlo Objectives estimator with (VIMCO). Standard errors were estimated using 100 bootstrapped samples and are shown parentheses.
DATA
(N,M)
BEAST
VBPII0
VBPI20
LOOR
REP
VIMCO
DSI
(27, 1949)
-7114.26(0.19)
-0.3())00)
0.36(0.13)
-2.29(0.15)
-1.83(0.21)
0..5(0046)
DS2
(29, 2520)
226566.42(0.222)
0.16(0.24)
0.01(0.20)
-0.76(0.14)
-0.14(0.43)
-0.37(0.29)
DS3
(36, 1812)
-33787.62 62((0.36)
0.41(0.12)
0.38(0.13)
-3.66(0.53)
-1.91(0.99)
2.63(00.50
DS4
(41, 1137)
-13506.05(0.3))
0.03(0.53)
0.46(0.43)
2.48(0043)
-0.47(1.21)
-1.23(0.23)
DS5
(50, 378)
-8271.26(0.99)
1..70(
5.69(0.48)
0.29(1.82)
-4.01(0.288
0.94(2.08)
DS6
(50, 1133)
6744.31
0.76(0.20)
0.32().35)
33.9600.34)
-3.26(0.60)
2.72(0.37))
DS7
(59, 1824)
-3733388800666
(..7().26)
00.24(0.17)
-2.73(0.30)
2.82(0.31)
10.42(00.70)
DS8
(64, 1008)
-8650.20000.0
- 0.82(0..7)
0.47(0.64)
-3.28(0.99)
-4.95(0.47)
-2.8800.60)
DS9
(67, 955)
4072.66(0.53)
5.32(0.31)
4.12(0.46)
-3.12(1.21)
5.79(0.74)
-.60000.4)
DS10
(67, 1098)
10102.65(0.63)
0.88(0.20)
1.44(0.22)
5.38(0.42)
-3.98(1.14)
6.82(0.49)
DS11
(71, 1082)
-6272.57
-18.79(0.41)
16622((0.46)
-6.79(0.89)
-7.31(0.71)
-9.62(1.46)
COV
(72, 3101)
-7871.61(0.74)
-39.08(0.58)
-33.26(0.76)
-611.8((L800
-374.62(0.48)
214.25(00.22)
18

---

Variational Inference with Products over Bipartitions
Table Estimated evidence lower bounds for variational inference methods nats). ELBOs were estimated using importance sampling on 1,000 random samples from each variational distribution. Our VIPR methods beat the VBPI baseline on half of the datasets. Dataset names, method acronyns, and conditions match Table Standard errors were estimated using 100 bootstrapped samples and are shown in parentheses.
DATA
MM
VBPII0
VBPI20
LOOR
REP
VIMCO
DSI
(27, 1949)
7157 (1.7.)
-7158.18(0.66)
-7159.56(0.10)
-7159.54(0009)
-7161.60(0.00)
DS2
(29, 2520)
-26573. (0.23))
-65733 (00)3)
26569.56(0.06)
26569.50 (0.69950(0008))
26570.7740013))
DS3
(36, 1812)
-33793.961 (0.8809600200)
-33794.75(0.28)
-37794
-33794.77(0.07)
-33796.5300.15)
DS4
(41, 1137)
-13541.391 (3541.39(33.12)) (33.12)
-13613.68
13512 (0.11)
13512.60
-13513
-13513.4100.14 (0..144)
DS5
(50,378) 378)
-281. (0.(0.26)) (0.26)
-8298.64 -898.66455977)
8279 (2.9.93(0.111)
-2880.35
8282.0030.177)
DS6
(50,1133)
-6711.77 77 (0.22)
-6752.60(0.21)
-6744366
(7755.9(00144)
6756.10(0.21)
DS7
(59, 1824)
37331 (.)22)
-37331
37333.36(0.19))
37332.04
37352.10(0.422)
DS8
(64, 1008)
8657. (0..30))
-8658.8300.22)
(662.26(0.16))
8661.88(0.16)
8664..4(00.26)
DS9
(67, 955)
---088.8400..0.
-4091.21
4085.61 (0.18))
-082722
-4090.52
DS10
(67, 1098)
11011181((0...
-10112.80(0
-10114.76(0.15)
-10115.1600.88
110119. 70(0.26)
DS11
(71, 1082)
63329.37(9.90
-6559.2 (554.24(511.99)
6289.60(0.17)
6289..0(00188)
6294.31
COV
(72, 3101)
-8100.96(10966)
-7933.440.998
8189.82(0.21)
-211.1100.21)
-8087.43(0.00)
19

---

Variational Inference with Products over Bipartitions
DS1
DS2
(nats)
Escimated MALL mnats)
-7175
-7200
-25700
-7275
-28800
10P120
LOOR
-20000
LOOR
VWCCO
VICCO
-7300
10'
10"
10"
10
10
10
10"
'0'
10"
lleratinn
heralion
DS3
DS4
D63
DS4
Estimated (natc)
-33800
Estimated MLL (nato)
-18900
-33000
-33900
-13000
15000
24400
19PG5
15000
107120
-4460
19400
14000
WWCC
10
DS5
Beration
DS6
herstoon
DS5
D96
Esiimaton MLL (yaat)
-2800
IEPTE
-7100
10'
10'
DS7
DS8
-07260
Estirated MLLL (yatt)
-8000
-37750
100000
-38000
-1000
-30000
LOOR
LOOn
-35750
12000
10'
10'
10'
t0'
16¹
10'
10
10'
DS9
honatinn
DS10
horation
DS10
-10000
-1000
Estinated MLL (nats)
100000
-10000
10
10
10
10
10'
10
DS11
Ne leration
COV
horation
DS11
COMID-P 19 Dataset
Estimated (mats)
Estirated NILL (nato)
-5002
10P10
-6400
V2P113
-10000
-8300
V8000
1977
-11000
Iteration
oe
Figure 5. Trace plots for all datasets. Trace plot of estimated marginal log-likelihood vs. iteration number (i.e., parameter update number). Marginal log- likeliiiehodd was estimated using 500 importance samples VBPI and importance samples for VIPR methods.
20

---

Variational Inference with Products over Bipartitions
B.5. Computational Complexity Results
We provide an additional plot for the results described on the MS datasets in Section 5.1. Computing the variational density of VBPI is linear in the number of taxa, but normalizing the SBN scales with the number of parameters. Therefore, the computational complexity of VBPI scales with the number of parameters. Figure below further illustrates this by plotting the slope of a log-log curve against the number of taxa, where the y-axis represents the computational complexity of each algorithm.
2.2
VBPI20
2.0
LOOR
VIMCO
Alcg Alog (ntaxa)
1.6
1.4
0
100
200
300
400
500
ntaxa
Figure 6. Slope of the logarithum of seconds-per-iteratin vs. the logarithm of the number of taxa. Each VI method was run rurfor 1,000 iterations on subsets of COVID-19 dataset. The y-axis corresponds to computational complexity of the algorithm as a function number of taxa (i.e., corresponds to linear complexity, 2 corresponds quadnatic complexity, etc.)
21

---

Variational Inference with Products over Bipartitions
C. Gradient Estimators for 96
C.1. The REINFORCE Estimator
Working from the definitions in Section 3.3.1, the leave-one-out REINFORCE estimator for VIPR is derived as follows:
V VOL(Ø) =
w
V @
log qsh (()))
(28)
K k=1
(k) = So(T (8), (k)) f((h))
(29)
f(-k) K
Σ fo(T!0) ((0)
(30)
kpk
(T(h), ,(()) 94
(31)
C.2. The Reparameterization Trick
We continue the derivation of the reparameterization trick for VIPR. Working from Equation 9 in Section 3.3.2, we proceed by summing the tree structures and then integrating over 2, (0), the space of all values of Z that are consistent with given the parameters Φ. This yields the following:
((0) = Σ "L 1) log
p ((Y,go(Z)) (5( (Z))
dZ.
(32)
Note that the region of integration (6) depends upon @, so interchanging the integral and the gradient introduces some error due to the Leibniz integral rule. Nonetheless, we proceed with the interchange and approximate the full gradient as follows:
p P(Y, 9ol (Z))
]
VoL(Ø) = Ez
V log
(33)
46 (Z))
Finally, we define a biased estimate of VoL(o) as the following:
g(,, g0 (Z(k))
(oL())
V log
(34)
K
20 Z (k))
(k)
N):
(35)
Y gy (2
can be calculated using automatic differentiation software such as Autograd (Maclaurin et al., 2015) or Py PyTochh (Paszke et al., 2019).
As with the LOOR estimator, the gradient V & log
(g (y)
C.3. The VIMCO Estimator
We derive the VIMCO Estimator for VIPR described in Section 3.3.3. For our model, the l-sample ELBO (Mnih and Rezende, 2016) is defined as follows:
p(T(M), t((), v((i)) (kk), t(t))
(Ø) = Egg
log
(36)
K
Here, (7(2), (()) ~ 90 for = 1,..., K. This is the objective function used by Zhang and Matsen IV (2024) to perform VBPI. When using the K K-sample ELBO objective from Equation (36), the VIMCO estimator is an analogous gradient estimator the LOOR estimator for the single-sample ELBOO, and is defined follows:
22

---

Variational Inference with Products over Bipartitions
K
VoLx(0) ~
- -w0(h) (0)
Vo log log (
(37)
k=
So(T(h) ((a)))
T(())
(38)
E (t(T
(6) i (0) = (x) - -log
Σ (T( t(0)) + 85-0)
(39)
K
Pyk
K
1
(k(0) = -log
(40)
K
k=1
p(-1)
K 1 Σ Ea(())) (00)
(41)
(T( ,t(()) 9s.
(42)
23


Reviews:
[{"review_id": "xDFWs07exJ", "reviewer": "Reviewer_en52", "summary": "This paper introduced a new method, VIPR, for phylogenetic inference. This new method greatly improves the computational efficiency without sacrificing accuracy compared with the traditional MCMC based method. The new method derives a closed-form density of the distribution over the entire tree space based on coalescent times and single-linkage clustering. This study proposed a new variational distribution based on coalescent time and single linkage clustering, which makes the computation more efficient. Experiments on benchmark dataset and one empirical dataset shows comparable accuracy and improved computational efficiency.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review: please find detailed responses below.\n\n*How do parameter numbers influence the time complexity of VBPI?*\n\nIn [Zhang and Matsen, JMLR2024], as the number of taxa grows, the number of parameters grows with the number of trees in the SBN. There is no closed form for the number of trees; it depends on the MCMC algorithm and posterior concentration. We empirically calculated the number of parameters in VBPI vs the number of taxa on datasets simulated with \"ms\" [Hudson 2002] with 1,000 sites.\n\n\\# of tree structure parameters:\ntaxa | VBPI | VIPR |\n---|---|---\n8 | 4 | 56\n16 | 44 | 240\n32 | 55 | 992\n64 | 3,826 | 4,032\n128 | 29,939 | 16,256\n256 | 127,217 | 65,280\n512 | 319,533 | 261,632\n\nComputing the variational density of VBPI is linear in the number of taxa, but normalizing the SBN scales with the number of parameters. We will add this experiment to the Appendix of the camera-ready copy.\n\n*Why should VIPR have a lower computational complexity?*\n\nIn our experiments, VIPR attains accurate marginal log-likelihoods estimates in fewer parameter updates than VBPI (Appendix B, Figure 5). VIPR has O($N^2$) parameters, and the number of parameters in VBPI can be larger than that if the SBN support is large (see Table above).\n\nWe improved our code using *scipy.cluster.hierarchy.linkage* and streamlined our phylogenetic likelihood function. We performed new speed comparisons for all methods on simulated datasets with varying numbers of taxa. See our response to Reviewer 3 for results.\n\n*What is the impact of log-normal branch-lengths? Any way to relax this?*\n\nAfter running BEAST, we plotted histograms of pairwise log-coalescent times across sampled trees for some of the datasets. In most cases these histograms looked normal, motivating our log-normal branch lengths. We will include a some of these histograms as a supplement. VIPR can incorporate any branch length distribution with continuously differentiable density. We will consider flexible branch-length distributions in future work.\n\n*Sequence divergence is not included.*\n\nWe calculated pairwise Hamming distances between each taxa for each dataset (dropping sites with missingness). Values in parentheses are standard deviations:\n\nDS | Hamming distance/\\#sites\n---|---\n1 | .040(.017)\n2 | .214(.057)\n3 | .230(.051)\n4 | .138(.055)\n5 | .192(.041)\n6 | .056(.029)\n7 | .203(.069)\n8 | .082(.031)\n9 | .025(.014)\n10 | .070(.026)\n11 | .082(0.053)\nCOV | .008(0.003)\n\nWe will add this Table to Appendix B in the camera-ready.\n\n*Include simulated datasets with better control on tree depth, sequence divergence, mutation rates, etc.*\n\nThis simulation is an excellent idea. We aim to do this in a follow-up paper.\n\n*Jukes-Cantor may be over-simplified vs. complex evolutionary models*\n\nWe agree that Jukes-Cantor is a simplified assumption. We aim to include K2P [Kimura 1980] and GTR [Rodriguez 1990] in a follow-up paper. (Note that Zhang and Matsen JMLR2024 only consider JC.)\n\n*Add more baseline methods like RAxML, neighbour-joining, and Vaiphy? The experiment does not cover uncertainty estimation.*\n\nThese methods do not apply specifically to variational inference over ultrametric trees. RAxML and Neighbour-joining do not provide estimates of the marginal likelihood. Vaiphy is for multifurcating trees. We have already mentioned VaiPhy and we will add more about non-Bayesian methods in the introduction:\n\n\"Phylogenetic inference can also be performed using non-Bayesian methods, including RAxML, Neighbour-joining, and Phyloformer. Phyloformer uses deep learning to construct pairwise representations of evolutionary distances between taxa. Phyloformer then uses pairwise distances to construct a  tree using a neighbor-joining algorithm similar to the method described here. Non-Bayesian methods do not provide estimates of marginal likelihood, which are useful for model selection.\"\n\nRegarding uncertainty estimation, we aim to add posterior predictive checks of tree length and clade support on simulated data to the Appendix for the camera-ready.\n\n*How would VIPR handle non-ultrametric trees?*\n\nA natural approach for non-ultrametric trees with our framework is to extend our strict clock models to relaxed clock models. Preliminary calculations suggest it is possible to do so at the cost of roughly twice as many variational parameters compared to the VIPR variational family.\n\n*How robust is VIPR when handling noisy or high divergent dataset?*\n\nVIPR struggled most with the COVID-19 dataset, where genomes varied relatively little (see Hamming distance Table above). We assume that highly divergent datasets would also be challenging. Future studies can quantify how VI methods such as VIPR are affected by noisy or divergent datasets.\n\n*What is the limitation of VIPR on dataset with 100+ taxa?*\n\nVIPR's empirical computation time per iteration is approximately linear in the number of taxa (see response to Reviewer 3). Future work can apply VIPR to larger datasets.", "claims_and_evidence": "This paper claims that the new method VIPR relaxes the dependency on MCMC subroutines and achieves better efficiency on phylogenetic inference. Despite lower computational complexity, this new method achieves comparable accuracy with the golden standard Bayesian phylogenetic inference methods. \nThe claims are supported by experiments on benchmark datasets and one empirical dataset, SARS-CoV-2, comparison across baselines, including BEAST, the gold standard MCMC-based method for phylogenetic inference, and approximate the true posterior distribution, and VBPI, a recent variational inference method for phylogenetic inference using subplot Bayesian networks but still rely on MCMC for tree sampling. \nThe experiment was conducted on 11 standard benchmark datasets with a wide range of taxa numbers and sequence length, and one empirical dataset, SARS-CoV-2, to evaluate real-world dataset with rapid evolving speed. \nThe results show that VIPR performs comparable with two baselines on MLLs and ELBO. The running time of the experiments shows that VIPR has a time complexity of roughly O(N^2). \nOne minor concern is that VBPI shows roughly the same computational complexity on empirical dataset. Could the author elaborate more on how does the parameter numbers influence the time complexity of VBPI? Why VIPR should have a lower computational complexity?", "methods_and_evaluation_criteria": "VIPR is proposed based on the previous studies of Bouckaert (2024), and Zhang and Matsen IV (2024) with improvements on the scalability and computational efficiency. \nVIPR does not rely on MCMC sampling like the other traditional methods. It directly models the distribution over the tree space. Compared with VBPI, this new method uses a variational distribution over distance matrix. This derives a differentiable variational distribution over the tree space, makes it possible to apply efficient gradient estimation for faster and more stable inference. \nVBPI directly optimizes the coalescent times/branch lengths, relaxes the limitation of to Bouckaert (2024) method matrix representation approach on the ability of tree representation. \nThe methods are evaluated on benchmark datasets and empirical dataset. The benchmark datasets covers a wide range of taxa numbers and sequence length, representing a range of complexity. \n\nVIPR assume a log-normal variational distribution. What could be the impact of this assumption? Any chance to relax this assumption to achieve better flexibility on the inference? \nAnother minor issue is that the sequence divergence is not included for the datasets. It would be helpful to get a rough sense of how difficult are those datasets and what is the impact on the method performance. \nThe author can also consider to include simulated dataset with better control on tree depth, sequence divergence, mutation rates, etc.  \nThe paper only considered Jukes-Cantor model, which may be over-simplified. Could the author consider more complex evolutionary models such as GTR?", "theoretical_claims": "This paper defines the tree space and shows its variational distribution covers the entire tree space. VIPR variational family enables gradient-based optimization. Proposition 1 shows a closed-form solution for the density function of trees. The probability density function over trees looks good. The derivations of gradient estimators look correct. \nThe theoretical claims are mostly valid and proofs look correct.", "experimental_designs_or_analyses": "The experiment design is reasonable. \nThe datasets covers a relatively wide range of difficulty levels. Including 11 standard benchmark datasets and 1 empirical dataset. \nThe method performance is compared against the MCMC-based golden standard method, BEAST, and a latest variational bayesian phylogenetic inference method, VBPI. \nThe metrics used for evaluation are valid. \n\nCould the author add more baseline methods such as faster heuristic methods like RAxML? \nThe experiment does not cover the uncertainty estimation.", "relation_to_broader_scientific_literature": "VIPR improves the previous methods introduced by Bouckaert (2024), and Zhang and Matsen IV (2024) with better scalability and computational efficiency. Current MCMC-based methods are limited to small datasets <100 taxas due to high computational cost. With the new method, it could enable Bayesian method on much larger dataset. \nThis new method also enables better integration with machine learning pipelines. Traditional phylogenetic inference methods are not differentiable. VIPR provides a differentiable method that is compatible with deep learning pipelines.", "essential_references_not_discussed": "To achieve a comprehensive landscape of phylogenetic inference studies, the author should consider to discuss other tree inference methods, such as maximum likelihood methods and distance based heuristics.", "other_strengths_and_weaknesses": "This new method shows good novelty. VIPR addresses an important limitation of previous methods, removes dependency of MCMC sampling. VIPR makes the Bayesian tree inference more scalable for larger datasets. \nThe paper is well-structured, with a clear motivation, theoretical proof, and experimental results.\nOverall, the paper shows a strong contribution to Bayesian phylogenetic inference method. \n\nThe paper is relatively simplified on the theoretical assumptions. Please consider to expand the methods to more complex evolutional model to better fit the real-world applications. The paper can be stronger by adding discussion on uncertainty estimation, comparison with other mainstream phylogenetic inference methods, such as RAxML, neighbor-joining, and Vaiphy.", "other_comments_or_suggestions": "I have listed my comments in previous sections.", "questions_for_authors": "I'm curious about the following questions:\n1. How would VIPR handle non-ultrametric trees? \n2. How robust is VIPR when handling noisy or high divergent dataset? \n3. What is the limitation of VIPR on extreme large dataset, for examples dataset with 100+ taxa?", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "0okuTKSJnU", "reviewer": "Reviewer_R16S", "summary": "This paper proposes a variational Bayesian phylogenetic tree analysis method using a matrix representation of tree structures. Phylogenetic tree analysis is one of the important analytical techniques used to estimate the developmental process and diffusion pathways of a target, and is more and more in demand in formulating future preventive measures, for example, for recent infectious disease pandemics. Conventional Bayesian phylogenetic tree analysis faces several challenges. One is the use of Markov chain Monte Carlo methods in many models and algorithms, whose efficiency, both theoretically and empirically, is not yet clear. Another is that many of them do not explicitly include the fusion time in the phylogenetic tree in their models, and as a result, those methods cannot properly capture the ultrametric nature of the tree structure. As a way to solve these two problems, this paper proposes a model and its inference method that can properly reflect ultrametric by explicitly modeling coalescence time using a tree structure representation that has high affinity to variational methods. The effectiveness of the proposed method is demonstrated with 7 data frequently used in many recent Bayesian phylogenetic analyses and the more practical SARS-Cov-19 data.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We appreciate the thoughtful suggestions below, and hope that we have addressed your comments sufficiently.\n\n*One minor concern is that a similar tree-structured matrix representation has been studied independently in another paper [Bouckaert2024] very recently.*\n\nAs we discuss in our literature review, we have two primary significant improvements compared to [Bouckaert2024]: 1) we use optimization-based VI to maximize the ELBO, and 2) we provide a closed-form density in Proposition 1. Using this density means we do not have to restrict the tree space. In addition, our framework allows us to use a relatively general class of variational distributions for the pairwise distances, whereas Bouckaert require log-normal distributions in order to incorporate a covariance matrix $\\Sigma$.\n\n*Would it be difficult to see the reduction in performance when restricted to only restricted trees within the framework of the proposed method?*\n\nWe believe that this comment refers to selecting an ordering of taxa similarly to [Bouckaert2024], and then re-running VIPR while restricted to the \"cube space\" from [Bouckaert2024] consistent with the ordering. This is an excellent idea to investigate the effect of restricting tree space. \n\nWe have now constructed the maximum clade credibility (MCC) tree from BEAST using our gold standard MCMC run, selected an order from the MCC tree, and then calculated the percentage of tree topologies from the BEAST gold standard that are within the \"cube space\" implied by this ordering. This process estimates the percentage of the posterior that is impossible to reach using the restricted tree space from [Bouckaert2024]:\n\nDS | \\% of MCMC trees outside cube space\n---|---\n1 | 29.2\n2 | 15.2\n3 | 76.8\n4 | 79.7\n5 | 98.0\n6 | 94.7\n7 | 69.9\n8 | 42.7\n9 | 99.9\n10 | 84.6\n11 | 99.9\nCOV| 99.9\n\nThese results are striking, but [Bouckaert2024] mentions that CubeVB may struggle on high-entropy posteriors in their discussion. We will add this Table and discussion to the Appendix of the camera-ready.\n\n*It is not easy to intuitively understand what improvement VIPR has over the variational distribution of the standard mean-field approximation.*\n\nOne of the key challenges for variational inference over phylogenetic trees is that using a mean-field approximation is not straightforward. There are two main reasons for this.\n\nFirst, we can only apply a mean-field approximation after decomposing the distribution of the tree as a product over cliques of random variables. There is no standard way of doing this for trees. In Matsen IV (2024) this is done by forming a subsplit Bayesian network with one node per subtree appearing in the MCMC samples used to initialize the support. Our novel decomposition (Proposition 1) is another way of decomposing the distribution of the tree as a product over coalescent times. This results in $\\mathcal{O}(N^2)$ parameters, improving upon the worst-case performance of Matsen IV (2024), in which the number of parameters could be super-exponential in the worst case. We then use our novel decomposition for the mean-field approximation (using state-of-the-art techniques for gradient evaluation and optimization: autograd, VIMCO, REINFORCE and the Reparameterization Trick).\n\nSecond, for ultrametric trees we cannot assume independence between coalescent times, lest the resulting tree violate the ultrametric constraint. To overcome this challenge, we form our our variational family to approximate the matrix of pairwise coalescent times (the matrix bold T). We then map T to ultrametric trees using single-linkage clustering.\n\n*A quantitative comparison with another related study might have strengthened its persuasiveness.*\n\nWe chose VBPI for a baseline comparison because it was the only VI-based method for ultrametric trees that we are aware of in the literature. We did not include [Bouckaert2024] because it does not rely on optimization, so we could not include it in our trace plots of marginal log-likelihood vs iteration number. As you have suggested, restricting our method (and the BEAST gold-standard) to the same restricted tree space as [Bouckaert2024] is a valuable experiment to isolate the effects of optimization versus unrestricted tree space, and we aim to complete this experiment in a follow-up paper.\n\nThank you for catching the title formatting error&mdash;we have now corrected it.", "claims_and_evidence": "The main claim of this paper is that using a matrix representation for the representation of tree structures has two benefits: (1) ultrametric measures can be captured and (2) a differentiable variational representation that avoids MCMC sampling, whose theoretical analysis and empirical goodness for mixing time is not yet clearly known.\n\nOne minor concern is that a similar tree-structured matrix representation has been studied independently in another paper [Bouckaert2024] very recently. However, the authors, in fairness, identify differences and improvements over prior work in Section 2.5.", "methods_and_evaluation_criteria": "This paper uses a variational phylogenetic tree representation that utilizes a matrix representation of the tree structure and derives an inference algorithm using three choices of loss functions. Through experiments, the proposed method is compared to previous state-of-the-art phylogenetic tree analysis methods (including the most related and recent one [Zhang&Matsen, JMLR2024]). The evaluation criteria used are the marginalized likelihood and ELBO in terms of learning and prediction performance.\n\nOne minor concern to me is whether the ablation study made it difficult to quantitatively assess the improvement of the proposed method from the other related study, [Bouckaert2024]. Would it be difficult to see the reduction in performance when restricted to only restricted trees, as in Literature A, within the framework of the proposed method?", "theoretical_claims": "The key theoretical contribution of this paper is that through the matrix representation of the tree structure, the variational distribution obtains an easy-to-handle closed-form representation, as described in Proposition 1. \n\nThe key theoretical contribution of this paper is that through the matrix representation of the tree structure, the variational distribution obtains an easy-to-handle closed-form representation.\nI may not yet properly understand the empirical benefit of this theoretical result; as shown in Appendix 1, I can see that this representation does indeed lead to an easy-to-handle closed-form expression. On the other hand, it is not easy to intuitively understand what improvement this has over the variational distribution of the standard mean-field approximation. Any help from the authors in this regard would be greatly appreciated.", "experimental_designs_or_analyses": "Experiments are conducted on seven datasets that have been used expressively as benchmark data in recent Bayesian phylogenetic tree analyses, as well as on the SARS-Cov-19 data for more practical applications. The evaluation by the marginalized likelihood and ELBO also reflects recent trends, and the experts feel that the improvement in performance is clearly reported.\n\nThe comparison method seems convincing enough, as it is extremely up-to-date. On the other hand, as discussed in the “Methods” section, a quantitative comparison with another related study A might have strengthened its persuasiveness.", "relation_to_broader_scientific_literature": "As evoked by the recent pandemics, phylogenetic tree analysis is one of the machine learning tasks that has received particular attention in recent years. This paper is not intended to bring any new insights from a scientific point of view, but this technology is expected to contribute to the development of computational biology through the development of general-purpose machine learning.", "essential_references_not_discussed": "This paper provides a comprehensive discussion of Bayesian phylogenetic tree analysis, from its historical development to the latest developments in recent years. In particular, the relationship between the proposed methods and the challenges and room for improvement are carefully discussed in fairness to recent related research.", "other_strengths_and_weaknesses": "(Editing)", "other_comments_or_suggestions": "I was just a little concerned as to whether the title conforms to the format specified for the conference. Considering the margins involved in some of the figures, this is not overly space-hungry and may not be a problem for the draft stage for peer review. However, it may be a good idea to have it corrected in the camera-ready version if accepted.", "questions_for_authors": "One minor concern to me is whether the ablation study made it difficult to quantitatively assess the improvement of the proposed method from the other related study, [Bouckaert2024]. Would it be difficult to see the reduction in performance when restricted to only restricted trees, as in Literature A, within the framework of the proposed method?\n\nThe key theoretical contribution of this paper is that through the matrix representation of the tree structure, the variational distribution obtains an easy-to-handle closed-form representation.\nI may not yet properly understand the empirical benefit of this theoretical result; as shown in Appendix 1, I can see that this representation does indeed lead to an easy-to-handle closed-form expression. On the other hand, it is not easy to intuitively understand what improvement this has over the variational distribution of the standard mean-field approximation. Any help from the authors in this regard would be greatly appreciated.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "eHp91SKwkK", "reviewer": "Reviewer_uJ3o", "summary": "This paper targets variational inference of ultrametric phylogenetic trees and proposed a method called VIPR.\nAlthough many efforts have been paid in the field of machine learning based variational phylogenetic infernece, very few researchers consider this on ultrametric trees.\nVIPR sample a ultrametric phylogenetic tree by executing single linkage clustering on a distance matrix which is learnable and parametrized by log-normal distribution.\nA main contribution of this paper is the density formula of the abovementioned tree distribution (Proposition 1), and this allows training VIPR with gradient-based methods.\nThe authors validate the effectiveness of VIPR on DS1-11 and Cov-2 benchmarks and compare VIPR to the sota method VBPI.\n\n## update after rebuttal\n\nFirst of all, I would like to thank the authors for their detailed response that addresses my questions well. I have one further question on the convergence speed of VIPR. The author shows that it converges faster than other baselines (e.g., VBPI). However, this is a design choice in VBPI that uses the anneal schedule to encourage exploration over the tree space. Or in other words, there is a trade-off between the convergence speed to a good likelihood value of trees and the coverage of the high posterior region. I would expect the convergence speed of VIPR to slow down a bit with a similar annealing schedule, but end up with a better coverage over posterior trees (this would potentially improve the performance of VIPR in terms of marginal likelihood estimation). It would also be interesting to show how close is the approximate tree topology distribution provided by VIPR to the ground truth posterior from MCMC.\n\nOverall, I really like the idea of constructing a variational distribution over ultrametric trees from a distribution of pairwise distances through single-linkage clustering. It would be interesting to explore how to make it flexible enough for complicated tree posterior, which seems a bit challenging, as it is nontrivial how the correlation between pairwise distance would translate into the correlation between tree topologies. That said, I have updated my review accordingly.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your insightful review: see our responses below.\n\n*For inference accuracy, VIPR's trees' likelihoods lagged behind VBPI. For inference speed, no results report computation time.*\n\nThe aspect where VIPR shines is in time (or number of parameter updates) to attain an approximation error, as shown in Figure 5 of Appendix B. For example, on DS1 VIPR-LOOR achieves an estimated marginal log-likelihood (MLL) of -7175 after \\~100 iterations, while VBPI10 takes \\~1000. For DS5, VIPR-LOOR achieves an MLL of -8300 after \\~200 iterations, while VBPI10 takes \\~5000.\n\nFigure 5 compares iterations, and we agree that we should also report computation time. Thus, we ran all VI methods on simulated datasets with varying numbers of taxa for 1,000 iterations and reported computation time (see our response to Reviewer 1 for the simulation procedure). \n\nSeconds/1,000 iterations:\n\ntaxa | VBPI10 | VBPI20 | LOOR | REP | VIMCO\n---|---|---|---|---|---\n8 | 24 | 45 | 55 | 74 | 55\n16 | 43 | 88 | 110 | 150 | 112\n32 | 94 | 192 | 234 | 314 | 240\n64 | 192 | 381 | 475 | 633 | 473\n128 | 500 | 913 | 1,016 | 1,383 | 1,018\n256 | 1,560 | 2,395 | 2,150 | 2,952 | 2,162\n512 | 6,958 | 8,780 | 5,014 | 6,822 | 5,060\n\nOne iteration is one parameter update. Our method is \\~ twice as slow as VBPI per iteration for 8 taxa, but it scales better and outperforms VBPI for 512 taxa. We will add the results to Figure 4 and the Appendix. \n\nWe improved our code since submission (see our response to Reviewer 1). VIPR's primary computational bottleneck is now the phylogenetic likelihood, which takes $\\mathcal{O}(NM)$ time.\n\n*How is VBPI used for ultrametric trees?*\n\nZhang and Matsen IV [ICLR 2019] applies to general additive phylogenetic trees, but the follow-up paper [Zhang and Matsen IV, JMLR2024] extends the approach to ultrametric trees in Sections 6 and 7. The github repository (https://github.com/zcrabbit/vbpi-torch) contains code for ultrametric trees (in the directory \"rooted\") that we use in our experiments. We will make this more clear in the manuscript by providing a Github repository link and referencing the sections within the JMLR paper.\n\n*GeoPhy should be considered*\n\nGeoPhy is similar to VIPR in that it uses a tree construction algorithm on a distance matrix, but it is used for unrooted trees while we focus on ultrametric trees.\n\n*Figure 3(c) (VBPI is better) contradicts with Table 2 (VIPR-VIMCO is better).*\n\nFigure 3(c) reports the marginal log-likelihood, but Table 2 reports the ELBO. For the COVID dataset, VBPI is better in marginal log-likelihood, and VIPR-VIMCO is better in ELBO.\n\n*ARTree should be discussed.*\n\nWe will change Section 4.2 (line 264, column 2): \"We compare VIPR to the VBPI algorithm as implemented by Zhang and Matsen IV (2024), which uses MCMC runs to determine likely subsplits in an SBN.\"\n\nWe will also add ARTree to the introduction (line 26, column 2): \"For example, ViaPhy (Koptagel et al., 2022) uses a gradient-free variational inference approach and directly sample from the Jukes and Cantor (1969) model, GeoPhy (Mimori and Hamada, 2023) uses a distance-based metric in hyperbolic space to construct unrooted phylogenetic trees, and ARTree (Xie and Zhang, 2023) uses graph neural networks to construct a deep autoregressive model for variational inference over phylogenetic tree structures.\"\n\n*Phyloformer is similar to VIPR*\n\nWe will mention Phyloformer in the introduction: see our response to Reviewer 1.\n\n*The title page violates the ICML format*\n\nThank you for catching the title formatting error. This arose due to a copy/paste mistake. As Reviewer 2 noted, our Figures are not too space-hungry, and we could resolve this by cutting down the length of the Results and Discussion as you suggest, and improving the location of the Figures.\n\n*Ultrametirc trees should be (re)defined*\n\nWe will adopt the suggested definition (line 74, column 1):\n\n\"We focus on ultrametric trees, in which the leaves of the trees are all equidistant from the root. We denote our ultrametric trees with a rooted, binary tree topology tau and a set of coalescent times ...\"\n\n*The $N_e$ in the prior distribution is unclear*\n\nWe expanded section 2.3:\n\n\"We use the Kingman coalescent (Kingman 1982) as the prior distribution on the trees. This coalescent process proceeds backward in time with exponentially distributed inter-event intervals, and coalescent events occurring at rate $\\lambda_k = \\binom{k}{2}/N_e$, where $k$ is the number of taxa and $N_e$ is the effective population size, a parameter which governs the rate at which species coalesce. We fix $N_e = 5$ in our experiments. At each coalescent event, a pair of taxa are chosen to coalesce into a single taxon uniformly at random over all pairs of taxa. ... \"\n\nWe fixed the \"PhylogGFN\" typo, thank you for catching this.", "claims_and_evidence": "Although this paper presents a sound methodology for VIPR, its potential advantages (e.g., inference accuracy or speed) are not revealed by the experiments.\nFor the inference accuracy, the likelihoods of VIPR's inferred trees lagged behind those from VBPI, as suggested by the Figure 3 (b,c) and Table 2.\nFor the inference speed, no results directly reports the computation time of different methods.", "methods_and_evaluation_criteria": "The proposed methods and/or evaluation criteria make sense.", "theoretical_claims": "I've read the proof for Proposition 1 but not checked it carefully.", "experimental_designs_or_analyses": "- As the variational family in VBPI is applied to the general class of additive phylogenetic trees, the authors should clearly explain how they use VBPI to infer ultrametric trees.\n- GeoPhy (NeurIPS 2023; https://arxiv.org/pdf/2307.03675) considered a NJ algorithm on distance matrix for constructing phylogenetic trees, similarly to VIPR. It should be considered as an important baseline in terms of inference speed and accuracy.\n- Figure 3(c) (VBPI is better) contradicts with Table 2 (VIPR-VIMCO is better).", "relation_to_broader_scientific_literature": "Many prior works considered variational phylogenetic inference, e.g., VBPI, GeoPhy, but few of them considered inferring ultrametric ones.\nThis paper makes contribution in this sense.", "essential_references_not_discussed": "- The claim \"The VBPI baseline requires MCMC runs to determine likely subsplits (i.e., evolutionary branching events).\" in Sec 4.2 is not accurate. A notable progress of VBPI is ARTree (NeurIPS 2023; https://arxiv.org/abs/2310.09553), which does not rely on subsplits and should be discussed.\n- Phyloformer (https://www.biorxiv.org/content/10.1101/2024.06.17.599404v1) constructs a phylogenetic tree with a neighbor joining algorithm on pairwise representations. This idea is similar to VIPR.", "other_strengths_and_weaknesses": "Weaknesses:\n- The title page violates the ICML format (one-column titles and missing author information). I suggest the authors cut down the length of Results and Discussion to creating some space for the title page.\n- I think ultrametirc trees should be defined by \"the leaves of the trees are all equidistant from the root\", and the authors' definition in Sec 2.1 seems somewhat misleading.\n- The authors does not clearly explain the $N_e$ in the prior distribution in Sec 2.3.", "other_comments_or_suggestions": "There is one typo, Line 32: PhylogGFN.", "questions_for_authors": "I have no other questions fo the authors.", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This article propose a variational inference approach to the challenging problem of Bayesian phylogenetic inference. The key idea is proposing a variational family related to trees from hierarchical clustering, and because the density admits a simple closed form, it avoids much of the costly computation that existing methods are subject to. Reviewers all lean positive and while they did not discuss after some prompting, my own read agrees with their positive sentiment and find that this is a quality paper and the ideas will be interesting to the ICML readership.
Decision: Accept (poster)