Improved Variational Bayesian Phylogenetic Inference using Mixtures
Authors: Ricky Molén
Oskar Kviman
Jens Lagergren
Venue: Accepted by TMLR ()

Abstract: We introduce VBPI-Mixtures, an algorithm aimed at improving the precision of phylogenetic posterior distributions, with a focus on accurately approximating tree-topologies and branch lengths. Although Variational Bayesian Phylogenetic Inference (VBPI)—a state-of-the-art black-box variational inference (BBVI) framework—has achieved significant success in approximating these distributions, it faces challenges in dealing with the multimodal nature of tree-topology posteriors. While advanced deep learning techniques like normalizing flows and graph neural networks have enhanced VBPI's approximations of branch-length posteriors, there has been a gap in improving its tree-topology posterior approximations. Our novel VBPI-Mixtures algorithm addresses this gap by leveraging recent advancements in mixture learning within the BBVI domain. Consequently, VBPI-Mixtures can capture distributions over tree-topologies that other VBPI algorithms cannot model. Across eight real phylogenetic datasets and compared to the considered benchmarks, we show that VBPI-Mixtures result in lower-variance estimators of the marginal log-likelihood and smaller KL divergences to an MCMC-based approximation of the true tree-topology posterior.

Published in Transactions on Machine Learning Research (11/2024)
Improved Variational Bayesian Phylogenetic Inference using Mixtures
Ricky Molén
rickym@kth.ss
KTH Royal School of Technology Science for Life Laboratory
Oskar Kviman"
okviman@kth..se
KTH Royal School Technology Science for Life Laboratory
Jens Lagergren
jensl@kth.s
KTH Royal School of Technology Science for Life. Laborotory
Reviewed on OpenReview: https: /ttps://openpeen opw.revee.....eee.oo net/ forum? FBLM+
Abstract
We introduce VBPI-Mixtures, algorithm aimed at improving the precision of phyloge- netic posterior distributions, with focus accurately approximating tree-topologies and branch lengths. Although Variational Bayesian Phylogenetic Inference (VBPI) a state-of- the-art black- box variational inference (BBVI) framework has achieved significant success in approximating these distributions, faces challenges in dealing with the multimodal na- lure of tree-lopology posteriors. While advanced deep learning techniques like normalizing flows and graph neural networks have enhanced VBPT's approximations of branch-length posteriors, there has been gap in improving its tree-topology posterior approximations. Our novel VBPI-Mixtures algorithm addresses this gap by leveraging recent advancements in mixture learning within the BBVT domain. Consequently, VBPI-Mixt ures can capture distributions over tree-topologies that other VBPI algorithuss cannot model. Across eight real poylogenetic datasets and compared to the considered benchmarks, we show that VBPI- Mixtures result vorirnvariance estimators of the marginal log- likelihood and smaller KL divergences to an MCMC-based approximation of the true tree-topology posterior.
1 Introduction
Phylogenetic inference has wide range of applications in various fields, such as molecular evolution, epi- demiology, ecology, and tumor progression, making it essential tool for modern evolutionary research. Bayesian phylogenetics allows researchers to reason about uncertainty in their findings about the evolution- ary relationship between species.
The posterior distribution over phylogenetic trees given the species data is, however, challenging to infer, since the latent space is a Cartesian product of the discrete tree-topology space and the continuous branch- length space. Furthermore, the cardinality of the tree-topology space grows as a double factorial of the number of species (taxa), making the marginal likelihood computationally intractable in most interesting problem settings.
For over two decades, Markov Chain Monte Carlo (MCMC) approaches have been the go-to approaches for Bayesian phylogenetic analysis, where the MrBayes software (Huelsenbeck & Ronquist 2001 has been
"Equal contribution.

---

Published in Transactions on Machine Learning Research (11/2024)
<
(a)
(b)
(c)
Figure 1: Visualization of samples from the tree-topology posterior using a 1,000,000,000 iterations long MCMC run on (a) DS4, (b) DS7 and (c) DS8. Nodes represent unique tree-topologics and are colored basod cluster assignments, illustrating the multimodality of the tree-lopology posterior. The size of a node is determined by its sampling frequency, which is why nodes with low frequency appear black. More details in Sec.
particularly popular. However, random walk Metropolis- Hastings MCMC methods (Huelsenbeck & Ronquist,
2001 Höhna et al., 2016) rely on local operations to explore the trec-topology posterior, limitation which is known to require long MCMC runs in order to visit posterior regions which are separated by low-probability
tree topologies (Whidden & Matsen IV 2015). Sequential Monte Carlo methods for Bayesian phylogeny
2015 2020) have been proposed to avoid these local operations, but the resampling mechanism can filter out important trees in early steps of the algorithm as well as cause degeneracy, necessitating many particles.
[Bouchard-Cöte et al.) 2012) Wang et al.
More recently, variational inference (VI) has been applied to Bayesian phylogenetics. In general, VI is often promoted over sampling-based approaches in high- nigh-dimensional problems as a variational approximation of the posterior is obtained from optimization, making VT less vulnerable to the curse of dimensionality. However, in practice, it can be challenging to do VI without utilizing sampling. For example, in Koptagel et al. (2022) coordinate-ascent update equations are derived in the phylogenetic setting, but these are evaluated using importance sampling.
In blx VI (Ranganath et al. 2014), the gradients are instead taken of Monto-Carlo integrated csti- mates of the objective, lypically the evidence lower bound (ELBO), using samples laken from the variational approximation. This approach has been successfully applied in the Variational Phylogenetic Bayesian In-
ference (VBPI; Zhang & Matson IV (2019) framework, along with its extensions (Zhang) 2020 2023). In these extensions, more complicated branch-length approximations have been proposed using normalizing flows (NFs; Rezende Mohamed (2015) and graph-neural nets (GNNs; Kipf & Welling (2017)).
BBVI allows the practitioner to learn posterior approximations without deriving update equations or closed- form gradient formulations, the samples are taken from distribution that is commonly known to concentrate on high-probability regions of the posterior, resulting in a learning procedure that does not sufficiently explore the posterior distribution. This was addressed for continuous distributions in Ruiz et al. (2016), where samples were instead drawn from an overdispersed proposal distribution. In discrete hicrarchical models, similar to the targets approximated by VBPI, however, insufficient exploration may result in low-level states not being properly modclod. As the trco-topology posterior known to typically be multimodal with many "subpeaks" (Whidden & Matsen IV 2015) (visualized in Fig. [ it is thus of significant importance to encourage the tree-topology approximation to explore the posterior.
We propose VBPI-Mixtures, novel combination of two recent advances, VBPI and mixture learning in BBVI (Kviman et. al.) 2023). The mixture components cooperatively explore the tree-topology posterior during learning, and the increased flexibility of the VBPI-Mixtures results in approximations that can model posteriors intractable for the vanilla VBPI (sce See. 3). Using a toy experiment where we design a compli- cated hicrarchical categorical target distribution and pick a variational family such that it does not contain the target distribution, we first show that the mixtures of these sub-optimal variational distributions achieve
2

---

Published in Transactions on Machine Learning Research (11/2024)
smaller Kullback-Leibbe divergences to the targets and smaller importance-weigltt variances than single- component approximation which uses more Monte Carlo samples to estimate the objective. We then apply VBPIIMixtures on eight popular real datasets, visualize the learned tree-topology approximations and com- pare numerically with MCMC "golden runs", illustrating the joint exploration of the tree-topology space by the mixture components. Additionally, we derive a mixture-specilic formulation of the VTMCO estimator MMiih & Rezendo 2016) variance reduction technique using control variates to estimate the gradient of the importance weighted ELBO (Burda et al. 2015 Our contributions can thus be summarized as follows:
We propose VBPI-Mixtures, a novel algorithm for Bayesian phylogenetics.
We show that mixtures of subsplit Bayesian nets (SBNs) can approximate distributions that a single SBN cannot, making persuasive case for VBPI-Mixtures.
We derive the VTMCO estimator for mixtures.
We visualize two-component mixture of SBNs on real data, verifying that the components jointly explore the tree-topology space.
Experimentally, we achieve stato-of-the-art results on eight popular real phylogenctics datasets in lerms producing lower-variance estimalors of the marginal ligelikel,oood, showing that mixtures of SBNs can accurately capture tree-topology posteriors that non-mixture SBNs cannot, and produce KL divergences to MCMC-based approximation of the true trec-topology posterior that decrease with S.
2 Background
Let ee F(r)} denote the set of branch lengths for topology, т and F(r) the set of edges in т. Furthermore, the data, X = (X.,... [x} MMNN are an observed sequences of length M on the N leaves of the phylogenetic tree. Each entry in Xm, is character in the alphabet Ω, e.g., Xm,n (A,C, G, T) if DNA sequences are considered. The posterior distribution over leaf-labeled phylogenetic trees,
p(B, x|) P[X|T, (),,B))po (BB)),
(1)
(XX
is intractable due to its normalizing constant, p(X), i.e., the marginal likelihood.
Although p(X) is intractable, the three terms in the generative model in Eq. | can be computed efficiently. Typically, P(T) is uniform distribution over the (rooted or unrooted) tree-topology space, and the branch- length prior is an exponential distribution with rate X, such that p(B]r) = Ile de -Mb(e) The likelihood, p(X B), can be evaluated in lincar (in N) time using the standard dynamic programming algorithm proposed by Pelsenstein 2003).
Let a clade be a non-empty subset of the set of the N leaf labels, X. A subsplit a partition of this clade into two lexicographically orderod, disjoint clades, while a split simply a root subsplit-a bipartition of X. Furthermore, primary subsplit pair (PSP) is a subsplit conditioned on split-a tripartition of X. In Appendix A additionally give brief introduction to Bayesian phylogenetic inference for machine learning researchers.
Variational Inference in Bayesian Phylogenetics The Vl-based approach to Bayesian phylogenetics is to approximate T() using a simpler distribution, (B, = qe (B)r)a.()) Generally, in VI, the approximations are learned by maximizing the evidence lower bound (ELBO),
((X) = =E (B,.) log p(X|T, B)x( (B
(2)
implicitly minimizing the reverse Kullback-Leibler (KT.) divergence from p(B,7|X) to (B,
3

---

Published in Transactions on Machine Learning Research (11/2024)
Subsplit Bayesian Networks Given set of tree topologies, T it is straightforward to form look-up table of all subsplits in T. The SBN uses the look-up table to define support over possible tree topologies and learns the probabilities of the subsplits in the table via, for example, stochastic optimization. As the look-up table contains probabilities of subsplits, it is referred to as a conditional probability table (CPT). When the CPT has been learned, the SBN provides a tractable probability distribution over tree topologies from which it is possible to sample. See Zhang & Matsen TV (0188) or Zhang & Matsen IV (2024) for the original, more in-depth accounts of SBNs.
Regarding the pre-computed set of tree topologies, T, that are given input to the SBNs, it is worth noting that, although T is a subset of the set of all possible tree topologies, this should not be considered as a major constraint on the support of the SBNs. Support for this argument has been provided in [Zhang & Matsen IV (2024) (see Fig. 5 in the reference), where they showed that the number of pre-computed trees is close to the total number of topologies in the support of a. 10 billion MCMC (MrBayes (Huelsenbeck & Ronquist 2001)) coverage of UFBoot close to 100% of the MrBayes topology posterior.
run. Also, the
VBPI In VBPI ZZang & Matsen IV 2019 the posterior approximations are learned by maximizing multi-sample (Burda et al. 2015 version of ((X),
1 K log
p(X|T* )p(5p 17 pp(+)) (B ) (T+)
((X:K) K) = Egg , no. (8.r)
(3)
K k=
where (T( is SBN with learnable CPT, o, and (BIT) is multivariate LogNormal distribution with a diagonal covariance matrix, such that
qu (BIT) = II (( eEB(T)
(4)
Two different parameterizations of p(e, T) and Te, have previously been proposed. The simpler approach is to let u(e,r) = -/// and Tee,T) = / where e/T denotes a split of I in edge e. The parameters v/r and v/e/- are shared among all tree topologies where c/T exists, resulting in amortized mapping from a tree topology to the parameters of go (B))). Additional local information about the given T can be added into the paramctcrization of the approximation by using PSPs,
ple, r) = VefT
Σ of ( = =V//T + Σ R kEefT
(5)
iEefT
where e // r denotes the set of PSPs neighboring to the split e/T, and of is a learnable parameter associated with the i-th pair.
Multimodality of the tree-topology posterior In Whidden & Matsen IV (0015) modes are referred to as clusters of MCMC samples that are densely grouped in the tree-topology space and have high posterior density compared to their neighbors. They proposed a method for detecting and quantifying peaks by calculating the reversible subtree pruning and regrafting distance between topologies, combined with the MrBayes MCMC posterior probability. By applying this method, it was identified that cortain datasots had a high number of modes (e.g. DS1, DS4, DS5, DS6, DS7), which shows the complexity of the tree-topology space. In Fig. we visualize the multimodality of the posterior on three datasets.
Mixtures in Black-Box VI Learning mixtures of approximations in (black-box; Ranganath et al. (2014) VI Nalisnick et al., 2016 Morningstar et al., 2021 Kviman et al. 2023 is compelling off-the-shelf tech- nique to increase the flexibility of a variational approxination. Mixtures can be applicd to any variational approximation, including an NE-based approximation or one for discrete latent variables, with little overhead.
The objective function, the ELBO for mixtures, is estimated by sampling from each mixture component in a stratified manner (Monningstar al. 2021), or via multiple importance sampling techniques (Kviman et al.
1T is in practice obtained from some efficient tree-topology sampling algorithm, typically UFBoot Minh al. 2013
4

---

Published in Transactions on Machine Learning Research (11/2024)
A
B
A
B
B
A'
B
2
3
5
6
3
5
4
5
2
3
6
Figure 2: Mixtures of SBNs increase the flexibility of the tree-topology approximation. For instance, they can exactly capture a target distribution that assigns its probability uniformly across Ti and 72, leaving zero probability to 73 and T1. Meanwhile, this is not possible for a single SBN. See Example 3.1 for details.
2022 why the objective is often reforred to as MISELBO. Maximizing MISELBO encourages the mixture components to cooperatively cover the larget distribution, which is thought to be the key ingredient for their success in density estimation tasks (Kviman al.) 2023). In the next section, formulate MISELBO for VBPI and explain why mixtures of SBNs are beneficial for exploring the complex trec-topology space.
Importantly, the MISELBO objective shares similarities with the objective used in the mixture of experts (MoE) literature (Sh et al. 2019), however, note that in MoEs the mixture components are conditioned on different data types (e.g., images or texts), whereas here the data used for conditioning is identical for all components. As such, there a crucial distinction between the MISELBO objective and the one from the MoE literature.
ARTree In Xie & Zhang (2023): the authors employ an autoregressive sampling where actions are sampled of where a Icaf should be inserted into graph starting with three leaves. When all lcaves in the taxa are inserted they have sampled one topology. The model utilizes fixed order of insertions such that the topology can only be achieved with single unique order of decisions of where the insertions of leaves happen. Additionally, ARTree is an alternative to SBNs by leveraging GNNs to encode topological features, thereby eliminating the constraint of pre-sampled topologies inherent to SBNs.
PhyloGFN PhyloGFN Zhou et al. 2024 employs deep generative flow networks Bengio et al. (2023) to sample lopologies based on a set of actions of joining leaves and subtrees. The generative flow network setup makes possible to handle the non-ume unique set of actions to sample the topology as opposed to the constraint used in Xie & Zhang (2023). This model constructs phylogenetic trees using a bottom-up approach with transformer networks and is not constrained by pre-sampled topologies.
3 Variational Bayesian Phylogenetic Inference using Mixtures
Here we present our proposed method, VBPI-Mixtures. We derive the VIMCO estimator for learning mixtures of SBNs, and show how to combine mixtures of branch length models with an expressive NF model. We start by providing the MISELBO formulation for VBPI with importance samples,
$
p(X B )p(p )p)(?"") K
1 Σ
((X:KSS)
Egis .1 wow. (B.77) log K
(6)
S
k= 1 Det (Bk Jgo, , (T
To evaluate L(X; S)S,, approximate the s-th expectation by Monte-Carlo integration using simulations
from Ie..d. r). Although there exist works where the mixture component weights are learned (Roeder
et al. 2017 Morningstar et al. 2021), , we constrain ourselves to uniform mixture weights as in Kviman et al.
2023
Mixtures promote exploration Note that maximizing Eq. 6 w.r.l. the variational parameters will pro- mote the components to diversify and explore regions the latent space that correspond to regions measured by p. That is, samples drawn from component that result in high values in the numerator will indicate that the s-th component covers a part of the latent space corresponding to high-probability region, attracting the other components to this region, while the denominator penalizes similarities among
5

---

Published in Transactions on Machine Learning Research (11/2024)
the components, promoting the other components to spread out and find other important parts of the la- tent space. This exploratory behavior is crucial in black-box VI as the ELBO is only evaluated in sampled (visited) latent variables (states).
More specifically, the samples are proposed from the same distribution that attempts to maximize the ELBO. Consequently, for 5 = 1, there is a risk that, for a fixed K, less probable regions will not be sufficiently explored during learning and thus will be poorly modeled. In fact, [Zhang & Matsen IV (2019) showed that the vanilla VBPI does not necessarily benefit from more importance samples during training in terms of resulting marginal log- likelihood scores. Fortunately, the MISELBO objective offers a promising solution, as the mixture components are promoted to spread out and efficiently explore the multimodal phylogenetic posterior.
Mixtures increase the flexibility of the approximations A mixture of LogNormal pdfs is clearly more flexible than a single LogNormal pdf. For mixtures of SBNs, this is also true but may be less clear. An SBN constructs a tree by stochastically partitioning the available clades. The partition of a clade is sampled independently of the partitions sampled in the other clades. Mixtures of SBNs allow for modeling correlations in the sampling of the partitions, and thus increase the fiexibility of the approximation. We explain this feature with simple example, which trivially generalizes to larger trees and also applies to unrooted trees.
Example 3.1 Consider the four rooted free topologies in Pig. where A, A' and B' are sour different
subtrees for the Iwo clades (in red), and target densily, such that p(T)) = p(T2) = 0.5. A unāformly weighted minture of two SBNs can easily approzimate exactly by letting q (I(TI) = 1 and (2( (T=(T))=11 = resulting (2(2)) (=22==0 = 0.5. However, as a single SBN, 9, samples A or and or B' independently, and in order to achieve q(T)), 4(22) > O, will have to assign a non-zero probabilily Specifically, say ((A)= = a and, consequently, ((A)) aa while BB) = B and, consequently, q(B)) 1-3. It follows that q(ra) = a(1 - B) and q(Ta) (1 15--jj a)B. Finally a3 q(Ti) ((T2) (1 a))113)) implies a + B, which, in turn, implies that q(r3) = a² and q(Ta)= = B2. That is, all four trees will either have probability 1/4 under q, or one of T3 and 74 will have a higher probability than Tn and 72. So, single SBN can only yield distribution that is very different from p.
in (11) 0+
all four sublrees.
The above example exemplifies that there are tree-topology distributions that cannot be modeled using a single SBN, but can be modeled by a mixture of SBNs. There is no converse example, as a mixture can trivially model a single SBN by letting qi (T( = (T) for all r. In Appendix B we construct an example that shows that these conflicting tree-topology posteriors can indeed occur for real DNA data.
3.1 VIMCO for Mixtures of Tree-Topoloy Approximations
Here we derive the VIMCO estimator of Eq. G Although the notation in this section is specific Bayesian phylogenetics, our result is applicable to any mixture approximation. We purposefully follow the derivations in Mnih & Rezende (006)) closely.
3.1.1 Gradient Analysis
The gradients of Eq. 6 are studied first. Let
p(Bbk xhe X) 1 E Ovis (B3)) )gs (To,(TT)) (Th)²²
f(x, 13th (b)
(7)
8
==
and LKK log R Ek f(x, Ba,,B,, +,), where By T are simulated from de. (B,T). Note that sαx Bk + and so also i.K. a function of the SBN parameters for all mixture components, i.c., P1, &s. However, we omit these as arguments to the function in order to avoid cluttered notation. Furthermore, the sets of SBN parameters are disjoint (L.e. the mixture components do not share weights).
6

---

Published in Transactions on Machine Learning Research (11/2024)
We are interested in the gradient of Eq. components, say 1
6 with respect to the SBN parameters for of the mixture
V (( Ve
Σ
Eg. (8,7) [
(8)
==1
The full derivations are given Appendix C and lead to the following expression
0: (()
5 gonno (8,+) (B,r)
L
Φe log qo. (r+)
(9)
s
K
s Σ j=1
$ qv. (B,77) (B,7)
w
V
4e log
S
90 (B (35111)000 Jao, (T=)
We make three important observations, (i) for S == 1, retrieve the gradients used to train VBPI, (ii) the second term is negated, maximizing it promotes component i + s to be dissimilar from component s, diversifying the mixture. Finally, (dii) the first term is merely scaled (by L/S) version of the corresponding term in the $ = 1 case. Connecting to obscrvation (dit), we conclude that extending the VIMCO estimator to S > 1 cannot be trivially achieved without our derivation provided above.
Furthermore, the analysis of the gradients of the importance weighted lower bound in (Mnih & Rezende (016)-using notation, L(X; ()):K)>< K)-applics here, too. That is, the gradients in the second term in Eq 9 are multiplied by normalized weights, ensuring that the norm of the weighted sum over all K gradients is not greater than the norm of the largest term in the sum. This means that o; will be updated mainly according to gradients based on simulations scored highly by f, while mitigating the impact of gradients from lower-scoring simulations.
In the first term, on the other hand, all gradients are multiplied by the same i, indicating that the gradients of high-scoring simulations will not receive more weight than low-scoring ones, causing high variance and slow learning.
3.1.2 The VIMCO Estimator
As concluded above, the second term in Eq. (9) well-behaved, and we do not need variance-reduucion techniques to use it for learning in practice. The first term, however, requires attention in order to facilitate efficient learning.
Fortunately, as the first term is merely a scaled version of the corresponding term in the VIMCO estimator
when = 1, we can directly apply the localized learning signal strategy from Mnih Rezende (016 to obtain the VTMCO estimator for 521 1,
V ((X;KK S)
kk [()
Vo,
log qo. (7:"
(10)
5
K
s
W V
log 5
Σ Ivs (B, (B% Mo, (T!))
S
3=1
where Th B ~ ao..s (T, and T,, B ~ vo (T, B). Here, i - is the local learning signal for sample k defined as
1
Σ (2,77 BB8 f (,, ,,,))) k'tk
i [1) LK log K
(11)
with
x. T B- being an estimator f (x, The BB), typically the geometric mean Mnih Rezende 2016 Zhang& & Matsen IV (2019); Zhang (0220). f (ar -. BBT)) exp ( log TŞ ,B
In Appendix D we show that there is no bias introduced the VIMCO estimator by letting the variational distribution be uniform mixture distribution.

---

Published in Transactions on Machine Learning Research (11/2024)
4 Experiments
In Sec. we argued that a single-component approximation will struggle to properly model all parts of the target distribution when learned with black-box VI. Below, in Sec. 4.1 experimentally verify this claim, and, furthermore, confirm that mixture components collaborate in order to jointly cover the target density, resulting in more accurate approximations and efficient exploration.
We then, in Sec. 1.2 demonstrate that the increased model flexibility and promotion of exploration translates into lower-variance estimators of the marginal lielliiel and smaller KL divergence to an MCMC- -Csee approximation of the true tree- veo-topo posterior. We also visualize representations of VBPI-Mixtures on real data. Code for all experiments provided al. Github
4.1 Exploring Discrete Two-Level Hierarchical Model using Black- Box VI
Here, examine how mixtures explore discrete hicrarchical target distribution when learned via black- box VI. We construct this experiment lo analyze the capacities of discrete mixture models versus non- mixture models when the variational family of each component does not contain the target distribution. The experiment is thus aligned with Example 3.1.
Using a two-leveled hierarchical model of categorical distributions as the target distribution, p(z1, 22) = ((22)2)))(21, we infer the parameters of the variational approximations by minimizing KT. 1 E do, (21, 22) (2) i.e. the reverse KT. The CPD p(22|=1) is categorical distribution with ne categories, conditioned the sampled category in the previous level, 21, and p(21) is categorical dis- tribution with in categories. The family of variational approximations assumes independence between and 22, such that the s-th component is defined as qo, (21, = gøz (22) (z1), where oz are the learnable category probabilities over the n2 categories of component s. Hence a single component is insufficient to capture the posterior accurately. Meanwhile, a mixture of these approximations will be flexible enough given a sufficiently large S.
The parameters of p are drawn from Dirichlet distribution with all concentration parameters equal to 0.5, and the learning rates when inferring variational distributions are chosen based a grid search a separate target distribution. All o, are initialized uniformly over the categories. We compare an S = 10 mixture model which uses Z = 1 Monte Carlo estimates to estimate the objective function, with an S = 1 and Z = 10 model. This is to equate the number of samples used during learning.
We learn the approximations using either == or K 10 importance samples. That is, when K the distributions are inferred by minimizing the reverse KT., and when K 10 they are inferred by minimizing the importance-woightee reverse KL. The accuracies of the distributions inferred with K = are measured in terms of the forward KL (see Fig. 3). As an additional metric of accuracy, the distributions inferred using K = 10 are evaluated also based on the variance of their importance weights (shown in Fig. 4bb. That is, we sample 10000 importance weights, wi = p[ 2/)/,2))/3 E, (s) where (+), ~ 3 Σ. 9o. (2f,2)) and compute the variance of these weights. Note that the total number of samples used is the same for all models (if S then = 10, and vice versa). The variance metric assesses the aptness of employing the distributions as proposal distributions (lower is better).
Regardless of choice of learning rate, the non-mixture models are not able to accurately capture the posterior. The mixture models can, on the other hand, accurately approximate the target distribution. In figures 3 and we report the mentioned metrics over iterations when n1 = 5 and n2 = 10. All models are trained for 60k epochs, with decaying lcarning rates by 0.9 every 10k epochs.
4a 4b
4.2 Posterior Approximations using Real Data
We performed experiments on eight datasets (Hedges et al. 1990
Garey et al. 1996
Yang & Yoder 2003
Henk et al.] 2003) Lakner et al.) 2008 [Zhang & Blackwell 2001
Yoder & Yang
2004
Rossman et al.
2001) which we will refer to as DS1-8. These are popular datasets for evaluating Bayesian phylogenetics
methods, and, as in Zhang& & Matsen TV (2019) Zhang 2020); Moretti et al. 2021
Koptagel et al. (022
8

---

Published in Transactions on Machine Learning Research (11/2024)
m= n1=10
1.2
1.0
5=1 (<=1,2==1011,mitt0..00000
0.8
S=1 (K=1,2=10, lr_init==00000555
5=1 (K=1, |r_Init=00000001255))))
0.6
5=10 (K=1, 2=1, Ir_init=0.1)
0.4
5=10 (K=1, 2=1, Ir_init=0005))
5=10 (K=1, Z=1, lr_nit:0.02)
0.2
0.0
10000
20000
30000
40000
50000
60000
Figure Forward KL divergences from the target visualized defined in Sec. 4.1 the mixture approx- imations with S 10 and Z Monte Carlo samples to estimate the objective function, and to single component approximations with Z = 10. On the z-axis are the number of training iterations. When using higher initial learning rates than those shown in the plot, the KL-curves of the single-component approxi- mations diverged.
n n_1=5,n_ n_2=00
n_1=5, n_2=10
5-
S=1 (K=10, Z=10, ___ini==000000 0005)
S=1 (K=10, Z=10, Ir l_ii==0000005)
S=1 (K=10, Z=10, Ir fr_int==000000 00025)
S=1 [K=10, Z=10, Ir. l_li==00000055
1.2
S=1 (K=10, Z=10, ir_init=0.000111
S=1 (K=10, Z=10, Ir k__nit=00000125)
5=10 [K=10, Z=1, Ir_ 11
4-
5=10 (K=10, Z=1, Ir. 1_init=0..1)
1.0
S=10 (K=10, Z=1, Ir 1_init=0.05)
S=10 (K=10, Z=1, Ir v_init=0.05)
5-10 (K-10, Z-1, Ir_ (_iinitt00055)
5-10 (K-10, Z- Ir. 1__jnnit0.022
0.8
3
0.6
2-
0.0
0-
10000
20000
30000
40000
50000
60000
0
10000 20000 30000
40000
50000 60000
(a) Forward KL divergences.
(b) Imporlance weight variance with 10000 importance weights.
Figure 4: Evaluations of the accuracies of the approximations with respect to the target distribution described
in Sec. 4.1
Zhang 2023); Zhang & Matsen IV 2024), we focus on learning the approximations of branch-length and cree-topology distributions. Following the referenced works, we assume the exponential branch-length prior has rate 10 and a uniform prior over all unrooted trecs (sec Soc. 2 for details about the generative modcl). The substitution model is the Jukes- Cantor 69 model (Jukes al, 1969), and the candidate trees. T, are gathered from ten replicates of 10000 ultrafast maximum likelihood bootstrap trees (Minh et. al. 2013 The

---

Published in Transactions on Machine Learning Research (11/2024)
(a) DS5
(b) DS6
(c) DS8
Figure 5: Visualization a uniformly weighted S 2-component mixture of SBNs on (a) DS5, (b) DS6 and (c) DS8, where each node corresponds unique tree-topology. The upper row shows the distribution of five million sampled tree topologies from the first component, where node, т, is colored blue if qo. (r) > qa, (T), orange otherwise. Vice versa for the lower row. The size of node is determined by its sampling frequency, which is why nodes with low frequency appear black. The components clearly spread out, exploring different parts of the space.
implementation is based on the code provided by Zhang & Matsen IV 2024), and we trained all VBPI models during 400,000 iterations, using the same hyperparameter settings as Zhang & Matsen IV (2019): Zhang (020)). Based on the study in Zhang & Matsen IV (2024), we let K = 10 during training. Additionally, 3BNNbasen models use the same pre-computed topology sets, T.
4.2.1 Visualizing the Explorative Behaviour of Mixtures of SBNs
To graphically confirm the power of employing mixtures of SBNs, in Fig. 5 visualize representations of the learned tree-topology posterior approximations for subset of the different real datasets. The subset selected based on the datasets where the explorative behavior of the approximations is most clear. The representations for the other datasets, along with more implementation details, are included in the Appendix The SBNs correspond to VBPI-Mixtures without NFs. The components (upper vs. lower row) have jointly explored the space, partly specialixing on disjoint sets tree-topologies, verilying that the MISELBO objective promotes coordinated exploration of the discrete latent space, as discussed in Sec. 3
E.2
4.2.2 Quantitative Evaluation of the Approximations
We quantitatively evaluate the approximations by, first, computing their KL divergences to the true posterior, and by, secondly, benchmarking their marginal log-likelihood estimates. The results are averaged over five independently trained models with different parameter initializations. All MrBayes (Huelsenbeck & Ronquist 2001) results were produced using ten million long MCMC runs with four chains, sampling every 100 iterations. Our methods are denoted Mixs, MixNF,s, representing VBPI-Mixtures with PSPs or NFs, respectively. Mixtures that employ NFs share flow models, as described in Kviman al. (023).
Statistical distances the trec-topology posterior Here, we compare statistical distances the tree-topology posterior obtained from MrBayes as described above. In Table ! the KL divergence from the posterior to the approximations, i.e., KT.(((7|X)]lee(()) is computed, where 4s represents mixture of
10

---

Published in Transactions on Machine Learning Research (11/2024)
Table 1: Illustrating the impact of VBPI-Mixtures in terms of (L(p(T/X)]laa Lower is better.
DS1
DS2
DS3
DS4
DS5
DS6
DS7
DS8
VBPI-NF
0.0726
0.0110
0.0540
0.2093
2.2117
1.2842
0.2544
0.6018
Mix
0.0631
0.0059
0.0475
0.0965
2.0337
1.0883
0.1183
0.5199
Mix
0.0598
0.0051
0.0377
0.0769
1.9526
1.0461
0.0847
0.4567
Negative marginal log-likelihood estimates on DS1-8, All VBPI methods use 1000 importance samples, and the results are averaged over 100 runs and three independently trained modcls. All models are trained with one Monte Carlo sample (Z) and ten importance samples (K), except for those models with names notated with specific Z- or K-values. Following Zhang & Matsen IV (2019); Zhang (2020)
Table 2:
2023): ;
Zhang & Matsen IV (2024). we bold font the results with the lowest standard deviations (shown in
parentheses). Details in Sec.
4.2.2
Mixtures monotomically improve with S.
Data
DSL
DS2
DS3
DS4
DS5
DS6
DS7
DS8
Taxa
27
29
36
41
50
50
59
64
# Sites
1949
2520
1812
1137
378
1133
1824
1008
VBPI with Mixtures
VBPI
7108.50
06367.70(0.00) 09) 33735,
13330,0310. 23) (214.80(0.30)
6724,59(0.53)
97332.12(00.55
8652.39(0.7)
VBPIK=a0
7123.98(0.19)
26367.71
33735.09 (0(0088) 13330.07(0.222 (0. 8215.01(0.599)
6729.16(0.60)
37334.43(0.4) (0.46) 8654.73(0.61
VBPTz-3
7108.92(0.20)
26367.70(0.11
23735.09(0.1)
13330.066 (600000(0000 24) 8214.72(0.45)
(0.5.27(0.50)
87332.15(0..
5653.40(0.52
Mixs=2
7108. 44(0 12)
26367.
06) 33735. 10(0.07) 13330 00(0. 17) 8214. (0.36)
6721.54(0.31)
37332 022.2.010..22 21) 6651.68(0.11)
Mixs=s
7108.42(0..1)
26387.71
33735. (0(0.06) 13329. (7) 8214.73
(0.285)
37332.03 18) 8650.83
VB VBPI with NF's and Mixtures
VBPIxr
7108. 15)
76367.72(0
33735.10(0.07 10(0.07) 13330 00000.000323
7201070(00.1
6721.0(0 15)
37332 01 27) 8650 (0.46)
VBPINE =a
7108.37(0.16)
26367.7U(0.05. 33735.1100..06 13329, 9010.19) (30.3)
6724.38(0.32 32)
37331.9200 92....920.220 26) 8650. (0..52))
VBPINE.
7108.37(0.17)
26367.70(0.03)
33735.09(0.04
13329.94 94..111)) (00..11)221.55.. 6214.
67.3.37(0.39)
37331-9(00. (9. (0.18) 8650.51
MixgF,s-2
7108.40;0.10)
(6.67
13735.10(0.05 13329. 95 15) 6214.62 (0.26)
6724.44(0.322)
37331 (6( 19) 8650.56(0.39
Mix
7108.40(0.09
26367.70(0.03) (3533509900.04) 13329.94 (1.1))44.01))...338.3 (214.56(0.22)) 6724.40(0.25)
37331.96(0.2) $650.54
Scores from Zhang & Matsen 1V (2019), Zhang (2023), Xie & Zhang (2023), Zhou et al. (2044
MrBayes..
7108.42(0.18.
26367: 44Q060)) 13630 (0.-8)
6721.0700
80) 37332 rot (2 42) 8049 18,088811777 1.75)
GGNN
7108. 10(0.19) 26367. (3007, 10) 33735.11 (0.099)))))1 (3.1999)( 8211.67(0.36)
6721.3800
12)
37332. (3.32.03(0.300 30) 8650 68 18)
EDGE
7108.41(0.144) 26367.73(0.07)
33735.12(0.. 13329 13329.94(0019)) 94( (214.64(0.388)
6724.37
40)
37332.0440.22)
8650.65 (5.0655(0.45)
ARTree
(1(0.1 66367.71(0.07)
38735.09(0.09)
13329 94(0.17) 8211.59
672133
16)
37331
6650.61(0.00 60.488)
PhyloGFN
7108.95(0.06)
26.66.90(0.26)) 33735.60(0.35) 13331 83(0.19) (2.5.15(0.20) 6739.68(0.54)
37359.96(1.14) (601.7)
SBNs, or a single SBN, from VBPI-NF. Lower is better, and, notably, VBPI-Mixtures consistently produce KT, divergences across all datasets that monotonically decrease with S.
Marginal log-likelihood estimates In terms of marginal log- likelihod estimates, we benchmark our methods against the existing VBPI algorithms: VBPI with PSP parameterization (Zhang & Matsen IV
2019), VBPI-NF (Zhang) 2020) with ten ReaINVPs (Dinh et al. EDGE and GGNN). Additionally, we compare our results with the stepping-stone (SS; Xie et al. 2011
2016). and VBPI-GNN [Zhang] (2023
method applied to MrBayes as well as to PhyloGFN (Zhou et al. 2024) and ARTree (Xie & Zhang 2023
The results are given in Table 2 All models have been trained using a single Monte Carlo sample and ten
importance samples except for those marked with Z and K respectively. Following Zhang & Matsen IV (2019) Zhang (2020 2023): Zhang & Matsen IV 20024), we bold font the results with the lowest standard deviations. Rewarding low-variance estimates is motivated, as they imply, for instance, more reliable Bayesian model selections for downstream tasks. We can see that the PhyloGFN has achieved good performance regarding standard deviation on three of the datasets. However, PhyloGEN does seem to struggle with producing competitive mean marginal log-likelihodd. For VBPI-Mixture, increasing the number of mixture components results in significant improvements in terms of lower standard deviations (on all datasets) and higher mean marginal log- logelikel scores (especially apparent on the more complex datasets, e.g. DS5-8).
Tradc-off between marginal log- likeliood and standard deviations As mentioned above, in the VBPI literature, the larget metric in the marginal log-likelihood estimation task is frequently considered to be the standard deviation. Lee, the lower the standard deviation, the better, with less focus on the mean marginal likeliiel results. The reason for this might be that the mcan marginal log-likelihood scores on DS1-8 are quite saturated, which naturally instead shifts the focus to the standard deviations of the estimates. Tn Xie & Zhang (023), they alternatively motive the standard deviation criterion by pointing
11

---

Published in Transactions on Machine Learning Research (11/2024)
+7.10Be3
+2.6367e4
+3.3735e4
+1.33299e4
0.125
[
0.1
0.100
0.70
0.075
0.050
DS5
D66
D57
DSB
+8.214e3
+6.724e3
+3.733e4
+8.65e3
2.1
0.75
0.75
0.50
0.50
0.25
1,9
0.50
1.8
Mix-F
VBPI-NF
MixNF
MixNF
VEPI-NF
MixNF
VBPI-F
5m3
2=3
Figure 6: Illustrations of the mean and standard deviations of the (negative) marginal log- log-likelhhoo scores from two of the models. The plots are complementary to the scores in Table 2 where VBPI with = 3 can achieve lower mean scores, while our VBPI-Mixture with $ = consistently produces smaller or similar standard deviations.
critically demonstrating cases
the importance sampling domain, stating that the variance of the estimator good reflection of the accuracy between the importance sampling distribution and the target distribution.
Nonetheless, observe that the mean log- liglliies can be slightly improved when using Z 3 Monte Carlo samples, we additionally plot the mean log-likelihood vs. standard deviation trade-off between our S = 3 VBPI-NE-Mixture and the Z = 3 VBPI-NF. This provides the practitioner with more information when deciding which method to employ. Please see Fig.
Compute infrastructure and run times Most computations have been conducted on an AMD EPYC 7742 where two cores have been used per rum. Final runtimes are shown Table 3 The table shows the joint run time for both training and testing. Also worth noting is that multiple mixture components also multiply the number of particles, so the majority of the time increase due to the Felsenstein pruning algorithm for the likelihood model evaluation, which grows linearly.
Finally, and crucially, the code used was not optimized for run time, and so wall-clock time is not apt metric for comparisons.
Table 3: Compute time reported in minutes for fitting the model well evaluating the marginal likelihood and continuous estimate of ELBO while training every 5000 iterations using 1000 samples with a single particle. The Ufboot2 was run on 19-13900k, the command used for a single repetition for DS1 //iqtree2 -s fasta -bb 10000 -ubt -m JC69.
Data
DS1
DS2
DS3
DS4
DS5
DS6
DS7
DS8
VBPI
735.4
813.95
962.67
1087.98
1273.57
1299.55
1606.03
1660.18
Mixs=2
1984.47
2172.18
2598.17
2888.58
3322.77
3389.5
4127.52
4234.98
Mixs=3
3224.57
3425.35
4056.07
4493.87
5614.58
5604.43
6705.18
7004.28
VBPI-NE
900.32
953.8
1150.78
1278.27
1585.87
1523.02
1869.
1974.17
MixNE s=2
2385.38
2442.75
2893.15
3300.35
3800.88
3728.75
4647.45 48
4754.22
Mix
3669.8
3939.73
4576.08
5180.75
6206.52
6220.4
7568.85
7776.37
Ufboot2
0.19
0.28
0.27
0.18
0.11
0.15
0.46
0.13
5 Conclusion
We introduced VBPI-Mixtures, novel algorithm that increases the flexibility of the phylogenetic posterior approximation by utilizing recent. advances in mixtures for black-box VI. We showed that mixtures of SBNs
12

---

Published in Transactions on Machine Learning Research (11/2024)
approximate distributions that single SBN cannot, making persuasive case for VBPI-Mixtures. Experimentally, achieve state-of-the-art results on eight popular real phylogenetics datasets in terms of producing lower-variance estimators of the marginal log-likelihood, show that mixtures of SBNs can accurately capture tree-topology posteriors that non-mixture SBNs cannot, and produce KL divergences to an MCMC-based approximation of the true tree-topology posterior that decrease with S.
The drawbacks of this framework primarily lie in the increased computational time 4.2.2 compared to the baselines and the dependency on pre-computed topologies, both of which create higher barrier to entry for downstream tasks. However, are confident that the compute time can be significantly reduced with a greater focus on parallelization.
References
Yoshua Bengio, Salem Lahlou, Tristan Delou, Edward J. Hu, Mo Tiwari, and Emmanuel Bengio. Gflownet foundations. Jowrnal of Machine Learning Research, 24(210):1 55, 2023. URL http://pieirrrppppppreere/
v24/22-0364.htm html
Alexandre Bouchard-Côté, Sriram Sankararaman, and Michael Jordan. Phylogenetic inference via sequen- tial monte carlo. Systematic biology, 614)55579- 2012.
Yuri Burda, Roger Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. In International Conference on Learning Representations. PMLR, 2013.
Laurent Dinh, Jascha fock-Dic,,,ein, and Samy Bengio. Density cstimation using real nvp. arXiv preprint arXiv:1605.08888 2016.
Joseph Felsenstein. Inferring Phylogenies. Sinauer, October 2003.
R Garey, T J Near, M R Nonnemacher, and others. Molecular evidence for acanthocephala as a subtaxon of rotifera. Journal of Molecular, 1996.
Aric Hagberg, Pieter Swart, and Daniel S Chult. Exploring network structure, dynamics, and function using networkx. Technical report, Los Alamos National Lab. (LANL), Los Alamos, NM (United States), 2008.
S B Hedges, K Moberg, and others. Tetrapod phylogeny inferred from 18S and 28S ribosomal RNA sequences and review of the evidence for amniote relationships. Mol. Biol., 1990.
Daniel A Henk, Alex Weir, and Meredith Blackwell. Laboulbeniopsis termitarius, an ectoparasite of termites newly recognized as a member of the laboulbenionyyeels. Mycologia, 95(4):561155) 163661166,,,- 2003.
Sebastian Höhna, Michael J Landis, Tracy A Heath, Bastien Boussau, Nicolas Lartillot, Brian R Moore, John P Huelsenbeck, and Fredrík Ronquist. Revbayes: Bayesian phylogenetic inference using graphical models and an interactive modelsspecification language. Systematic biology, 65(4):726 736, 2016.
John P Huelscnbock and Fredrik Ronquist. Mrbayes: Bayesian inference of phylogenctic trees. Bioinformat- ies, 17(8):754 755, 2001.
Thomas H Jukes, Charles R Cantor, and Others. Evolution of protein molecules. Mommalian protein melabolism, 3:21-132, 1969.
Thomas N Kipf and Max Welling. supeviiee classification with graph convolutional networks. ICLR, 2017.
Hazal Koptagel, Oskar Kviman, Harald Melin, Negar Safinianaini, and Jens Lagergren. Vaiphy: variational inference based algorithm for phylogeny. In Advances in Neural Information Processing Systems, 2022.
Oskar Kviman, Harald Melin, Hazal Koptagel, Victor Elvira, and Jens Lagergren. Multiple importance sampling elbo and deep ensembles of variational approximations. In International Conference on Artificial Intelligence and Statistics, pp. 10687-10702. PMLR, 2022.
13

---

Published in Transactions on Machine Learning Research (11/2024)
Oskar Kviman, Ricky Molén, Alexandra Hotti, Semih Kurt, Victor Elvira, and Jens Lagergren. Cooper- ation in the latent space: The benefits of adding mixture components in variational autoencoders. In International Conference on Machine Learning, pp. 18008 18022. PMLR, 2023.
Clemens Lakner, Paul van der Mark, John P Huelsenbeck, Bret Larget, and Fredrik Ronquist. Efficiency of markov chain monte carlo tree proposals in baycsian phylogenctics. Syst. Biol., 57(1):86-103, February 2008.
Bui Quang Minh, Minh Anh Thi Nguyen, and Arndt von Haeseler. Ultrafast approximation for phylogenetic bootstrap. Mol. Biol. Evol, 30(5):1188-11 May 2013.
Andriy Mnih and Danilo Rezende. Variational inference for monte carlo objectives. In Maria Florina Balcan and Kilian Q Weinberger (eds.), Proceedings of The 33rd International Conference on Machine Learning, volume 48 Proceedings of Machine Learning Research, pp. 2188-2196 New York, New York, USA, 2016. PMLR.
Antonio Khalil Moretti, Liyi Zhang, Christian A Naesseth, Hadiah Venner, David Blei, and Itsik Pe'er. Vari- ational combinatorial sequential monte carlo methods for bayesian phylogenetic inference. In Uncertainty Artificial Intelligence, pp. 971 071-981. PMLR, 2021.
Warren Morningstar, Sharad Vikram, Cusuh Ham, Andrew Gallagher, and Joshua Dillon. Automatic dif- ferentiation variational inference with mixtures. In International Conference on Artificial Intelligence and Statistics, pp. 2250-3258. PMLR, 2021.
Eric Nalisnick, Lars Hertel, and Padhraic Smyth. Approximate inference for deep latent gaussian mixtures. In NIPS Workshop Bayesian Deep Learning, volume 2, pp. 131, 2016.
Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. In Artificial intelligence and statistics, pp. 814- 822.. PMLR, 2014.
Danilo Rezende and Shakir Mohamed. Variational inference with normalizing flows. In International con- ference on machine learning, pp. 1530 1538. PMLR, 2015.
Geoffrey Roeder, Yuhuai Wu, and David K Duvenaud. Sticking the landing: Simple, lower-variance gradient estimators for variational inference. Advances in Neural Information Processing Systems, 30, 2017.
Amy Y Rossman, John M McKemy, Rebecca A Pardo-Schultheiss: and Hans-Josef Schroers. Molecular studies of the bioncctriaccac using large subunit rDNA sequences. Mycologia, 93(1):100-110, January 2001.
Francisco JR Ruiz, Michalis Titsias, and David M Blei. Overdispersed black-box variational inference. Uncertainty in Artificial Intelligence, 2016.
Yuge Shi, Brooks Paige, Philip Torr, et al. Variational mixturre-f-exxpets autoencoders for multi-modal deep generative models. Advances in neural information processing systems, 32, 2019.
Liangliang Wang, Alexandre Bouchard-Cott and Arnaud Doucet. Bayesian phylogenetic inference using combinatorial sequential monte carlo method. Journal of the American Statistical Association, 110(512): 1362-1374, 2015.
Liangliang Wang, Shijia Wang, and Alexandre Bouchard-Côté. An annealed sequential monte carlo method for bayesian phylogenetics. Systematic biology, 69(1:1555-183, 2020.
Chris Whidden and Frederick Matsen IV. Quantifying meme exploration of phylogenetic tree space. Systematic biology, 64(3):472-45) 491, 2015.
Tianyu Xie and Cheng Zhang. Artree: A deep autoregressive model for phylogenetic inference. arxie, 2023.
Wangang Xie, Paul O Lewis, Yu Fan, Lynn Kuo, and Ming-Hui Chen. Improving marginal likelihood estimation for bayesian phylogenetic model selection. Systemalic biology, 60(2):150-160, 2011.
14

---

Published in Transactions on Machine Learning Research (11/2024)
Z Yang and A D Yoder. Comparison of likelihood and bayesian methods for estimating divergence times using multiple gene loci and calibration points, with application to radiation of cute-looking mouse lemur species. Syst. Biol., 2003.
Anne D Yoder and Ziheng Yang. Divergence dates for malagasy lemurs estimated from multiple gene loci: geological and evolutionary context. Mol. Ecol., 13(4):757-77) April 2004.
Cheng Zhang, Improved variational bayesian phylogenetic inference with normalizing flows. Advonces in neural information processing systems, 33:18760-18771, 2020.
Cheng Zhang. Learnable topological features for phylogenetic inference via graph neural networks. ICLR, 2023.
Cheng Zhang and Frederick A Matsen IV. Generalizing tree probability estimation via bayesian networks. Advances in neural information processing systems, 31, 2018.
Cheng Zhang and Frederick A Matsen IV. Variational bayesian phylogenctic inference. In ICLR, 2019.
Cheng Zhang and Frederick A Matsen IV. A variational approach to bayesian phylogemetic inference. Jownnal of Machine Learning Research, 25 115):115556 2024. URL http:///paar..or//pppprrss/////
N Zhang and M Blackwell. Molecular phylogeny of dogwood anthracnose fungus (discula dostructiva) and the diaporthales. Mycologia, 2001.
Ming Yang Zhou, Zichao Yan, Elliot Layne, Nikolay Malkin, Dinghuai Zhang, Moksh Jain, Mathieu Blanchette, and Yoshua Bengio. PhyloGFN: Phylogenetic inference with generative flow networks. The Twelfth International Conference on Learning Representations, 2024.
A A Brief Introduction to Bayesian Phylogenetic Inference for Machine Learning Researchers
This introduction aims to briefly explain the basic concepts required to understand the generative model provided in the main text.
Phylogenetic trees capture evolutionary relationships among species and provides valuable insights into life's evolutionary history. Within this domain, phylogenics are often depicted as bifurcating tree graphs, where nodes represent common ancestors, and branches (edges) signífy evolutionary events and genetic distances between species. This framework enables an understanding of species relatedness, ancestry, and the evolutionary processes governing life's diversity.
Bayesian phylogenetic inference builds upon this framework by applying Bayesian statistical methods to infer the evolutionary history. It allows for probabilistic approach to model uncertainty and variation, considering prior beliefs about evolutionary parameters and updating these beliefs as new data is incorporated. It is common to use DNA or protein sequences as data since it describes different attributes of the species, and the edges represent mutation between species. Through sampling from a posterior distribution and utilizing tools like Markov Chain Monte Carlo (MCMC) methods or variational inference, Baycsian phylogencties offers robust and nuanced view of evolutionary relationships, integrating multiple sources of information and providing rigorous statistical foundation for evolutionary hypotheses.
Mainly, two latent variables are regarded as important in Bayesian phylogenetic inference. First, the tree topology, T, a binary tree with the observations assigned to its leaves. The tree-topology space grows as (2n - 33!! for rooted and (2n - 5)!! for unrooted trees, where n is the number of leaves (observations/tar /taxa). Furthermore, each odge, e, of the topology associated with a positive continuous variable, the branch lengths, b(e). The Cartesian product of discrete and continuous spaces makes inference in phylogenetics a challenging task.
15

---

Published in Transactions on Machine Learning Research (11/2024)
B Conflicting Tree Topologies
Here we construct a realistic scenario where DNA sequences induce conflicting tree-topologies in the posterior that cannot be modelled by the vanilla SBN. On the other hand, they can be captured by VBPI-Mixtures.
It is well-known that DNA data sometimes has confficting signals. Here construct toy example to demonstrate how T1 and 72 have higher posterior support than тз and T1. For simplicity, use the connection between lower parsimony score and higher likelihood when branch lengths are short (to appear in the appendix). First, consider that the leaves in Fig. 2 have DNA sequences with nucleotides at sites and j specified in the table below (the first and second columns represent and j, respectively).
Table 4: Nucleotide assignments to sites and j in the observations.
Sites (1)
(
3
(4)
(5)
(6)
A
G
j
A
G
A
We compute the parsimony scores for each clade (A, A', B, B') and the cost of joining two clades to form
TI:-., 74- The cost of transitioning from one nucleotide (A, C, G, T) to another is 1.
In the first table below, we start by calculating all possible parsimony scores for each clade. Here X-Y denotes the parsimony score, Y, if the nucleotide site i or j in the ancestral sequence in the root is X. Bolded is the best (lowest) score.
Table 5: Parsimony scores for the four different subtrees
MP
A= ((1,))3)
B (=,
A' (1.(23))
B' =4,
6th
A-2,C-116--8T---
A-2,C-3.--11-3 A-11C--8G--1,--
A-11C-1,6-2T-2 A-2--1,6-88---
A-11.--2.-111--- A-2,C-3,6-1,--33
j:th
A-11,C-1,6-2T--2
Table 6: Parsimony scores for all possible tree topologies
MP
7AA
A(A'A
T(AA
Ta (A' AB)
6th
A 1.CC--3,C-.TT--
A A -3
A A-3.C-3.C 44
3-3,C-3.CC-33T- A 3-3.--3. 14
jith
A 2 33 33T 14
A 3,,C-33,6-3.TT11
A-3,C-33G >4
Now we can see that Ti and give us a belter parsimony score, demonstrating that the scenario exemplified in Fig. 2 can occur in biological data when there are conflicting signals.
We direct readers who want learn more about parsimony scores to Chapter in Felsenstein (2033)
C
Gradient Derivation
Here show the full derivations of the gradients w.r.t. d.. For completeness, recall that
PBB T X)
[x, Be, 7b)
(12)
3 E 1=1 av, (B )ot)) (T+)
and i = log K Ek f(x, Bhg T.), where Bb, - are simulated from qv..d, T).
That is, we are interested in the gradient of Eq. w.r.t. the SBN parameters for one the mixture components, say i,
V oe L(X;K,) = V..
Σ E
yvs (B,7) [
(13)
5
s=1
16

---

Published in Transactions on Machine Learning Research (11/2024)
There are two cases to take into account in the sum, either i or i + s. Starting with i and using the product rule,
V 0i s Bgg.K k " (,,,,,
[,,4
V
Σ T.i
1:K 2o (T! (T.
B -n, (B) [1,5
(14)
S
Egg K m0 (8)771xg
V Vo. (T( )
(15)
8
ne (B) V vo 1: + go. (Ti 1:K VV
k
(16)
8
Recalling the identity that Vego(z) = 36 V @ log (o we start by rewriting the first term inside the expectation in Eq.
16 as
i Vo Io. (0. = (T! 1: ) Vo. log qo: (T!
(17)
k
q6. (T. (T:
1: L
V d. go, (T,),
(18)
R=1
and then the second term
K
Vo. 4o. (T. (T..:K)) VV. P.K (T! I:KK ) Vo.
log
1
Σ fx, ,
(19)
K
m V d: log (x, B5,
(20)
K
1
1:K -0. (T,
log
qv, (B. (B: De,, (T²")
(21)
s
k=1
-=1
where /.
flx.88 x!)
Exchanging the two terms Eq.
(16 with Eq. (18 and 211
respectively,
Σ x- 18 "
we get
K
E Vo. s
10, (B,7)
[
E
0, (B,7)
i
Vo. log qdh (T+*
(22)
s
k=
E go (B,r)
w
V d, log
1
Muy (B)))=))ss (+:)
(23)
5
For # s we may move the gradient operator into the expectation directly and reuse the derivation of Eq.
21
E 8
yos vs (B,7)
[ik]
=
E Boris, vs (8,7)
Vo, L.
(24)
K
1
E
Eggn. (B,r)
tht V de log
1
Σ
Iw (B% (++)]
S
5
h=1
==1
Considering both cases, we return to Eq.
13
V oe L(X;
=
Buorio (B..r)
E
i
V de log Zo. (T+)
S
k=
(25)
EE Eg (B,r)
Σ
2 v
log
1 Σ
(B. (3,))))0,
s
17

---

Published in Transactions on Machine Learning Research (11/2024)
which is the expression for the gradient we need in order to apply the VIMCO estimator (see Sec. 3.1.2
D Unbiasedness of the VIMCO Estimator Applied to Uniform Mixtures
Recall that the VIMCO estimator applied to uniformly weighted mixture of variational distributions is given in Eq. 110 We want to analyze the expectation of the R.H.S. of this equation. Note also that 01, Us are disjoint sets of parameters, and so they are updated via the gradients independently.
We now expand Eq.
100 using Eq. and apply expectations:
K
E Rer.o (B,r)
L
K V
log (o!T!)
(26)
K
E Socies (B,r)
Σ
log
Σ
(1, , )8 + (Z, ,,,,))))))))
Vo log qo, (T+)
(27)
5
K
k=1
9
s
E
dy ny (B,r)
w Vo. d.
log
qis (B= (T-1)
(28)
5=1
Now, using the REINFORCE identity, we can rewrite the first and the third term above (for more details, you can go backwards in our gradient derivations from Eq.
22 show that:
K
K
8 (+,)
E B,nn 5
(8,r)
iik
de log lo (T,))
E Re (8,r)
Σ
th
Vo.
log
Σ ==1
qos (B,
(29)
S
k=1
Vo. 8 Le (B.r)
i
= V ((X;
(30)
From this equality, and by subtracting the second term in Eq. 26 from Vo, ((X;KK, get
1 V ((X: - B (B.r)
Σ log
1 K
Σ byyk
E
B (2, V log tog (T.7)
(31)
where the second term is the control variate. The control variate is, however, independent of the mixture formulation (except for the 1/S scaling which does not affect the unbiasedness of the term). As such, since VTMCO is an unbiased estimator of the gradient (Mnih & Rezende 2016), then so our version of it.
E
Additional Experimental Results and Implementation Details
E.1
The Two-Level Hierarchical Model
Rumning a grid search for all five algorithms (S 1,, using n, 5 and n2 10, we found that the learning rates 0.01, 0.1, 0.1, 0.2, and 0.25, respectively, were optimal. That is, these achieved the smallest KL divergences. For larger learning rates, S did not converge converge to worse KL divergences. The optimal learning rates found in the grid searches were used in all subsequent experiments.
In Fig. we visualize the KL curves as functions of the number of training iterations. For all configurations of K, π1 and 72, the S model performs best. The pattern- S converges slower and to worse KL divergences than 2ohoss also when all models use the same number of importance samples, K (shown in Fig
E.2 Visualization Details and More Plots
We employed the software as in Whidden & Matsen IV (0015), iec rSPR was used to determine the distances between the topologies. Additionally, adopted the methodology for cluster creation. This involved
18

---

Published in Transactions on Machine Learning Research (11/2024)
1=5, 2-10
1-3-00,n 2-20
14
S-1 100
s-1 (K-20)
5-2 10-100
5-2 (K-20)
Smd 100
142
Se3 (K+b)
5m4 12-
S-4 (K-5)
1.-
3-5 X--
3.0
5-3 (S-1)
0.8
0.-
2.4
2300
5000
1500
20000
12:00
:3000
17500
20000
0
2500
5000
7500
10000
12500
25000
17300
20000
(b)
1-.100
1_1_1000,
2000000
16
14
24-
IK-5)
[cm-4]
12
1U
08-
04
02-
2500
5ooe
1500
20000
12500
15000
37505
20300
0.
2506
5000
7500
10350
12500
15070
17:00
20000
(c)
(d)
Figure 7: KL curves for the two-level hierarchical model using different configurations of K, ni and n2.
assigning the most probable peak to cluster and subsequently assigning all unassigned trees to the same cluster if their distance from the peak tree within one standard deviation below the mean distance of all unassigned trees. This iterative process continued until all topologies were assigned, or until eight clusters were reached.
For graph creation, we employ the Graphviz layout known as Scalable Force-Directed Placement (SFDP), in conjunction with the NetworkX library Hagberg et al. (2008). The clusters are represented by different colors, and the size of each node is determined by the normalized sampling frequency Fig g Morcover, edges between the topologies are only displayed if their distance is exactly one. It's worth noting that the iSPR distance measure is utilized, which counts the number of changes similar to the approach used in the MCMC method.
To ensure that the visualization focused on the most credible information, we imposed a constraint by limiting the nodes to the 95% most credible set. This ensured that only the most reliable nodes were included. Additionally, to manage computational resources effectively, we set a maximum limit of 4096 nodes for the graphs.
A similar approach was used for Fig. with the main difference being that we sampled from the components and displayed the joint set of topologies. The colors in this figure were based which component sampled the topology the most. This representation was considered an approximation of the posterior.
F Limitations
We use SBNs and branch length models form our mixture approximations. This introduces a larger number of model parameters. In our current implementation, the training time prolonged did
19

---

Published in Transactions on Machine Learning Research (11/2024)
(a)
(b)
(c)
(d)
(e)
Figure 8: Visualization of samples from the tree-topology posterior using a 1,000,000,000 iterations long
MCMC run on (a) DS1, (b) DS2, (c) DS3, (d) DS5 and (e) DSG. Nodes represent unique tree-topologies and are colored based on cluster assignments, illustrating the multimodality of the treo-topology postcrior. More details Sec. 2
not parallelize the parameter updates of the parameters of the mixture components. This can, however, be done, in order to heavily decrease training time.
Additionally, using shared parameters for the mixture components can also be utilized, if the practitioner is running on limited memory budget. However, shared weights would violate the assumption of disjoint SBN parameter sets in Appendix D Nonetheless, devising clever modifications to reduce the number of parameters and training limes for mixtures in black-box VT is an exciting future research field.
G Broader Impact
Bayesian phylogenetic inference algorithms are crucial for researchers to reason about uncertainty in their evolutionary findings. Variational inference algorithms provide a compelling alternative to MCMC-based algorithms as parametric approximation obtained. This implies that VI, and VBPI specifically, can be used in settings where the application of MCMC is less straight for instance in out-of-distribution detection, evaluation on held-out data. Also, as have shown in our experiments in this paper, model evaluation can be more robust when using VI over MCMC, resulting in smaller variance of the estimator of the marginal log-likelinn likelihood. This is an important feature for downstream tasks.
20

---

Published in Transactions on Machine Learning Research (11/2024)
(a) DS1
(b) DS2
(c) DS3
(d) DS4
(e) DS7
Figure 9: Visualization of uniformly weighted S 2-component mixture of SBNs on (a) DS1, (b) DS2 and (c) DS3, (d) DS4 and (e) DS7, where each node corresponds to unique tree-topology The upper shows the distribution of five million sampled tree topologies from the first component, where node, is colored blue if 40, (T) > 4og (T), or orange otherwise. Vice versa for the lower row. The size of node is determined by its sampling frequency, which is why nodes with low frequency appear black. The components clearly spread out, exploring different parts of the space.
21


Reviews:
[{"review_id": "iGNTo9j8rn", "reviewer": "Reviewer_9pfv", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "I have reviewed the previous submission of this paper, and I’m pleased to see that most of my earlier comments have been addressed, including the inclusion of running times, clarification on the use of pre-computed topologies, and softening of the language. Given these improvements, I would support acceptance.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors propose a mixture model of tree topologies to enhance parameter learning in Variational Bayesian Phylogenetic Inference (VBPI). Additionally, they derive a VIMCO estimator for mixtures.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "I would suggest adjusting the position of the legend in Fig. 3 or removing the bounding box for better clarity.", "broader_impact_concerns": "N/A", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "CsMeUBBe7y", "reviewer": "Reviewer_5pL4", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The authors have addressed the concerns raised in the previous review round effectively in this revised version of the paper. The modifications have significantly improved the clarity and addressed the key issues I highlighted earlier. Based on these improvements, I believe the paper is now in good shape for acceptance.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper proposed a variant of variational Bayesian phylogenetic inference that uses mixture of SBNs instead of a single SBNs for tree topology posterior approximation. The authors show that certain distribution of tree topologies can not be properly learned via a single SBN. Stochastic gradient estimators (following the VIMCO method) have been derived. Experiments on a toy hierarchical model and real data Phylogenetic posterior estimation tasks demonstrate the effectiveness of the proposed method.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "No requested changes.", "broader_impact_concerns": "No broader impact concerns are needed so far.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "KKpUFjnlbf", "reviewer": "Reviewer_LrA3", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Strengths\n\n* The paper is written clearly and the derivations seem correct.\n* The case for using a mixture variational family is well-motivated via Example 3.1.\n* The method is simple to implement and the authors plan to share the code to reproduce it.\n* The evaluation on the toy model illustrates the advantages of mixtures well, thus supporting Example 3.1. \n\n## Weaknesses\n\nThe authors have addressed the major weaknesses from the previous submission. Below are some less critical points that would help to improve the paper.\n\n### Additional metrics\n\nIn Table 1 and Figure 3, the authors use forward KL divergence, $\\text{KL}(p(\\tau \\mid X) \\mid\\mid q\\_{\\phi}(\\tau))$, to measure the accuracy of the variational approximation. (Figure 3 doesn't specify which KL is used but it seems to be the forward KL in the provided code---this should also be clarified in the caption.) But, the objectives optimise an _importance-weighted_ version of a reverse KL. This means, that the variational distribution should be used as importance sampling proposal distribution, e.g. by sampling from $q\\_{\\phi}(\\tau)$ and then resampling using the importance weights $p(\\tau \\mid X) / q\\_{\\phi}(\\tau)$, instead of treating $q\\_{\\phi}(\\tau)$ as a direct approximation of the posterior $p(\\tau \\mid X)$ [1]. As such, it would be useful if the accuracy of the variational distribution would be assessed using some importance sampling metrics, such as the variance of the importance weights which should be low if $q$ is a good proposal.\n\n* In Table 1, it would also be good to include the results for $\\text{VBPI}\\_{Z=3}$, which seem to be very close to $\\text{Mix}\\_{\\text{NF},S=3}$ in terms of MLL, and thus would be good to know their posterior accuracy.\n\nReferences:\n\n[1] Cremer et al (2017). Reinterpreting Importance-Weighted Autoencoders.\n\n### Clarity\n\n* In the introduction, the authors state that \"Using a toy experiment [...], we first show that the mixture components specialize in different parts of the solution space\". However, the re-worked experiment section 4.1 does not directly show this, so the claim must be adjusted. \n* In Section 3, first sentence of the \"Mixtures promote exploration\" paragraph states: \"Note that minimizing the denominator in Eq. (6) corresponds to diversifying the mixture distribution\". However, due to the log-sum in the objective, minimizing _just_ the denominator does not make much sense (because it doesn't factor into an entropy term like in the standard ELBO). Hence, the statement needs to be clarified.\n* Sections 3.1.1 and 3.1.2 implicitly assume that the parameters of the mixture components are disjoint, as explained in Appendix D. This assumption should also be highlighted in the main text. \n\t* Moreover, In Appendix F, the authors suggest that using shared parameters for the mixture components could be used to reduce memory budget. But, this would violate the assumptions before and should also be highlighted.\n* The observation in Section 3.1.1 \"(ii) as the second term is negated, maximizing it corresponds to diversifying the mixture distribution w.r.t. $\\phi_i$,\" is unclear. It is not clear what it means to \"diversity the mixture w.r.t. $\\phi_i$\", or why exactly does maximizing the term correspond to any kind of diversification. If there are prior work on this, then a precise reference might suffice. \n* The second paragraph in Section 4 claims that \"the increased model flexibility and promotion of exploration translates\ninto better marginal log-likelihood estimates\". But the results mostly rely on improvements to the standard deviation of the MLL, so this claim should be adjusted.\n\n### Minor typos:\n\n* Page 2, first sentence of section 2. $X$ is not yet introduced, hence confusing.\n* Page 2, after eq 2, \"reverses\"->\"reverse\"\n* Section 4.1, last paragraph, second sentence. Two typos. \"This [is] to equate ... samples use[d] during learning.\"\n* Section 4.1. last paragraph \"Figures 3\" -> \"Figure 3\"\n* Figure 7: (b) is the same as (a). Should probably be D2?\n* Appendix D, third sentence \"Not\"->\"Note\"", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_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": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper describes a mixture approach for variational Bayesian phylogenetics, where the goal is to perform probabilistic inference of the posterior over evolutionary topologies. Using a simple example, the authors motivate that  existing work that uses a single subsplit Bayesian network (SBN) may not accurately approximate this posterior. To resolve this, the authors propose using mixtures of SBNs, which increases the flexibility of the variational posterior and may be able to better capture the true posterior. The authors also extend existing analytical results on discrete variational inference with importance-weighted bounds (VIMCO) to variational mixtures.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Please see the weaknesses above.", "broader_impact_concerns": "None of note.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper is a resubmission to TMLR. The three reviewers found that the submission addressed all the raised concerns, and therefore the paper is ready for acceptance. See details of the previous submission [here](https://openreview.net/forum?id=yuhG3VHK5f).

The paper is well written, and its technical correctness has been significantly improved with respect to the previous version.
Decision: Accept (poster)