VaiPhy: a Variational Inference Based Algorithm for Phylogeny
Authors: Hazal Koptagel
Oskar Kviman
Harald Melin
Negar Safinianaini
Jens Lagergren
Venue: NeurIPS 2022 Accept (2022)
Keywords: Phylogenetic Inference
Variational Inference
Bayesian Inference
Combinatorial Sequential Monte Carlo

Abstract: Phylogenetics is a classical methodology in computational biology that today has become highly relevant for medical investigation of single-cell data, e.g., in the context of development of cancer.  The exponential size of the tree space is unfortunately a formidable obstacle for current Bayesian phylogenetic inference using Markov chain Monte Carlo based methods since these rely on local operations. And although more recent variational inference (VI) based methods offer speed improvements, they rely on expensive auto-differentiation operations for learning the variational parameters. We propose VaiPhy, a remarkably fast VI based algorithm for approximate posterior inference in an \textit{augmented tree space}. VaiPhy produces marginal log-likelihood estimates on par with the state-of-the-art methods on real data, and is considerably faster since it does not require auto-differentiation. Instead, VaiPhy combines coordinate ascent update equations with two novel sampling schemes: (i) \textit{SLANTIS}, a proposal distribution for tree topologies in the augmented tree space, and (ii) the \textit{JC sampler}, the, to the best of our knowledge, first ever scheme for sampling branch lengths directly from the popular Jukes-Cantor model. We compare VaiPhy in terms of density estimation and runtime. Additionally, we evaluate the reproducibility of the baselines. We provide our code on GitHub: \url{https://github.com/Lagergren-Lab/VaiPhy}.

VaiPhy: a Variational Inference Based Algorithm for Phylogeny
Hazal Koptagel Oskar Kviman 12+ Harald Melin 12
Negar Safinianaini Jens Lagergren 12
School of EECS, KTH Royal Institute of Technology, Stockholm, Sweden
2Science for Life Laboratory, Solna, Sweden
"Equal contribution
(koptagel, okviman, okynnan aaraldme, negars, jengals (atsh.@kkhh...
Abstract
Phylogenetics is classical methodology in computational biology that today has become highly relevant for medical investigation of single-cell data, e.g.. in the con- text of cancer development. The exponential size of the tree space is, unfortunately, a substantial obstacle for Bayesian phylogenetic inference using Markov chain Monte Carlo based methods since these rely on local operations. And although more recent variational inference (VI) based methods offer speed improvements, they rely on expensive nuto-differentiation operations for learning the variational parameters. We propose VaiPhy, a remarkably fast VI based algorithm for approxi- mate posterior inference in augmented tree space. VaiPhy produces marginal log-likedn
liglikkei estimates on par with the stat-oofthhe-as methods on real data and is considerably faster since it does require auto-differentiation. Instead, VaiPhy combines coordinate ascent update equations with two novel sampling schemes: (i) SLANTIS, a proposal distribution for tree topologies in the augmented tree space, and (ii) the JC sampler, to the best of our knowledge, the first-ever scheme for sampling branch lengths directly from the popular Jukes-Cantor model. We compare VaiPhy in terms of density estimation and runtime. Additionally, we evaluate the reproducibility of the baselines. We provide our code on GitHub: attps://gil //g:t/bbttttt.t.// www Loom/Lagrrgee---- Latrrrrren--aaaaaPhyy
Introduction
Phylogenetic software has for long time been applied in biological research, and biological findings relying on phylogenetic trees are frequent. Moreover, due to the emergence of single-cell sequencing. phylogenetic inference for bifurcating and multifurcating trees can now be utilized in medical studies of how cells differentiate during development and how tumor cells progress in cancer.
When attempting to infer the posterior distribution over phylogenetic trees for [X] taxa, represented as sequences, the main obstacle is the exponential size of the tree space, i.e., the of candidate trees having those taxa as leaves. In Markov chain Monte Carlo (MCMC) based Bayesian phylogenetic inference, e.g., MrBayes [28] and RevBayes [20], this tree space is explored by performing local operations on the trees. In order to avoid local operations, Bouchard-Cône et al. [3] launched Sequential Monte Carlo (SMC) as an alternative methodology for Bayesian phylogenetic inference.
Probabilistic phylogenetic inference is machine learning problem, and the key to unlocking this erucial application be expected to be found in the machine learning methodology: in particular, considering the capacity of the VI methodology, [32, 1]. to deliver impressive performance gains for Bayesian inference, e.g., compared MCMC. Expectedly, such gains have recently been made within the realm of phylogenetic analysis.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

4; - =;
4: - 100
200
4 300
500
0.05
Figure Visualization of branch lengths, b, sampled directly from the Jukes-Cantor model. The figure shows samples for different numbers of estimated mutations along an edge, Di.j. when the total number of sites are M = 2000.
For instance, inference based on the computationally costly CAT model, which allows different equilibrium frequencies across sequence sites, has been made substantially more efficient by the application of VI, [5]. Also, [11] applied VI to render the analysis of sequence evolution model more efficient. They consider the efficiency improvements obtained by VI, in particular using the generic STAN software [11], when analyzing single fixed tree topology. Interestingly, the VI approach can also be used to analyze a large number candidate phylogenetic trees. In [38], so-called subsplit Bayesian networks (SBNs) were introduced in order to represent the posterior over given set of candidate trees. This approach has subsequently been developed further by introducing an edge (a.k.a. branch) length distribution [39] and making the edge length distribution more expressive by taking advantage of normalizing flows [37]. In particular, [39] obtains posterior distributions that improve on the marginal likelihood obtained by using MrBayes, [21], in combination with the improved Stepping Stone (SS) method, [8].
The SBN approach evades the exponentially sized tree space by considering a given set of candidate trees and the partitions the given taxa defined by those candidate trees. noteworthy that the VI methodology carlier has been successfully applied to bipartite matching and a few other combinatorial problems [2]; the phylogeny problem, however, involves both combinatorial and continuous parameters.
Multifurcating trees are common in biology, and biologists are well familiar with them. For instance, in clone tree that represents the evolutionary relationships between the subclones of tumor, nodes may have more than two children. Even in cases with no biological support for multifureating trees, biologists are used to working with multifurcating trees. E.g., the consensus method used to analyze the posterior of a MrBayes run may produce a multifurcating tree.
Contributions: We present novel framework for phylogenetic inference that produces marginal
likelihood estimates on par with the oftte----hhearr VI methods while being up to 25 times faster than VBPI-NF on real datasets. The framework is on the top level a mean-field VI algorithm and is comprised of the following three novel methods:
SLANTIS is a proposal distribution that proposes multifurcating spanning trees. Sampling from SLANTIS is performed following an approach suggested by [6] for sampling a perfect matching in bipartite graphs: is fast, and the sampled tree's compensation weight is easily computed without additional cost of time.
The JC sampler is novel scheme for sampling branch lengths directly from the Jukes- Cantor (JC69) model [22]. This a significant contribution, as no previous sampling scheme is known for this task and the JC model arguably is commonly used substitution model in phylogenetics.
VaiPhy is a VI algorithm for obtaining posterior approximations over the set of branch lengths, the ancestral sequences, and tree topologies, given the leaf set observations. VaiPhy utilizes SLANTIS and the JC sampler order to compute importance weighted evidence lower bound (IWELBO; [4]) and its update equations. Compared to existing VI methods in phylogenetics, VaiPhy extends the target distribution by considering an auxiliary variable, the ancestral sequences, and it does not require any input consisting of precomputed trees.
We demonstrate the flexibility of our framework by using its learned parameters to parameterize two new proposal distributions for the Combinatorial SMC (CSMC) algorithm. We compare CSMC

---

method, 6-CSMC, with the other two main SMC based methods, outperforming them in terms of marginal likelihood estimation on real data.
Finally, we contribute to the field of Bayesian phylogenetic inference by benchmarking the most prominent methods in the field seven popular real datasets.
Background
In the Bayesian phylogenetic setting, interested in the posterior distribution over phylogenies pairs of tree topologies and sets of branch lengths - given observations on the leaves,
Pa (T, po( B))
(1)
(x
where r, BB and X denote the tree topology, the set branch lengths and the observations, respectively. More specifically, B contains a single branch length per edge, e, in the set of edges т, i.c., B = (b(e): E E. Unfortunately, pe (T, B| is intractable, so we need to approximate it using Bayesian inference. Before going into how approximate this quantity, we first provide some crucial notation.
In phylogenetics, one is often concerned with bifurcating X-trees. These trees have labeled leaf nodes, and their internal vertices have either two (root a rooted tree) or three degrees. Meanwhile, multifurcating trees' internal vertices can have higher degrees. We use M, E, and N to denote the sequence length, the set of nucleotides, and the number of vertices, respectively. Further, we let |X] be the number of given taxa, V and F. be the set of vertices and edges, while (V, F.) is acyclic graph.
Let Z be the probabilistic ancestral sequence of node i position m, and denote the model parameters, such as the branch length prior parameters or the parameter choices in the substitution model. Then, the evolutionary model along edge, (i,j), is po (Z a: 2Z B) exp(((ee Qu where Q is the rate matrix, parameterized according to the JC69 model, and o,8 B E Σ are nucleotides. The likelihood is defined as
(x(X = Σ 21,
ZEEM
and can be computed effectively using Felsenstein's pruning algorithm [9].
Proposed Framework
Existing Bayesian phylogenetic approaches often perform inference in a sample space over bifurcating X-trees and branch length sets, S2. In contrast, we propose to work in an augmented sample space, A. Firstly, we work in a tree topology space for multifureating trees with labeled internal vertices. Secondly, we extend the target distribution in Eq. (1) by considering ancestral sequences on the internal vertices as an auxiliary variable, Z.
Apart from that, our augmented space offers higher modeling flexibility (clearly, since Ω C A); it enables posterior estimation over Z, which further lets us estimate the number of mutations between all vertices in the sample space. This is key insight that ultimately allows us to devise VaiPhy, mean-field coordinate ascent VI [19] algorithm.
3.1 VaiPhy
In the augmented space, A, we wish to approximate the posterior distribution (T( B, Z(X) x Po( (X, r),
using a factorized variational distribution
(3)
3

---

where [( = II I(A)) (Z)()). Here, [A) denotes the set of internal vertices in the augmented space. We will use ( to compute dij [0, M), the expected number of mutations the edge j), for all possible edges. Notice that the maximum number of observed mutations between two vertices is the number of sites, M.
As standard in classical VI, we alternatively update the distributions Eq. (3) order to maximize the evidence lower bound (ELBO),
log (0X
C E (T,,8
v(B(X))(T)X))((/x)
over iterations. Each iteration, we also update Φ, following [12]
Øis = F
(5)
083
sj where Fa.33 the estimated number of times nucleotide mutates from to 3, defined as Fig
M ((ZT" xxX"" "a (z ")
The update equations for the variational distributions are derived in a mean-field fashion. Starting with X)) we find that
log (Z,lx) = 1B (X,
(6)
Σ
b(,,)))
where = denotes equal up to an additive constant, YN. (1) is the set of sequences of node i's neighbors in 7, and -i is the set of latent nodes except node i. Continuing with the two other approximations, we have
log 4" (T|X) $ By(B.2(x) [log (x, (,,,,
(7)
$ (y(Bxx)
Σ
a(Zr, 2,1 221
(U)) 21.2
and
(b(i,)))) 1 Eg (011,
[log (x,,
(8)
Eg(r X)
Σ (02,,7711))0gg900222,, Zi,Zs
+ log (b(i,
The exact computation of the expectations in Eq. (6) and (8) involves summing over all possible trees, which computationally intractable. Luckily, we can easily approximate the expectation using importance sampling, and our proposal distribution, (T),, is outlined in the subsequent Sec. 3.2. We direct the reader to Appendix for the derivations of the above variational distributions and the details of the algorithm.
To obtain good starting points for the variational distributions, we initialize with trees obtained using the Neighbor-Joining algorithm [14]. The details of the initialization are provided in Appendix A.2. Alg. a high-level algorithmic description of VaiPhy. The proof regarding the natural gradient presented in Appendix I.

---

Algorithm 1 VaiPhy in pseudocode - i is the total number of iterations and v) is the natural gradient step size
1: Input: Ø, 77,
2: for iter = :....: do
3: Approximate log X) using 5p(T) (Eq. (6))
4: ((2.) + (1 7) Zi) nd" (n(().)
Pij 6 -o P [Eq. (5))
5:
6:
((T + (X) {Eq. (7)
7:
((3lX) 9" XB using so (T) (Eq. (8)
8: end for
9: return qTT B, x
3.2 SLANTIS
We propose the method Sequential Look-Ahead Non-local Tree Importance Sampling (SLANTIS), which enables sampling of a tree, т, from distribution, (T( SLANTIS has support over multifur- cating trees with labeled internal vertices and is parameterized by Φ. To describe the algorithm, we next introduce some notation.
Using Φ, SLANTIS immediately constructs W, an N x N matrix describing all pairwise edge weights based on the substitution model,
Wi
Σ Sy (b(i.)) Qas) + 0,Bay
Σ (M (1y)) (b(i,)) Qae), 0,ftomm
(9)
where Q is the rate matrix of the JC69 model and a and are nucleotides. Moreover, let be graph with vertex set {N} with edge weights specified by W. Let Gr GAA where A is the (V, E) of the leaves, i.e., the graph induced by the internal vertices of our trees. The pseudocode of the algorithm is presented in Appendix B. The main idea, which is borrowed from [6]. is to order the edges of the graph and for each edge in the order, make a choice whether to include in the sampled tree or not. The choice is made based on Bernoulli trial, where the two probabilities are proportional to the weight of the maximum spanning tree (MST) containing the already selected edges and the currently considered edge and that of the MST containing the already selected edges but not the currently considered edge.' The action associated with the former outcome is to include the edge in the sampled tree, and the other outcome is not to include it.
To facilitate fast computation of the required MSTs, first consider the edges of MST T1 of G. After the edges of Ti, the edges of an MST Tr of G Ti can be found. In general, T, is an MST of G -== and after the edges of - the edges of T, can be found. In the algorithm, we repeatedly take advantage of the property that for any unique cycle, and the removal of its heaviest edge different from e produces the minimum weight sub- tree of T e containing e. the considered edge in T and we choose to include it in the sampled tree, T, we recursively update log (T) by adding the Bernoulli log probability corresponding to the edge. See Fig. for an illustration.
T' and edge e not in T' T U contains
3.3 The JC Sampler
The JC69 model is arguably one of the most frequently used substitution models in Bayesian phylogenetics. However, as far maximum likelihood estimates of branch lengths. Here, we observe that the JC69 model can be seen i..dd Bernoulli likelihood over sites and leverage the Deta-Bernoulli conjugacy for obtaining normalized distribution over branch lengths in variational inference. Below we show how to sample from this distribution and evaluate its likelihood using the change-of-we varriaalles formula.
aware, the model has only been used for obtaining
In [6], the candidate edges are selected from those that will make perfect matching, and the probability of each edge is uniform.
5

---

Given a branch length b(e) = b, the probability of the m th nucleotide nor mutating along edge (i.j) is, according to the JC69 model,
112, =12"
p = (1 + 4 1
(10)
As all mutations are equiprobable for this model, the probability of the nucleotide mutating is thus 1 p. giving rise to the following distribution when all M estimated mutations and non-mutations are observed and assuming mutations occur independently across sites
Bern (0t/p)
II Bern (oGjp)
II p 1l-oo (1 p) -p p))
(11)
where o is the expected number of mutations over the edge (i,j) at site m. Assuming an uninfor- mative Beta prior the Bernoulli parameter, p. and recalling that Beta is conjugate prior for the Bernoulli distribution, we get posterior distribution over p. Beta(p) M - &by + 1, tij + 1). Since we are ultimately interested sampling b, not p, we find a transform, f(-), mapping p to b, using Eq. (10);
1
3
f-1 (b)
+
=
(12)
4
4
f(p) = log
(13)
In all, sampling branch lengths from the JC model include two straightforward steps, (1) draw from the Beta posterior distribution, and (ii) pass it through the transform, b f(p).
Of course, we are also interested in evaluating the likelihood of the sampled branch length. To do so, use the change-of-variables formula; let
db)
(b) Beta(p|MM - D2s + 1, dig +1)
(14)
db
nd 0is p) uuse $b
(15)
B(M thij + ون + +1)
where B(-) is the Beta function. This provides an efficient way of evaluating the likelihood of branch length. In Alg. 2, we summarize the JC sampler with algorithmic description.
Algorithm The JC Sampler
1: Input: é, j.j
2: Sample p ~ Beta(p)
0iz + +1, 1) [Eq. (13))
3: Transform b == f()
4: Evaluate (b) (Eq. (14)
5: return b, (b)
f0
roject
accapt
accept e
elt) toe 35 [cos]
x(*)
elt) Fogl X .1 (spx:cap)
s(r) Fox Feu M (eps.cap.cer)
M (cos)
Figure 2: Propagation step of SLANTIS. The solid line is far, and bold green lines accepted edges. Dashed and dotted lines are the edges not in r. The blue color indicates the edge we consider (e), and the orange color indicates the alternative edge, either e* re (e))/WW/ + W()) or or = W (e))(W((e) - W W(e")). At the end of SLANTIS, the spanning tree with bold green lines is reported.

---

3,4 Assessing the Framework
The framework learss its parameters and distributions by maximizing the ELBO in Eq. (4). However, the IWELBO offers tighter lower bound on the marginal log-likeling leglikelihood. Using SLANTIS and the JC sampler, we can easily evaluate our framework using the IWELBO:
Σ po B TE)
L1 (.,,7))) log L 80 (Br, Tr)
(16)
where we may factorize the denominator as (B( = sy (b))) and (x(B)) = [eEEE()) 86 (b(e)). Note that have marginalized out the auxiliary variable, Z, in the target in Eq. (16), do not wish to compute L1 using our variational distributions.
4 Using VaiPhy to Parameterize the CSMC
As the tree topology space in A is considerably larger than the bifurcating X treee space, obtaining competitive loglliielhoodd (LLL) estimates not straightforward. Namely, the uniform prior over the tree space in A, (T) 1/(N - |X turns into heavily penalizing term. When compared to baselines in smaller sample spaces, this will become a considerable negative factor for VaiPhy. Compare this to the prior when considering unrooted, bifurcating X X-trees (which is the case for VBPI-NF): (T) = 1/(2|X - 5)!!, where X| is the number of taxa. As an example, for [X] = 10 and N 18, (r) distributes mass uniformly over 281 trillion trees, while the latter one merely considers million trees.
Inspired by [26], we instead seek to use VaiPhy to parameterize CSMC algorithm, thereby obtaining a procedure for projecting our trees to the lower-dimensional leaf-labeled tree space. Here, the prior is instead (r) = X/ - 3)!!.
Since is not constrained to multifurcating trees, we utilize the estimated number of mutations between vertices learned in A, also in the bifurcating X-tree space, Q. Here arises an interesting obstacle. In the CSMC algorithm, the internal vertices are unlabeled. Hence we are effectively blindfolded and need to estimate which entry in á to use for a given edge. We achieve this by scoring edges based on how frequently the corresponding split of the leaf-set X occurs in trees sampled from (T). See Appendix D.2 for more information.
At each rank p R). the CSMC proposes (i) forest, Tp- by merging the roots of two subtrees in T,-1, and (ii) a corresponding set branch lengths, B,. This done in parallel for K particles, resampling particles with probability proportional to their importance weights, w at every rank.' However, the reason to take advantage of a CSMC is that it allows to estimate the marginal log-likelihood in the bifurcating X-tree space using the well-known formula,
R
K
logpa (appa() (X) RJ log II K
p
(17)
We let é parameterize our two novel CSMC proposal distributions, r. (B, and (T,
The former constructs a mixture of JC samplers and evaluates the proposal's 's likelihood a multiple importance sampling fashion [7], while the latter forms categorical distribution over the roots in To-1. parameterized by a function of the likelihoods of trees sampled by SLANTIS. The in-detail descriptions can be found in Appendix D.
5
Experiments
5.1 Visualizing the JC Sampler
We conduct toy experiment, investigating how the estimated number of mutations along edge affects the behavior of the JC sampler. We would expect that greater number of mutations should
For more details and background of the CSMC algorithm, we point the reader to [26] and [33] for excellent descriptions.
7

---

generate a larger branch length. This is the case since the involved vertices are supposedly farther apart in terms of an evolutionary distance.
As demonstrated in Fig. 1, this advantageous behavior is indeed displayed. As Φi grows, the corresponding JC sampling distribution is shifted towards larger branch lengths. Interestingly, the JC sampler distributes its mass more symmetrically than the Log-Normal or Exponential distributions, which are typically used for modeling branch length distributions [33, 39, 37, 26].
follows from Eq. (13) that branch lengths grow infinitely large as p goes to 0.25. This phenomenon
in phylogenetics is referred to saturation. Saturation occurs for distantly related lineages and implies that the states of the nucleotides in these ancestral sequences are mutually independent. This not an interesting case, hence the -xhauusautivee parameter search of Φ.j in Fig. 1.
5.2 Density Estimation
Here we benchmark our methods, VaiPhy, and 6-CSMC, in terms of LL estimates on seven datasets, which we refer to as DS1-DS7 ([16, 13, 35, 17, 25, 40, 29]; in Appendix E, we provide additional information about the datasets). Additionally, to the extent possible, we have used the open-source code of our baselines to survey the benchmarks on these popular datasets. We only considered baselines that employ the JC69 model. Moreover, we assumed an Exp(10) prior over branch lengths and a uniform prior over tree topologies. This was aligned with the generative models used in all the original works of the baselines.
For all methods, we use the exact parameter configurations reported in the corresponding papers or specified in their available code. We let VBPI-NF [37] and VCSMC [26] train for the number of iterations as specified in their works. For MrBayes [28] with the SS method, we follow the protocol used in [39].
As benchmarking sophisticated versions of the phylogenetic inference CSMC [33], it seems relevant to include a vanilla CSMC. We define the vanilla CSMC as one that proposes branch lengths using the model prior and merges subtrees uniformly. There are no learnable parameters in the vanilla CSMC.
Indeed, there exist particle MCMC based approaches for parameterizing the CSMC, c.g., [33] or [34]; however, these have not previously been benchmarked on the datasets used here, and the latter uses the K2P substitution model. As such, we leave these comparisons for future work and constrain ourselves to VI based CSMCs.
Following [26], we give all CSMC methods K = 2048 particles. The parameters of VCSMC, VaiPhy, and @-CSMC are selected based on their best iteration, measured LL. Using these parameters, we compute the LL estimates averaged over ten random seeds. We run VaiPhy and d-CSMC for 200 iterations and evaluate them on Eq. (16) and Eq. (17), respectively. Vanilla CSMC is also evaluated using Eq. (17).
In Table 1, we provide the mean LL scores and standard deviations. On all datasets except DS2, our é-CSMC is the superior CSMC method (highlighted in red). It also produces significantly smaller variations its LL estimates compared to VCSMC and the Vanilla CSMC. On DS2, the Vanilla CSMC is surprisingly the winning CSMC.
Overall, MrBayes SS is the superior method these datasets, followed by VBPI-NE. VaiPhy, despite its heavily penalizing prior distribution over tree topologies (see Sec. 4), still outperforms VCSMC and Vanilla CSMC on all datasets but DS2, with extremely low-variance estimates of the LL.
5.3 Wall-Clock Time and Memory Comparisons
Being mean-field VI algorithm, VaiPhy truly enjoys the related benefits in terms of speed, which we demonstrate in this section. The experiments are performed on high-performance computing cluster node with two Intel Xeon Gold 6130 CPUs with 16 CPU cores each. Each node in the cluster has 96 GiB RAM.
We run VaiPhy for 200 iterations, and 128 trees are sampled using SLANTIS to compute the expectations in Eq. (6) and (8). To compute Lr., 3000 trees are sampled from SLANTIS after 200 iterations. The G-CSMC results include the aforementioned VaiPhy runs and 10 additional

---

Table 1: Real dataset LL estimates for the different methods. Our *-CSMC is, average, the best-performing CSMC method (red). Overall, MrBayes SS achieves the highest LL estimates (bold). *Not available due unresolved memory issues.
DATA SET
(#TAXA, 4SITES)
MRBAYsS
VBPLNF
(CCMC
VANILLA CSMC
VAIPHY
(-CSMC
DS1
(27. 1949)
-7033.45
0.15
7108
0.11
9180.34 170.27
4306.76 d: 166.27
7490.54
72290
7.25
D52
(29, 2520)
266363.55 : 0.35
26367 70 0.03
28700.7 4892.67
27884 37 d: 226.26
31203:44
3E
30568.49
11,34
DS3
(36, 18:2)
-33729.60
0.74
-33735 09
0.05
-37211 22 397.97
-35381.0 01 + 18
-33911.13
1011113::0 7E
-33798.06 6.62
DS4
(41, 1137)
-132922.37
1.42
-13329-99
2 0.11
-17106. 10 362.74
-15019.21 100.61
-13700 86
-13582.24
35.08
DSS
(90, 5781
-8191.96
0.71
-8214 61
0.40
-9499 65 1578.58 55
-8940 62 46.44
-8864 77 =00
-8367.51 5- 8.87
DS5
(40. 1133)
-6871.82
0.81
-6724
16..0.37
-9196.66
3046,
-8039 51 83.67
-7159
-901%.8%
99
D57
(64. 1008)
-$438.78
73,66
-8650
49 0.53
NIA"
-11013
57 113.49
-9460.21 21
1F
.9309.18 18.03
independent CSMC runs with 2048 particles. VBPI-NF runs for 400,000 iterations and uses 100,000 trees to estimate the lower bound after training. Each experiment is repeated 10 times. The wall-clock time comparison of the methods on DS1 is presented Fig. 3.
80000
Wal Walkk time (6e6)
Wall clock time (snc)
100000
20000
80000
Ve Phy
15000
60000
9-CMMC
VBR-NF
10000
40000
20000
5000
8
2
8
16
Number of cores
Number of cores
Figure 3: Wall-clock runtime results on DS1 data set with (2, (2,4,8, 16) cores. Left) Runtime results for VaiPhy, 6-CSMC, and VBPI-NF. Right) A zoomed-in version of the top figure showing runtime results for VaiPhy and şi D-C. Clearly, VaiPhy consistently faster than its benchmarks.
We were not able to run VCSMC on less than 16 cores due to the method producing memory issues on DS1. However, for 16 cores, their runtime was around 6840 seconds, with too small variations to show in a boxplot. VCSMC is thus faster than é-CSMC on 16 cores but still slower than VaiPhy.
VaiPhy is remarkably faster than other variational methods for each number of cores investigated. Although VBPI-NF produces impressive I.I. estimates, it is considerably slower than VaiPhy and the CSMC methods (all CSMCs are not shown in the figure). Additionally, using two cores, the wall-clock time required by VBPI-NF to arrive at VaiPhy's approximation accuracy was close to 11x that of VaiPhy DSI. For completeness, we compared wall-clock times with MrBayes, although this is not appropriate given that the software has been highly optimized over decades and is written in And, indeed, MrBayes is the fastest algorithm when allowed multiple cores. However, for single core on DS1, MrBayes ran for approximately 100 minutes while VaiPhy only needed an impressive 83 minutes on average for the same setup.
Finally, we discuss memory requirements. Our main variational baseline, VBPI-NE, used 8.4 million neural network weights and biases DS1, while the number of learnable parameters (i.e., the dimensions of & and x() in VaiPhy was 0.8 million; less than 10% of the former.
6 Conclusion
We have presented modular framework for Bayesian phylogenetic inference, building on three novel methods: VaiPhy, SLANTIS, and the JC sampler. At the heart of our framework lies ø, the estimated number of mutations matrix. Although & was learned in multifurcating tree space with labeled internal vertices, we showed that the evolutionary distances also served useful in the bifurcating X-tree space. demonstrated this by using & to parameterize two new proposal distributions for the CSMC, also introduced by us in this work; we referred to this new algorithm d-CSMC. Our 6-CSMC proved to the best-performing CSMC in terms of marginal aigliiklldood estimation on real datasets. Additionally, VaiPhy up to 25 times faster than its basclines on real data.

---

Speed is essential for making Bayesian phylogenetics a reality for future large data sets, in particular, single-cell data sets in medical investigations. Moreover, speed and the opportunity to leverage the VI methodology for phylogeny, as provided by our framework, i.e., our graph-based representation of the VI distribution over trees, Φ. opens up novel research opportunities. For instance, recently, ensemble and mixture models in VI have received a lot of attention [23, 27]. VaiPhy could obtain many ensemble components in parallel with low computational cost and thus has great potential to benefit from these findings. Moreover, more mixture components appear to imply significant improvements in terms of negative log-likelihood scores when learned by maximizing a multiple importance sampling- bbsen objective, MISELBO [24]. This motivates researching how VaiPhy could be fitted using a MISELBO objective.
For many potential applications, models other than the JC69 model will be more suitable for the evolution that takes place over the edges of the phylogenetic tree. However, our JC sampler provides great starting point. Reasonably, the JC sampler should serve as better base distribution for normalizing flows than those currently employed [37]. Furthermore, in this work, the JC sampler samples the branch lengths independent of the tree topology. We expect that tree topology dependent branch sampling will enhance VaiPhy's performance. This will be investigated in our future work.
VI can facilitate Bayesian analysis of models that integrate phylogeny with other tasks that usually are performed as pee-ppreeeseeeeeee pre-prceessig or post-pnn prosessocessing steps. The most important such example may very well be sequence aligmment. Although Bayesian phylogenetic inference is currently considered state-of-the-art this inference does not take the observed sequences associated with the taxa as input but a so-called multiple sequence alignment (MSA) of these, which is obtained by applying a heuristic to the sequences. This implies that the entire uncertainty related to the multitude of possible MSAs remains uncharacterized after a standard Bayesian phylogenetic inference, i.e., it does not affect the obtained posterior at all. For single-cell phylogenetic analysis of cancer cells, pre-processing may be even more involved. For instance, the so-called direct library preparation single-cell DNA data [36] for tumors require a pre-processing step in which the cells are clustered into cancer clones. For these clones, it is then desirable to perform a phylogenetic inference. Interestingly, the VI has already been launched for Bayesian clustering for this pre- processing step [30].
In conclusion, we have presented a VI algorithm for Bayesian phylogeny that, together with its building blocks, unlocks the potential of VI for fast Bayesian analysis in phylogeny and integrated phylogenetic analysis. The by now 25 year old methodology of MCMC based Bayesian phylogenetic analysis is certainly mature and optimized methodology with great accuracy. However, VaiPhy is faster than MrBayes of increasingly faster and more accurate VI based phylogenetic software. Moreover, it paves the way for integrated VI based phylogenetic inference.
single core, and the underlying methodology is bound to lead to sequence
Acknowledgments and Disclosure of Funding
First, we acknowledge the insightful comments provided by the reviewers, which have helped improve our work. This project was made possible through funding from the Swedish Foundation for Strategic Research grant BD15-0043, and from the Swedish Research Council grant 2018-05417 VR. The computations and data handling were enabled by resources provided by the Swedish National Infrastructure for Computing (SNIC), partially funded by the Swedish Research Council through grant agreement no. 2018-05973.
The authors declare no competing interests.
References
[1) D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. Journal of the American statistical Association, 112(518):859-87) 2017.
[2] A. Bouchard-Côtt and M. Jordan. Variational inference over combinatorial spaces. In Advances in Neural Information Processing Systems, volume 23, 2010.
[3] A. Bouchard-Cötéé S. Sankararaman, and M. Jordan. Phylogenetic inference via sequential Monte Carlo. Systematic Biology, 61(4):579-593, July 2012.
10

---

[4] Y. Burda, R. B. Grosse, and R. Salakhutdinov. Importance Weighted Autoencoders. In The International Conference on Learning Representations, 2016.
[5] T. Dang and H. Kishino. Stochastic Variational Inference for Bayesian Phylogeneties: A Case of CAT Model. Molecular Biology and Evolution, 36(4).msz020, 2019.
[6] P. Diaconis. Sequential importance sampling for estimating the number of perfect matchings in bipartite graphs: An ongoing conversation with laci. In Building Bridges II, pages 223-233. Springer, 2019.
[7] V. Elvira, L. Martino, D. Luengo, and M. F. Bugallo. Generalized multiple importance sampling. Statistical Science, 34(1):129-153: 2019.
[8] Y. Fan, R. Wu, M.-H. Chen, L. Kuo, and P. O. Lewis. Choosing among Partition Models in Bayesian Phylogenetics. Molecular Biology and Evolution, 28( (8)))5523553, 2011.
[9] J. Felsenstein. Evolutionary trees from dna sequences: a maximum likelihood approach. Journal of molecular evolution, 1766)336.-376 1981.
[10] R. W. Floyd. Algorithm 97: Shortest path. Commun, ACM, 5(6):345, jun 1962.
[11] M. Fourment and A. E. Darling. Evaluating probabilistic programming and fast variational Bayesian inference in phylogenetics. PeerJ, 711):2822, 2019.
[12] N. Friedman, M. Ninio, Pe'er, and T. Pupko. A structural EM algorithm for phylogenetic inference. Journal of computational biology: a journal of computational molecular cell biology,
9())3311-33,, 2002.
[13] J. R. Garey, T. J. Near, M. R. Nonnemacher, and S. A. Nadler. Molecular evidence for acanthocephala as subtaxon of rotifera. Journal of Molecular Evolution, 43(3):287-292, 1996.
[14] O. Gascuel. BIONJ: improved version of the NJ algorithm based simple model of sequence data. Molecular Biology and Evolution, 14(7):685-69 008855695,00 07 1997.
[15] S. Guindon, J.-F. Dufayard, V. Lefort, M. Anisimova, W. Hordijk, and O. Gascuel. New Algo- rithmss and Methods Estimate Maximun-Likelihood Likelihood Phylogenies: Assessing the Performance of PhyML 3.0. Systematic Biology, 59(3):307 05 2010.
[16] S. B. Hedges, K. D. Moberg, and L. R. Maxson. Tetrapod phylogeny inferred from 18s and 28s ribosomal ma sequences and review of the evidence for amniote relationships. Molecular Biology and Evolution, 7(6))607-633, 1990.
[17] D. A. Henk, A. Weir, and M. Blackwell. Laboulbeniopsis termitarius, an ectoparasite of termites newly recognized
member of the laboulbeniomycetts.. Mycologia, 95(4):561-564, 2003,
[18] D. T. Hoang, O. Chernomor, A. von Haeseler, B. Q. Minh, and L. S. Vinh. UFBoot2: Improving the Ultrafast Bootstrap Approximation. Molecular Biology and Evolution, 35(2):518-535 10 2017.
[19] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 2013.
[20] S. Höhna, M. Landis, T. A. Heath, B. Boussau, N. Lartillot, B. R. Moore, J. P. Huelsenbeck, and F. Ronquist. Revbayes: Bayesian phylogenetic inference using graphical models and an interactive model-specifican -pecifficatinn language. Systematic biology, (5(4))726-733)) 2016.
[21] J. P. Huelsenbeck and Ronquist. Mrbayes: Bayesian inference of phylogenetic trees. Bioin- formatics, 17(8):754-755, 2001.
[22] T. Jukes and C. Cantor. Evolution of protein molecules, pp. 21-132 in mammalian protein metabolism, vol. ii munro hn, 1969.
[23] O. Kviman, H. Melin, H. Koptagel, V. Elvira, and J. Lagergren. Multiple importance sampling elbo and deep ensembles of variational approximations. In International Conference on Artificial Intelligence and Statistics, pages 10687-10702. PMLR, 2022.
[24] O. Kviman, R. Molén, A. Hotti, S. Kurt, Elvira, and J. Lagergren. Learning with miselbo: The mixture cookbook. arXiv, 2022.
[25] C. Lakner, P. Van Der Mark, J. P. Huelsenbeck, B. Larget, and F. Ronquist. Efficiency of markov chain monte carlo tree proposals in bayesian phylogeneties. Systematic biology, 57(1):86-11)) 2008.
11

---

[26] A. K. Moretti, L. Zhang, C. A. Naesseth, H. Venner, D. Blei, and I. Pe'er. Variational combina- torial sequential monte carlo methods for bayesian phylogenetic inference. In Uncertainty in Artificial Intelligence, pages 971-981. PMLR, 2021.
[27] W. Morningstar, S. Vikram, C. Ham, A. Gallagher, and J. Dillon. Automatic differentiation variational inference with mixtures. In International Conference on Artificial Intelligence and Statistics, pages 3250-3258. PMLR, 2021.
[28] F. Ronquist, M. Teslenko. Van Der Mark, D. D. Ayres, A. Darling, S. Höhna, B. Larget, L. Liu, M. A. Suchard, and J. P. Huelsenbeck. Mrbayes 3.2: efficient bayesian phylogenetic inference and model choice across large model space. Systematic biology, 61(3): 2012.
[29] A. Y. Rossman, J. M. McKemy, R. A. Pardo-Schultheiss, and H.-J. Schroers. Molecular studies of the bionectriaceae using large subunit rdna sequences. Mycologia, 93(1):100- 10, 2001.
[30] N. Safinianaini, C. P. de Souza, and J. Lagergren. CopyMix: mixture model based single-cell clustering and copy number profiling using variational inference. bioRxiv, 2020.
[31] R. A. Vos, H. Lapp, W. H. Piel, and V. Tannen. Treebase2: rise of the machines. Nature Precedings, 2010.
[32] M. J. Wainwright and M. I. Jordan. Graphical Models, Exponential Families, and Variational Inference. Foundations and Trends@ in Machine Learning, 111-2):1-305, 2007.
[33] L. Wang, A. Bouchard-Côté and A. Doucet. Bayesian phylogenetic inference using com- binatorial sequential monte carlo method. Journal of the American Statistical Association. (106112):132-1377 2015.
[34] S. Wang and L.. Wang. Particle gibbs sampling for bayesian phylogenetic inference. Bioinfor- matics, 37(5):642-644 2021.
[35] 7.. 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 a radiation of cute-looking mouse lemur species. Systematic biology, 52(5):705-716, 2003.
[36] H. Zahn, A. Steif, E. Laks, P. Eirew, M. VanInsberghe, S. P. Shah, S. Aparicio, and C. L. Hansen. Scalable whole-genome single-cell library preparation without preamplification. Nature methods, 14(2):167-177) Feb. 2017.
[37] C. Zhang. Improved variational bayesian phylogenetic inference with normalizing flows. Advances in Neural Information Processing Systems, 33:18760-18771, 2020.
[38] C. Zhang and F. A. Matsen IV. Generalizing tree probability estimation via bayesian networks. Advances in Neural Information Processing Systems, 31, 2018
[39] C. Zhang and F. Matsen IV. Variational bayesian phylogenetic inference. In International Conference on Learning Representations, 2018
[40] N. Zhang and M. Blackwell. Molecular phylogeny of dogwood anthracnose fungus ( discula destructiva and the diaporthales. Mycologia, 93:355-365, 03 2001.
12

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
(b) Did you describe the limitations of your work? [Yes] We discuss the extension to
labelled multifurcating tree space and the requirement of the post-processing step that is required to make fair comparisons in bifurcating space.
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes] We show the conjugacy of JC sampler and we provide proof of natural gradient in the Appendix.
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either in the supplemental material as URL)? [Yes] We will provide the codes the supplementary materials.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] In the Experiments section we mention the choice of parameters. In Appendix, we provide code segments to show how we ran the benchmark methods.
(c) Did you report error bars (e.g., with respect to the random seed after running exper- iments multiple times)? [Yes] In Table we include the standard deviations and in Figure we include the error bars.
(d) Did you include the total amount of compute and the type of resources used (e.g.. type of GPUs, internal cluster, or cloud provider)? [Yes] We wrote the hardware specifications in the runtime experiments; however, didn't 't include the total amount of time required for all datasets.
4. If you are using existing assets (e.g., code, data, models) or curating/rellasing new assets...
(a) your work uses existing assets, did you cite the creators? [Yes] We cited all software and also provided GitHub links to the repositories in the Appendix.
(b) Did you mention the license of the assets? [N/A]
(c) Did you include any new assets either in the supplemental material URL? [Yes] The code is available in the supplementary material.
(d) Did you discuss whether and how consent was obtained from people whose data you're rou're
using/curating? [N/A]
(e) Did you discuss whether the data you are using/curating contains personally identifiable information offensive content? [N/A]
5. If you used crowdsourcing conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid participants and the total amount spent on participant compensation? [N/A]
13


Reviews:
[{"review_id": "lAALlBpw9dI", "reviewer": "Reviewer_aBtL", "summary": "The paper introduces two methods for inference of phylogenetic trees from sequence data. Taking a variational inference approach, the VaiPhy method is constructed using two new samplers: SLANTIS and JC. The ability to estimate the expected number of mutations over edges leads to a variant of the CSMC method denoted \\phi-CSMC. The methods are experimentally tested and compared to state-of-the-art methods on 7 datasets.", "questions": "* The sentence \"VaiPhy is remarkably much faster than all other methods for all numbers of cores (except for\n278 MrBayes in special cases, see further down)\" seems a little strong when it is followed by \"MrBayes is the fastest algorithm when allowed multiple cores.\" Can the first sentence can be expressed more clearly?", "limitations": "yes", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* introduction of several new techniques for variational inference in phylogenetic analysis\n* important, high impact application area\n* well-presented and technically sound paper (as far as I can evaluate)\n* while the two methods only outperform existing methods in some cases, the techniques show a path to future advances with potential to significantly improve the state-of-the-art\n\nWeaknesses:\n* variational inference has been used before for phylogenetic analysis, so in that sense the underlying idea is not fundamentally new. However, as far as I can judge, the proposed methods still presents significant novelty and originality", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Pa23lg4BYHX", "reviewer": "Reviewer_Mtqk", "summary": "The authors infer the posterior distribution over phylogenetic trees given only the leaves of the tree. To do so, the authors propose a method composed of three elements. A first one (VaiPhy) is a mean-field VI algorithm that returns the posterior distribution of the tree, internal vertices and set of branch lengths. VaiPhy uses the second element (JC sampler) to sample branch lengths and the third element (SLANTIS) to sample trees.  SLANTIS enables them to sample the branches of the tree one after the other, with probabilities proportional to the total weight of the maximum spanning tree. The JC sampler gives access to a normalised distribution over branch lengths, using the Jukes-Cantor model.\n\nThe authors compare VaiPhy to several other models, baselines and a VaiPhy-parametrised CSCM, using the IWELBO as metrics and wall-clock runtime. The authors show that among their proposed variants, VaiPhy-parametrised CSCM give the best results. While it is not better than MrBayes or VBNI-NF, it is significantly faster than VBNI-NF.\n", "questions": "I believe that the paper is overall well-written and structured. I have nevertheless a few clarification questions concerning the method and questions regarding the benchmark. Additionally, the method MrBayes SS seems to perform better than the proposed approach. Therefore, I think that the work is not ready for publication at NeurIPS in its current form but I would be happy to increase the score if the following questions were answered.\n\n*****\n\n\n\n### VaiPhy\nOne step in VaiPhy is to compute the expected number of mutations of an edge. The authors claim that the maximum number of mutations between two vertices is the number of sites. I believe that it is possible to have unobserved mutations, for example when the time step between two vertices is large compared to the mutation rate. I would be interested to know if this is something the authors are accounting for, and how or why not if it is not taken into account. \n\nVaiPhy seems to consider that the nucleotides are independent. However, it is possible that some mutations co-occur or are mutually exclusive. Could these first order interactions be taken into account in VaiPhy?\n\n### SLANTIS\nThe separation between previous work and proposed approach is unclear for SLANTIS. The authors could perhaps make a better distinction between the two in the document.\n\nDue to the probabilistic nature of SLANTIS, is it possible that the leaves are not always all included in the reconstructed tree? If yes, how does this affect the performance?\n\n### Parameterized CSMC\nI believe that using VaiPhy to parameterize CSMC instead of using directly VaiPhy is impactful for a small number of taxa but could be less so for larger number of taxa. Have the authors compared the two models when the number of taxa is larger than 100?\n\n### Experiments\nI think that the experiments are interesting but I would like to better understand how the proposed approach performs and scales in terms of number of taxa or sites for example compared to the baselines. \n\nFigure 1, which compares the number of mutations with branch length, is an interesting plot. I was wondering if it would be possible to also evaluate the reconstructed tree. This could be done on simulated data where the tree (and sequences of the ancestors) are known but only the leaves are used for phylogeny reconstruction. \n\n\n\n\n**Minor comments**\n- Line 4: change the word formidable for something less positive \n- Line 14: \"the, to the best\" -> phrasing to change\n- Line 22, 102: \"phylogentic\" -> phylogenetic\n- Line 76: define CSMC\n- Equation 6: define the subscript -i\n- Equation 7: the right hand side is missing the dependency in T\n- Equation 9: it is worth repeating what $Q_{\\alpha\\beta}$ is\n- line 146: Algorithm B -> algorithm 3\n- The sentence line 153-154 is very unclear, it is worth reformulating it.\n\n---------\nIn light of the authors' response to the reviewers, I am in favor of accepting the paper and increase the grade to 6.\n\n", "limitations": "I think that the authors could better address the limitations of their work in a few sentences.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Originality\nThe methods proposed by the authors is novel as it combines and adapts several techniques to a specific problem, which is phylogeny inference in single-cell data. \n\n### Significance\nPhylogeny inference is an important problem and is in need of approaches that scale to the size of current single-cell datasets for example. The authors do a step in this direction. However, MrBayes SS seems to both perform better and is faster than the proposed approach.\n\n### Clarity \nThe paper is overall well-written, apart from a few sentences / typos that are mentioned at the next section. \n\n### Quality\nWhile the overall quality of the paper is good, including maths and structure of the paper, I think that some points need to be better explained. They are the subject of the next section. \n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "PU9NJ1BbSzs", "reviewer": "Reviewer_vFFt", "summary": "This paper contributes to the existing Bayesian phylogenetic VI literature with three novel methods which, in combination, lead to better accuracy, lower memory requirements, and faster computation speed when compared to existing VI methods in the space. We now discuss each of the authors' three novel contributions in isolation. \n\nSLANTIS, a method for generating proposal trees in the multifurcating tree space, is the first contribution. As reported, this method is orders of magnitude faster than existing methods for generating proposal trees, taking far fewer iterations to approximate the relevant expectations than existing methods.\n\nJC sampler uses the computed matrix of expected mutations across branches to sample branch lengths exactly according to the Jukes Cantor model. Sampling branch lengths directly from the underlying substitution model instead of from an arbitrary prior just makes sense, and frankly I'm surprised that this is the first paper to do so.\n\nVaiPhy is the full variational inference method developed in this paper, utilizing both SLANTIS and JC sampler. VaiPhy enjoys significantly lower computational times than existing VI methods, while consistently producing better log-likelihood scores. Additionally, VaiPhy includes an auxiliary variable for ancestral sequences, which helps pave the way for potential future applications which may want to include ancestral data, such as that of fossils.\n\nBayesian phylogenetic VI is a new and growing methodology, so, while the methods in this paper do not produce the best results among all phylogenetic inference methods, they provide significant contributions towards the development of VI in phylogenetics.", "questions": "Section 3.4 states that \"the framework learns its parameters and distributions by maximizing the ELBO [... but] the IWELBO offers a tighter lower bound\" and \"using SLANTIS and the JC sampler we can easily evaluate our framework using the IWELBO.\" This statement is confusing: the former states that the ELBO is used in the framework, but the latter seems to imply that the IWELBO is used instead; which is true?\n\nHave you considered utilizing ancestral sequence data (e.g. from fossils) as inputs into the method, instead of featuring them exclusively as auxillary variables?\n\nWhy were convergence statistics for SLANTIS not included? When proposing a method which proports to require orders of magnitude fewer iterations than existing methods, convergence seems an important attribute to investigate.\n\nThe method improves on memory requirementns of previous VI implementations by utilizing significantly fewer learnable parameters, but what amount of memory does the method actually use at runtime? I.e. how many GiB of ram does the method utilize? This question is very relevant to the current limitations of the method.", "limitations": "Of all phylogenetic tree inference methods, this is neither the fastest nor the most accurate. The authors very honestly admit this, and though this is a limitation, it is should _not_ be considered a weakness of the paper, as  VI in Bayesian inference is still in its infancy, while their competitor MrBayes is comparatively a full grown adult. That being said, MrBayes should not have been excluded from the wallclock time vs. number of cores plots, as how far modern VI methods lag behind MrBayes is not irrelevant.\n\nA very typical assumption to make in phylogenetic likelihood computations is independence of branch lengths. Computations without this assumption are rarely feasible. This method calculates branch lengths in isolation, though, without respect to topology, which I expect is far too much of a simplifying assumption.\n\nLog-likelihood values for the method (as well as its competitors) are provided on 7 real-world datasets, but this metric is not actually very good for evaluating the performance of phylogenetic inference. Better log-likelihood scores do not _necessarily_ equate to a better estimate of the underlying tree, especially when 1) many simplifying assumptions are made in likelihood calculations, and 2) underlying models (such as the substitution model, JC in this case) may very well be misspecified. A better way to evaluate the performance of these methods is through a simulation study where the ground truth trees are known, and prediction accuracy can then be calculated via the Robinson-Foulds distance.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "## Originality\n\nWhile the classical VI methodology that underlies this paper's work is nothing new, SLANTIS and JC sampler are both completely novel methods. Additonally, ancestral sequence data has been included in previous Bayesian methods as part of the input, not just as auxillary variables. However, VaiPhy is the first VI-based method to take ancestral data into consideration, even if only probabilistically. Thus, this paper contributes a significant body of original material towards the improvement of VI-based phylogenetic inference.\n\n## Quality\n\nThe overall quality of the methods presented is very high. Further details are contained in other relevant sections.\n\n## Clarity\n\nThe methodology is explained very well; as someone familiar both with phylogenetics and variational inference I felt that I had a very good understanding of the mechanisms behind VaiPhy on my first read, which is not an easy task for a paper on something as computationally heavy as VI.\n    \nThe background section of the paper is devoid of any explanation of sequential Monte Carlo or conditional sequential Monte Carlo, which is a barrier for readers who are interested in the potential phylogenetic applications of the work, but unfamiliar with these concepts. Additionally, while SMC is at least given an introduction, CSMC is not (in fact the word \"conditional\" does not appear in the paper). While the paper is no doubt intended to appeal to an audience which is already familiar with variational inference, its applications are in phylogenetics, so these concepts should still be given an adequate introduction and brief explanation.\n\n## Significance\n\nVI is a developing field in phylogenetic inference. I do not know whether VI will end up leading the cutting edge in phylogenetics, but, irregardless, this paper's contributions certainly make that reality seem more possible.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "F53-8uER4NV", "reviewer": "Reviewer_k9KB", "summary": "Inferring the posteriors of phylogenetic trees is critical to evaluating uncertainty in tree inference. The most popular method for this is MCMC; however, in some situations effective sampling of the combinatorial tree space can be difficult, and MCMC does not return (straightforwardly) densities for individual trees, meaning that it can have difficulties representing diffuse posteriors. Variational inference, in principle, can solve these issues while decreasing the computational cost of inference. However current methods require lots of resources for training due to the need to perform automatic differentiation.\n\nIn this paper, the authors aim to improve the efficiency of variational inference. They make a number of contributions to do so:\n1) They introduce an efficient algorithm to sample trees conditional on branch lengths, substitution probabilities, and sequences at every vertex.\n2) For a Jukes-Cantor substitution model, they introduce an algorithm to sample branch lengths given the topology of the tree and sequences at every vertex.\n3) They introduce an algorithm for iteratively updating the posterior of the sequences at each vertex of the tree, the tree topology, and the branch lengths. They calculate the natural gradients for this algorithm to make training as efficient as possible.\n4) Finally, they also combine 1 and 2 to create a more efficient proposal for a tree sampling procedure based on stochastic monte carlo (SMC).\n\nThey demonstrate the efficiency of their method, taking substantially less time than another variational inference method and some sampling methods. They also demonstrate that their method in 4 indeed increases the efficiency of sampling.", "questions": "1) For the transformed branch length, of course when there are many sites and branch lengths are long this won't matter, what do you do when p<1/4? Do you account for p's that cannot give rise to valid branch lengths when calculating probabilities?\n2) Why did you compare to VCSMC instead of VNCSMC, the improved model? And did you perform variational inference on the parameters in phi-CSMC?\n3) Did you use the same branch length priors for all models (I believe other models have exponential branch length priors by default)?\n4) The likelihood estimation methods used for the VI and SMC methods are biased downwards and so not directly comparable. This lower bound is still useful to report for the SMC methods as it demonstrates more efficient sampling when it is higher, but this isn't a fair comparison to VSMC I believe which also optimises the parameters of its branch length priors and substitution matrix. However this isn't a fair comparison to MCMC which doesn't have a downward-biased marginal likelihood estimate. In other phylogenetics VI papers, samples are usually taken from the VI methods and their marginal likelihood calculated unbiasedly just as with the samples from MCMC. This could also be done for your method I believe, allowing evaluation without the inference of the sequences at internal nodes affecting the bias.\n", "limitations": "1) The mean field approximation to internal nodes may severely limit the ability of this model to fit certain posteriors.\n2) This method requires a very efficient way to sample branch lengths, which the authors demonstrate exists for the JC model, but may be more difficult to devise for more complicated substitution models.", "rating": 7, "confidence": 3, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n1) The paper proposes an extremely efficient method for variational inference of phylogenetic trees.\n2) In the process of devising this method, the authors propose a new method to sample trees topologies and brnch lengths that is of independent interest.\nThese points make it a substantial contribution to phylogenetic inference\n\nWeaknesses:\n1) It is obvious that classical variational inference methods are extremely inefficient (when first reported these methods were shown to be more efficient that MCMC based on \"number of iterations needed to reach a minimum accuracy\", however, as iterations come with different computational costs, the method the authors used to measure efficiency here is more sound.) However the number of iterations to run the classical method in this case seem to be arbitrary, or selected to be very large in the original publication to get to convergence. It would be more fair to compare compute time to log likelihood achieved.\nThis concern is more dire in the comparison between CSMC methods, the author's proposed method and MCMC - although not strictly necessary to demonstrate a contribution as one obtains different things from each methods.\n2) I am not 100% on the likelihood evaluations (see questions below).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "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": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept