ARTree: A Deep Autoregressive Model for Phylogenetic Inference
Authors: Tianyu Xie
Cheng Zhang
Venue: NeurIPS 2023 spotlight (2023)
Keywords: phylogenetic inference
autoregressive model
graph neural network
density estimation
variational inference
TL;DR: A novel deep autoregressive model for phylogenetic inference

Abstract: Designing flexible probabilistic models over tree topologies is important for developing efficient phylogenetic inference methods. To do that, previous works often leverage the similarity of tree topologies via hand-engineered heuristic features which would require domain expertise and may suffer from limited approximation capability. In this paper, we propose a deep autoregressive model for phylogenetic inference based on graph neural networks (GNNs), called ARTree. By decomposing a tree topology into a sequence of leaf node addition operations and modeling the involved conditional distributions based on learnable topological features via GNNs, ARTree can provide a rich family of distributions over tree topologies that have simple sampling algorithms, without using heuristic features. We demonstrate the effectiveness and efficiency of our method on a benchmark of challenging real data tree topology density estimation and variational Bayesian phylogenetic inference problems.

ARTree: A Deep Autoregressive Model for Phylogenetic Inference
Tianyu Xie', Cheng Zhang
School of Mathematical Sciences, Peking University
Center for Statistical Science, Peking University tianyuxie@pku. edu. cn, chengzhang@mathh pku. edu. cn
Abstract
Designing flexible probabilistic models over tree topologies is important for devel- oping efficient phylogenetic inference methods. To do that, previous works often leverage the similarity of tree topologies via hand-engineered heuristic features which would require pre-sampled approximation capability. In this paper, we propose a deep autoregressive model for phylogenetic inference based graph neural networks (GNNs), called ARTree. By decomposing tree topology into sequence of leaf node addition operations and modeling the involved conditional distributions based on learnable topological features via GNNs, ARTree can provide rich family of distributions over the entire tree topology space that have simple sampling algorithms and density estimation procedures, without using heuristic features. We demonstrate the effectiveness and efficiency of our method on benchmark of challenging real data tree topology density estimation and variational Bayesian phylogenetic inference problems.
topologies and may suffer from limited
Introduction
Reconstructing the evolutionary relationships among species has been one of the central problems in computational biology, with wide range of applications such as genomic epidemiology (Dudas
et al., 2017| du Plessis et al. 2021 Attwood et al. 2022) and conservation genetics (DeSalle & Ammaoo| 2004). Based on molecular sequence data (e.g. DNA, RNA, or protein sequences) of the observed species and model of evolution, this has been formulated as a statistical inference problem on the hypotheses of shared history, i.e., phylogenetic trees, where maximum likelihood and Bayesian approaches are the most popular methods (Felsenstein) 1981 Yang & Rannala 1997 Mau et al. 1999, Larget & Simon| 1999 Huelsenbeck et al.] (001), ) However, phylogenetic inference can be challenging due to the composite structure of tree space which contains both continuous and discrete components (e.g., the branch lengths and the tree topologies) and the large search space of tree topologies that explodes combinatorially the number of species increases (Whidden & Matsen IV
2013
Dinh et al., 2017).
Recently, several efforts have been made to improve the efficiency of phylogenetic inference algo- rithms by designing flexible probabilistic models over the tree topology space (Höhna Drummond] 2012: Larget 2013 Zhang & Matsen IV 2018). One typical example is subsplit Bayesian networks (SBNs) (Zhang & Matsen IV 2018)), which a powerful probabilistic graphical model that provides flexible family of distributions over tree topologies. Given a sample of tree topologies (e.g., sampled tree topologies from an MCMC run), SBNs have proved effective for accurate tree topology density estimation that generalizes beyond observed samples by leveraging the similarity of hand band-engineered subsplit structures among tree topologies. Moreover, SBNs also allow fast ancestral sampling and
Coorrspponding author.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

hence were later on integrated into a variational Bayesian phylogenetic inference (VBPI) framework to provide variational posteriors over tree topologies (Zhang Matsen [ 2019). However, due to the limited parent-child subsplit patterns in the observed samples, SBNs can not provide distributions whose support spans the entire tree topology space (Zhang Matsen IV 2022). Furthermore, when used as variational distributions over tree topologies in VBPI, SBNs often rely on subsplit support
estimation for variational parameterization, which requires high-quality pre-sampled tree topologies that would become challenging to obtain when the posterior is diffuse.
While SBNs suffer from the aforementioned limitations due to their hand-engineeree design, number of deep learning methods have been proposed for probabilistic modeling of graphs (Jin et al., 2018: You al.) 2018aj Cao & Kipf 2018; Simonovsky & Komodakis (018). Instead of using hand enginneeed features, these approaches use neural networks to define probabilistic models for the
connections between graph nodes which allow for learnable distributions over graphs. Due to the flexibility of neural networks, the resulting models are capable of learning complex graph patterns automatically. Among these deep graph models, graph autoregressive models (You al., 2018b) Li et al. 2018 Liao et al.) 2019 Dai et al. 2020; Shi et al. (020) are designed to learn flexible graph distributions that also allow easy sampling procedures by sequentially adding nodes and edges. Therefore, they serve an ideal substitution of SBNs for phylogenetic inference that can provide more expressive distributions over tree topologies.
In this paper, we propose novel deep autoregressive model for phylogenetic inference, called ARTree, which allows for more flexible distributions over tree topologies without using heuristic features than SBNs. With pre-selected order of leaf nodes (i.e., species or taxa), ARTree generates a tree topology by recursively adding new leaf nodes to the edges of the current tree topology, starting from star-shaped tree topology with the first three leaf nodes (Figure D The edge to which a new leaf node connects is determined according to a conditional distribution based on learnable topological features of the current tree topology via GNNs (Zhang 2023). This way, probability distributions provided by ARTree all have full support that spans the entire tree topology space. Unlike SBNs, ARTree can be readily used in VBPI without requiring subsplit support estimation for parameterization. In experiments, we show that ARTree outperforms SBNs on a benchmark of challenging real data tree topology density estimation and variational Bayesian phylogenetic inference problems.
2 Background
Phylogenetic likelihoods A phylogenetic tree commonly described by a bifurcating tree topology and the associated non-negative branch lengths q. The tree topology represents the evolutionary relationship of the species and the branch lengths q quantify the evolutionary intensity along the edges of T. The leaf nodes of 7 correspond to the observed species and the internal nodes of represent the unobserved ancestor species. A continuous time Markov model is often used to describe the transition probabilities of the characters along the edges of the tree (Felsenstein) (2004). Concretely, let Y = {YI...., Nur E Ω N M be the observed sequences (with characters in 12) of length M over N species. Under the assumption that different sites evolve independently and identically, the likelihood of Y given takes the form
M
M -I na 1=1
p(Y [, q) II p(
Pat. (guv), (u.))E
II
where a² ranges over all extensions of Y; to the internal nodes with at being the character assignment of node u (r represents the root node), E(T) is the set of edges of r, Que is the branch length of the edge (u, v) E E(r), (g) is the transition probability from character j to k through branch of length and 17 is the stationary distribution of the Markov model.
Subsplit Bayesian networks Let X' be the set of leaf labels representing the existing species. A non-empty subset of X' is called a clade and the set of all clades C(X') is equipped with total order > (e.g., lexicographical order). An ordered clade pair (W, Z) satisfying W n 2 = 0 and W > Z is called a subsplit. A subsplit Bayesian network (SBN) is then defined as Bayesian network whose nodes take subsplit values or singleton clade values that describe the local topological structures of tree topologies. For a rooted tree topology, one can find the corresponding node assignment of SBNs

---

by following its splitting processes (Figure in Appendix A). The SBN based probability of a rooted tree topology т then takes the following form
Psbu Pabu =) = p(S s1) II p(S. = = = sx,).
(2)
X=1
where St denotes the subsplit- singleton-clade-valeee random varaibles at node (node is the root node), Ty is the index set of the parents of node and {s: is the corresponding node assignment. For unrooted tree topologies. we can also define their SBN based probabílities by viewing them as rooted tree topologies with unobserved roots and integrating out the positions of the root node as follows: Psbn ( 7) = EeEE(+) Psbn( where 7 the resulting rooted tree topology when the rooting position on edge e. In practice. the conditional probability tables (CPTs) of SBNs are often parameterized based on sample of tree topologies (e.g., the observed data for density estimation (Zhang & Matsen IV 2018) or fast bootstrap/MCMC samples (Minh et al. 2013 Zhang 2020) for VBPI). As result, the supports of SBN-induced distributions are often limited by the splitting patterns in the observed samples and could not span the entire tree topology space (Zhang & Matsen IV 2022). More details SBNs can be found in Appendix A
Variational Bayesian phylogenetic inference posterior distribution takes the form
Given a prior distribution p(T, q), the phylogenetic
p(T, r p(Y q)((T x p(Y|T, ()p(T, q) p(Y)
Let (T) and (gl) be variational families over the spaces of tree topologies and branch lengths (T,
respectively. The VBPI approach uses ( 9) = (x(())))))) ( (g(T) to approximate the posterior p(T.q(Y) by maximizing the following multi-sample (K >> 1) lower bound
K
p( (1)
I. (6, xb = E
((rt.g 11 110 Qa.
log
Σ K
(4)
The tree topology distribution (T) is often SBNs which in this case rely on subsplit support estimation for parameterization that requires high-quality pre-sampled tree topologies and would become challenging for diffuse posteriors (Zhang & Matsen IV 2002). The branch lengths distribution (q)TT) can be diagonal lognormal distribution parametrized via heuristic features or learnable
topological features of v (Zhang & Matsen IV 2019 Zhang 2020 2023). Compared to the single- sample lower bound, the multi-sample lower bound in 44 allows efficient rrduued stochastic gradient estimators (e.g. VIMCO (Mnih & Rezende 2016)) for tree topology variational parameters. Moreover, using multiple samples would encourage exploration over the vast tree topology space, albeit it may also deteriorates the training of the variational approximation (Rainforth et al. (019). In practice, moderate K is often suggested (Zhang & Matsen IV 2022) See more details on VBPI in Appendix B
Graph autoregressive models By decomposing a graph sequence of components (nodes, edges, motifs, etc), graph autoregressive models generate the full graph by adding one component at a time, until some stopping criteria are satisfied (You et al. [2018b] Jin al., [0118] Liao et al. 2019). In previous works, recurrent neural networks (RNNs) for graphs are usually utilized to predict new graph components conditioned on the sub-graphs generated so far (You et al. 2018b). The key of graph autoregressive models is to find way to efficiently sequentialize graph structures, which is often domain- -specffcc.
3 Proposed method
In this section, we propose ARTree, a deep autoregressive model for phylogenetic inference that can provide flexible distributions whose support spans the entire tree topology space and can be naturally parameterized without using heuristic approaches such as subsplit support estimation. We first describe particular autoregressive generating process of phylogenetic tree topologies. We then develop powerful GNNs to parameterize learnable conditional distributions of this generating process. We consider unrooted tree topologies this section, but the method developed here can easily adapted to rooted tree topologies.
3

---

Ordinal tree topologies of higher rank
2
02 -
1 3 E
GNN
Sampling
GNN
Sampling
GNN
Sampling
1
DB
B B F Adding node
Leat
Interic
An overview of ARTree for autoregressive tree topology generation. The left plot the starting ordinal tree topology of rank 3. This tree topology is then fed into GNNs which output probability vector over edges. We then sample from the corresponding edge decision distribution and attach the next leaf node to the sampled edge. This process continues until an ordinal tree topology of rank N is reached.
Figure 1:
3.1 A sequential generating process tree topologies
To better illustrate our approach, we begin with some notations. Let Tn (Vn, E.) be tree topology with n leaf nodes and Vn, En are the sets of nodes and edges respectively. Note that [Val 2n and |E,, 2n 3 due to the unrooted and bifurcating structure of Tu. The leaf nodes in V, treated as labeled nodes and the interior nodes in V, are treated as unlabeled nodes. Let us assume a pre-selected order for the leaf nodes X' { (....... &N which is called taxa order for short by us. Now, consider a sequential generating process for all possible tree topologies that have leaf nodes X. We start with definition below.
Definition 1 (Ordinal Tree Topology). Let X (31,...,IN] be set of N(N 3) leaf nodes. Let Tn = (Vnz F,,)) be a tree lopology with n(n S N) leaf nodes in X. We say is an ordinal free topology of rank n, if its leaf nodes are the first n elements of X, i.e. Vn nX (I1....,nn]
We now describe procedure that constructs ordinal tree topologies of rank N recursively by adding one Ieaf node at a time as follows. We first start from T3, the ordinal tree topology of rank 3, which is the smallest ordinal tree topology and is unique due to its unrooted and bifurcating structure. Suppose now we have an ordinal tree topology Tn = (Vn, E,) of rank n. To add the leaf node #n+1 to Tu- we i) select edge e, (u, v) E, and remove from E,; ii) add edges (u, w) (w. v) to the tree topology; iii) add the leaf node Iz4-1 and an edge (w, xn+1) 1) to the tree topology. This way, obtain ordinal tree topology Ta+1 of rank n 1. Intuitively, the leaf node 2x++1 is added to the tree topology Tn by attaching it to an existing edge e,, E En. The position of the selected edge represents the evolutionary relationship between this species and others. After performing this procedure for n 3, .....NN 1, we finally obtain an ordinal tree topology TN of rank N. See Figure for an illustration.
node and two
During the generating process described above, the selected edges at each time step form a sequence D == (e3, ... -x-1). This sequence D of length N 3 records all the decisions have made for autoregressively generating a tree topology 7 and thus we call decision sequence. In fact, there is one-to-one mapping between decision sequences and ordinal tree topologies of rank N, which is formalized in Theorem [ Note that a similar process is also used in online phylogenetic sequential Monte Carlo (OPSMC) (Dinh al. 2016), where the leaf node addition operation is incorporated into the design the proposal distributions.
1} be the set of all decision sequences of length N -3 3 and T be the set of all ordinal free topologies rank N. Let the map
Theorem 1. Let = (D) D = (e3....,eN-- -NN-1), en E En.Y 3 S & N -
g: 5.D
D D
be the generating process described above. Then g is bijection between 'D and T.
According to Theorem [ for each tree topology T E T, there is a unique decision sequence given by g (T). We call this process of finding the decision sequences of tree topologies the decomposition

---

Algorithm 1: ARTree: An autoregressive model for phylogenetic tree topologies
Input: set X = {.'1...., .NN of leaf nodes.
Output: an ordinal tree topology 7 of rank N; the ARTree probability Q(T) of T.
73 (Vs, Eiz) < the unique ordinal tree topology of rank 3; do
for n = ....,
Calculate the probability vector 9n E R Sample an edge decision en. from Discrete (gn) and assume en = (u, v);
Fud using the current GNN model;
Create a new node wi;
En+1 +- -Fn (En((n)) U {(u, w); (w,, v), (w, (n+1))
Vn+1 & Vn U fn451+
Tn+11 f (Vn+1, Ent1);
end
T - TN:
Q(T) gole))aa(es)....p -N-
process. See more details on the decomposition process in Appendix that one can find g (T) in linear time. Lemma 1. The time complexity of the decomposition process induced by g-) is O(N).
C The following lemma shows
The proofs of Theorem and Lemma caan be found in Appendix D Based on the bijection defined in Theorem we model the distribution Q(D) over the space of decision sequences D instead of modeling the distribution (T) over T. Due to the sequential nature of D, we can decompose Q(D) as the product of conditional distributions over the elements:
N-1
Q(D) = II Qle.leg..... Cn-1))
(5)
n=3
In what follows, we simplify (c,, cg,...,, 1) 1) as ((nllenn) and let e<3 be the empty set.
3.2 Graph neural networks for edge decision distribution
By Theorem D the sequence een corresponds to sequence of ordinal tree topologies of increasing ranks ( Tn) (the empty set e<3 corresponds to the unique ordinal tree topology тз of rank 3). Therefore, the discrete distribution Q(en/ezn) in equation (5) defines the probability of adding the leaf node Tn+1 1 to the edge en of Tui conditioned on all the ordinal tree topologies (T3...., Tr) generated so far. In what follows, we will show step by step how to use graph neural networks (GNNs) to parameterize such a conditional distribution given tree topologies.
Topological node embeddings At the n-th time step of the generating process, we first find the node embeddings of the current tree topology on (Vn. En). which is a set f (u) E R E Vn) that assigns each node with an encoding vector in R.N. Following Zhang (023). we first assign one hot encoding to the leaf nodes, i.e.
[x:)) = 8.g, 1 < < n, sj
(6)
where is Kronecker delta function; we then get the embeddings for the interior nodes by minimizing the Dirichlet energy (( Tn) = E(,,)E. (w) (v) using the efficient two-pass algorithm described in Zhang 2023). One should note that the embeddings for interior nodes may change as new leaf nodes are added to the ordinal tree topologies, which is main difference between our model and other graph autoregressive models.
Message passing networks Using these topological node embeddings the initial node features, GNNs apply message passing steps to compute the representation vector of nodes that encode topological information of Tu, where the node features are updated with the information from their neighborhoods in a convolutional manner (Gilmer et al. 2017). More concretely, the 1-th round of message passing is implemented by
mi, (u, u) F Flrcssage (8%, ((),
(7a)
p/ (w) = -Fl Fipptating (fmL (u, u): E N())))
(7b)
5

---

where F message and F ((v) the neighborhood of the node v. In our implementations, the choices of Ett message and Ftt follow the edge convolution operator (Wang al. 2018), while other variants of GNNs can also be applied. The final node features of Ta are given by (J% (v) E Vn} after L rounds of message passing.
updating are the message function and updating function in the 1-th round, and
"updating
Node hidden states The conditional distribution Q((een) is highly complicated as it has to capture how Øn-1 be added to Tn based on how previous leaf nodes added to form the tree topologies. common approach is to use RNNs to model this complex distribution that strikes a good balance between expressiveness and scalability (You et al. 2018b] Liao al. 2019). In our model, after obtaining the final node features of Tu. a gated recurrent unit (GRU) (Cho et al. 2014) follows, i.e.
hn (v) = 1V)
(8)
where h. (v) is the hidden state of u at the n-th generation step and is initialized zero for the newly added nodes including those in 73- The node hidden states (hn (v); v E VaJ. therefore, contain the information of all the tree topologies generated so far which can be used for conditional distribution modeling.
Time guided readout We now construct the distribution Q(-leen) over edge decisions based on the node hidden states. As mentioned before, a main difference between our model and other graph autoregressive models is that the node embedding (v) of node o may vary with the time step We, therefore, incorporate time embeddings into the readout step which first forms the edge features (n(e) E IR of = (u, v) using
(e) = Fppooling (u) + bn, he (v) + bn),
(9a)
(e) Freadout (pn(e) + b,))
(9b)
where b,, is the sinusoidal positional embedding of time step n that is widely used in Transformers (Vaswani et al.) 2017), poooinng is the pooling function implemented elementwise maximum operator, and Freadout is the readout function implemented as 2-layer MLPs with a scalar output. Then the conditional distribution for edge decision
2-layer MLPs followed by an
Q(1lecn) ~ Discrete (ym), In softmax ((ru(e))eEF,,,
(10)
where probability vector yn E IR Ex.l for parametrizing Q(-leen) is obtained by applying softmax function to all the time guided edge features in equation (9b).
As the last step, we sample an edge en, E En from the discrete distribution in equation (10) and add
the leaf node Fn+1 to e,, as described in Section 3.1 This way, we update the ordinal tree topology from Tn of rank n. to Ta+1 of rank n + 1. We can repeat this procedure until an ordinal tree topology T TN of rank N is reached. The probability of r then takes the form
N-1
=
II Qolenleca), 2=3
(11)
where D is the decision sequence and & are the learnable parameters in the model. We call this autoregressive model for tree topologies ARTree, and summarize in Algorithm Note that equation can also be used for tree topology probability evaluation where the decision sequence D obtained from the decomposition process (Appendix C that enjoys a linear time complexity (Lemma Compared to SBNs, ARTree does not rely on heuristic features for parameterization and can provide distributions whose support spans the entire tree topology space as all possible decisions would have nonzero probabilities due to the softmax parameterization in equation (10). Although different taxa orders may affect the performance of ARTree, we find this effect is negligible in our experiments.
There exist other VI approaches for phylogenetic inference that also use unconfined models over the space of tree topologies. Moretti et al. (2021) proposed to sample tree topologies through subtree merging and resampling following CSMC (Wang et al. 2015), but employed a parametrized proposal distribution. Koptagel et al. (2022) proposed a novel multifurcating tree topology sampler named SLANTIS, which makes decisions on adding edges in a specific order based on a simply parameterized weight matrix and maximum spanning trees, and is integrated with CSMC to sample bifurcating tree

---

SBN-EM
ARTree
peak peak
peak peak
astimated proboblity
10-2-
estimated probsbiltt
divergence
20-5
10
20*2-
SAN-EM
SUN LM--LMMa
10a4
10-2
10
10-2
ground
Figure 2: Performances of different methods for TDE on DS1. Left/Middle: Comparison of the ground truth and the estimated probabilities using SBN-EM and ARTree, A tree topology is marked as an outlier if it satisfies log(estimated probability) divergence as a function of the sample size. The results are averaged over 10 replicates with one standard deviation as the error bar.
truth) >> 2. Right: The KI.
topologies (Ø- (B-C.MC). Unlike these methods, ARTree employs GNNs for autoregressive model that builds up the tree topology sequentially through leaf node addition operations, which not only allows fast sampling of trees. but also provides straightforward density estimation procedures. We demonstrate the advantage of ARTree over these baselines the experiments.
4 Experiments
In this section, we test the effectiveness and efficiency of ARTree for phylogenetic inference on two benchmark tasks: tree topology density estimation (TDE) and variational Bayesian phylogenetic inference (VBPI). In all experiments, report the inclusive KL divergence from posterior estimates to the ground truth to measure the approximation error of different methods. We will use "KL divergence" for inclusive KL divergence throughout this section unless otherwise specified. The code available at
Experimental setup We perform experiments on eight data sets which we will call DS1-8. These
data sets, consisting of sequences from 27 64 eukaryote species with 378 to 2520 site observations,
are commonly used to benchmark phylogenets MCMC methods (Hedges et al. 1990; Garey et al.
1996. Yang & Yoder [003 Henk et al. 2003 Lakner et al.] 2008; Zhang & B Blackwell] 2001]] Yoder &
Yang 2004 Rossman et al. 2001 Höhna & Drummond 2012] Larget 2013] W Whidden Matsen IV 2015). For the Bayesian setting,
focus the joint posterior distribution of the tree topologies and the branch lengths and assume a uniform prior on the tree topologies, an i.i.d. exponential prior
Exp(10) on branch lengths, and the simple JC substitution model (Jukes al. 1969). For each of these data sets, we run 10 single-chain MrBayes (Ronquist et al., 2012) for one billion iterations, collect samples every 1000 iterations, and discard the first 25% samples burn-in. These samples form the ground truth of the marginal distribution of tree topologies to which we will compare the posterior estimates obtained by different methods. All GNNs have L rounds in the message passing step. All the activation functions in MLPs are exponential linear units (ELUs) (Clevert et al. 2015). The taxa order is set to the lexicographical order of the corresponding species names in all experiments except the ablation studies. All the experiments are run an Intel Xeon Platinum 9242 processor. All models are implemented in Py PyTorch (Paszke et al.) 2019) and trained with the Adam (Kingma & Ba 2013) optimizer. The learning rate is 0.001 for SBNs, 0.0001 for ARTree, and 0.001 for the branch length model.
4.1 Tree topology density estimation
We first investigate the performance of ARTree for tree topology density estimation given the MCMC posterior samples on DS1-8. Following Zhang & Matsen IV [V] (2018). 8), we run MrBayes on each data set with 10 replicates of chains and runs until the runs have ASDSF (the standard convergence criteria used in MrBayes) less than 0.01 maximum of 100 million iterations. The training data
7

---

Table 1: KL divergences to the ground truth of different methods across 8 benchmark data sets. Sampled trees column shows the numbers of unique tree topologies in the training sets formed by MrBayes runs. The results are averaged over 10 replicates. The results of SBN-EM, SBN-EM-or are from Zhang & Matsen IV (2018).
KL divergence ground truth
Data set
# #Taxa
#Sites
Sampled trees
SBN-EM
SBN-EM-a
SBN-SGA
ARTree
DS1
27
1949
1228
0.0136
0.0130
0.0504
0.0045
DS2
2520
0.0199
0.0128
0.0118
0.0097
DS3
36
1812
43
0.1243
0.0882
0,0922
0.0548
DS4
41
1137
828
0.0763
0.0637
0.0739
0.0299
DS5
378
33752
0.8599
0.8218
0.8044
0.6266
DS6
50
1133
35407
0,3016
0,2786
0,2674
0.2360
DS7
59
1824
1125
0.0483
0.0399
0,0301
0.0191
DS8
1008
3067
0.1415
0.1236
0.1177
0.0741
sets formed by collecting samples every 100 iterations and discarding the first 25%. Now, given a training data set M = {Tm J we train ARTree via maximum likelihood estimation using stochastic gradient ascent. In each iteration, the stochastic gradient is obtained follows
VoL(O: = B
E b=1
V log Q
(12)
B where minibatch {Tm, I randomly sampled from M. We compare ARTree SBN baselines including SBN-EM, SBN-EM-o, and SBN-SGA. For SBN-EM and SBN-EM-a, we use the same setting previously done Zhang Matsen IV 2018) (see Appendix A for more details). In addition to these EM variants, gradient based method for SBNs called SBN-SGA is considered, where SBNs are reparametrized with the latent parameters initialized as zero (see equation (18) in Appendix B and optimized via stochastic gradient ascent, similarly to ARTree. For both ARTree and SBN-SGA, the results are collected after 200000 parameter updates with batch size B
The left and middle plots of Figure 2 show comparison between ARTree and SBN-EM on DS1, which has a peaky posterior distribution. Compared SBN-EM. ARTree provides more accurate probability estimates for tree topologies on the peaks and significantly reduces the large biases in the low probability region (the crimson dots). The right plot of Figure shows the KL divergence of different methods as function of the sample size of the training data. We see that ARTree consistently outperforms SBN based methods for all Ms. Moreover, as the sample size M increases, ARTree keeps providing better approximation while SBNs start to level off when M is large. This indicates the superior flexibility ARTree over SBNs for tree topology density estimation.
Table shows the KL divergences of different methods on DS1-8. We see that ARTree outperforms SBN based methods on all data sets. The gradient based method SBN-SGA is better than SBN-EM on most of the data sets because SBN-EM is well initialized (Zhang & Matsen IV] [V] 2018) and more likely to get trapped local modes. From this point view, the comparison between ARTrec and SBN-SGA is fair because they both use a uniform initialization that facilitates exploration.
4.2 Variational Bayesian phylogenetic inference
Our second experiment on VBPI, where we compare ARTree to SBNs for tree topology variational approximations. Both methods are evaluated the aforementioned benchmark data sets DS1-8 -8. Following Zhang & Matsen IV (2019), we use the simplest SBN and gather the subsplit support from 10 replicates of 10000 ultrafast maximum likelihood bootstrap trees (Minh et al. 2013). For both ARTree and SBNs, the collaborative branch lengths are parametrized using the learnable topological features with the edge convolution operator (EDGE) for GNNs (Zhang) 2023). We set K 10 for the multi-sample lower bound 4 and use the following annealed unnormalized posterior at the i-th iteration
o(,YT,g; p( Y ((.
(13)
8

---

30
10"4
SBN
-200
veriational ampoximatoon
20
KL divergence
ARtree
72360
20-1
36a
1000
ARTree(10)
ARThee(S)
S8N(10) SBN(5)
10-1
-3200
ARTroc"(10)
101
10
10
10"2
0
50
Iterations
ground truth
order index
Figure 3: Performances of ARTree and SBN tree topology variational approximations for VBPI on DSI. Left: the evidence lower bound (ELBO) as function of iterations. The numbers of particles used in the training objective in the brackets. The ARTree" method refers to ARTree without time guidance, i.e. ln2 for all n in the readout step. Middle: variational approximations vs ground truth posterior probabilities of the tree topologies. Right: KL divergences across 50 random taxa orders. The KL divergence of SBNs is averaged over 10 independent trainings.
Table KI. divergences to the ground truth, evidence lower bound (ELBO), 10-sample lower bound (LB-10), and marginal likelihood (MI.) estimates of different methods across 8 benchmark data sets. GT trees row shows the number of unique tree topologies in the ground truth. The marginal likelihood estimates are obtained via importance sampling using 1000 samples. The KL results are averaged over 10 independent trainings. For ELBO, LB-10, and ML, the results are averaged over 100, 100, and 1000 independent of &CCSMC are from Koptagel et al. (2022). For ELBO and LB-10, a larger mean better; for ML, a smaller standard deviation is better
respectively with standard deviation in the brackets. The results
Data
DSI
DS2
DS3
DS4
DSS
DS6
DS7
DSS
Taxa
27
29
36
41
50
50
59
64
Sites
1949
2520
1812
1137
378
1133
1824
1008
GT trees
2784
42
351
11505
1516877
809765
11525
82162
SBN
0.0707
0.0141
0.0554
0.0739
1,2472
0.3795
0.1531
0.3173
ARTree
0.0097
0.0004
0.0064
0.0219
0.8979
0.2216
0.0123
0.1231
LB-10 ELBO
SBN
-7100.200005 58/0.03
-33736.22/0.02)
-13331.83(0.03)
82177.8000001
-6728.65(0.06)
-37334.85(000)
-8655.05(0.05 (0.0505
-8655
AKTiee
-7110.990.04)
-263687710.077
133.36.1700.88)
-13331.82((005)
-8217.68/0004)
-6728.65(0.66
-37.334.8400133)
-665..0(0005))
SBN -710869X0.02)
-26367.5700022
-33755.26/0002
-13330.2910.00
-8215-2((0.01)
-6725.33(0.04)
-37332.58(0 58(0.03)
-8651.78(0.D)
ARTree -7108.68/0002)
-26367.8600822)
(53.33.25(0.02)
(13.30.27(0.33)
-8215.34(0.83)
-6725.33(0.44)
-37.332554(0033)
-6651.730.04)
d-CSMC 239003667.233
-30568 30568-49(31..M)
-33798.06// 99..866666666 62)
-13582244355006
-8367 (16.7.1(8887))
-7013.83(16.99)
N/A
-9209.18111 18/18.8
SBN -7108.410.155)
26367.71(00088
.33738.0900099)
-13329 -1352999400.00
-8214.6200.000
-6724.370..35)
-37331 97(0.28)
-8650.64400.5
ARTree
-1008.110.19)
266367.71(0.07)
-33735.0N0.09)
-13329.94(0.17)
-8214.59(0.34)
-6724.3700.66)
-37331.9(0027)
-8650.61(0.48)
where 3, min 1.0, 0.001 + i/H) is the inverse temperature that goes from 0.001 to 1 after H iterations. For ARTree, long annealing period H 200000 is used for DS6 and DS7 due to the highly multimodal posterior distributions on these two data sets (Whidden & Matsen IV 2015) and H 100000 is used for the other data sets. For SBNs, we set // 100000 for all data sets. The Monte Carlo gradient estimates for the tree topology parameters and the branch lengths parameters are obtained via VIMCO (Mnih & Rezende 2016) and the reparametrization trick (Zhang & Matsen IV 2019) respectively. The results are collected after 400000 parameter updates.
The left plot in Figure shows the evidence lower bound (ELBO) as function of the number of iterations on DS1. Although the larger support of ARTree adds to the complexity of training for tree topology variational approximation, we see that by the time SBN based methods converge, ARTree based methods achieve comparable (if not better) lower bounds and finally surpass the SBN baselines in the end. We also find that using fewer particles (K 5) in the training objective tends
We use standard deviation a criterion because in context importance sampling, the variance of an estimator reflects the approximation accuracy of the importance distribution to the target. Note that this is only reasonable if the importance distribution is proper (so that the marginal likelihood is unbiased), and hence does apply to G-CSMC even though it also has much larger variance.

---

to provide larger ELBO. Moreover, time guidance turns out to be crucial for ARTree, as evidenced by the significant performance drop when it is turned off. As shown in the middle plot, compared to SBNs, ARTree can provide a more accurate variational approximation of the tree topology posterior. To investigate the effect of taxa orders on ARTree, we randomly sample 50 taxa orders and report the KL divergence for each order in the right plot of Figure We find that ARTree exhibits weak randomness as the taxa order varies and consistently outperforms SBNs by a large margin.
Table 2 shows the KL divergences to the ground truth, evidence lower bound (ELBO), 10-sample lower bound (LB-100) and marginal likelihood (ML) estimates obtained by different methods on DS1-8. We find that ARTree achieves smaller KL divergences than SBNs across all data sets and performs on par or better than SBNs for lower bound and marginal likelihood estimation. Compared to SBNs, the ELBOs provided by ARTree tend to have larger variances, especially on DS2, DS3, and DS7, which is partly due to the larger support of ARTree that spans the entire tree topology space (see more discussions in Appendix
5 Conclusion
In this paper, we introduced ARTree, deep autoregressive model over tree topologies for phylogenetic inference. Unlike SBNs that rely on hand-engineered features for parameterization and require pre- sampled tree topologies, ARTree is built solely top of learnable topological features (Zhang) 2023) via GNNs which allows for a rich family of distributions over the entire phylogenetic tree topology space. Moreover, procedures and straightforward density computation, which make readily usable for tree topology density estimation and variational Bayesian phylogenetic inference. In experiments, ARTree outperforms SBNs on a benchmark of challenging real data tree topology density estimation and variational Bayesian phylogenetic inference problems, especially in terms of tree topology posterior approximation accuracy.
autoregressive model, ARTree also allows simple forward sampling
showed that
Acknowledgements
This work was supported by National Natural Science Foundation of China (grant no. 12201014 and grant no. 12292983), as well as National Institutes of Health grant AI162611. The research of Cheng Zhang was support in part by National Engineering Laboratory for Big Data Analysis and Applications, the Key Laboratory of Mathematics and Its Applications (LMAM) and the Key Laboratory of Mathematical Economics and Quantitative Finance (LMEQF) of Peking University. The authors are grateful for the computational resources provided by the High-performance Computing Platform of Peking University. The authors appreciate the anonymous NeurlPS reviewers for their constructive feedback.
References
Stephen W. Attwood, Sarah C. Hill, David M. Aanensen, Thomas R. Connor, and Oliver G. Pybus. Phylogenetic and phylodynamic approaches to understanding and combating the early SARS-CoV. 2 pandemic. Nature Reviews. Genetics, 23:547 562, 2022.
Nicola De Cao and Thomas Kipf. MolGAN: An implicit generative model for small molecular graphs, 2018. /tt.sttss///ariiiii org. httpr//ww.orr/a///8800..... wwwpaa/sss......... 11973
Kyunghyun Cho, Bart Van Merriënboer, Caglar Gulcehre, Dzmitry Bahdanau, Fethi Bougares, Holger Schwenk, and Yoshua Bengio. Learning phrase representations using RNN encoder-decoder for statistical machine translation. arXiv preprint urXiv: /406. 1078, 2014.
Djork-Aréé Clevert, Thomas Unterthiner, and Sepp Hochreiter. Fast and accurate deep network learning by exponential linear units (ELUs). arXiv: Learning, 2015.
Hanjun Dai, Azade Nazi, Yujia Li, Bo Dai, and Dale Schuurmans. Scalable deep generative modeling for sparse graphs. In International Conference on Machine Learning, pp. 2302-2312, 2020.
Rob DeSalle and George Amato. The expansion of conservation genetics. Nut. Rev. Genet,, 5(9): 702-712, September 2004, ISSN 1471 1471-0056. doi: 10. 10.1038//rg1433 1425. URL http:/// doi org/ 10. 1038/ 1038//nrg//2/ nrg4425
10

---

Vu Dinh, Arman Bilge, Cheng Zhang, and Frederick A Matsen IV. Probabilistic path Hamiltonian Monte Carlo. In Proceedings of the 34th International Conference on Machine Learning, pp. 1009-1018, 8, July 2017. URL http::ttpocceeeeiiiiiiie ppproeeddgg...r.pprrsss/////ail hattmmll....m
mlr ptress/v70/dinh//2
Vu C. Dinh, Aaron E. Darling, and Frederick A. Matsen IV. Online Bayesian phylogenetic inference: Theoretical foundations via sequential Monte Carlo. Systematic Biology, 67:503 217,2016.
Louis du Plessis, John T McCrone, Alexander E Zarebski, Verity Hill, Christopher Ruis, Bernardo Gutierrez, Jayna Raghwani, Jordan Ashworth, Rachel Colquhoun, Thomas Connor, Nuno R Faria, Ben Jackson, Nicholas J Loman, Áine O'Toole, Samuel M Nicholls, Kris Parag, Emily Scher, Tetyana I Vasylyeva. Erik M Volz, Alexander Watts, Isaac I Bogoch, Kamran Khan, COVID-19 Genomics UK (COG-UK) Consortiumt' David M Aanensen, Moritz G Kraemer, Andrew Rambaut, and Oliver G Pybus. Establishment and lineage dynamics of the SARS-CoV epidemic in the UK. Science, January 2021. ISSN 0036-8075, 1095-9203. doi: 10. 126/ science.abf2946. URL science.abf2944
https: htops://sccisn /ttpsst//ccinnnns http://ccencceceec sciencemag, http////eeee//aa//22///////////0
Gytis Dudas, Luiz Max Carvalho, Trevor Bedford, Andrew Tatem, Guy Baele, Nuno R Faria, Daniel J Park, Jason T Ladner, Armando Arias, Danny Asogun, Filip Bielejec, Sarah L Caddy, Matthew Cotten, Jonathan D' Ambrozio, Simon Dellicour, Antonino Di Caro, Joseph W Diclaro, Sophie Duraffour, Michael J Elmore, Lawrence Fakoli, Ousmane Faye, Merle L Gilbert, Sahr M Gevao, Stephen Gire, Adrianne Gladden- Young, Andreas Gnirke, Augustine Goba, Donald S Grant, Bart L Haagmans, Julian A Hiscox, Umaru Jah, Jeffrey Kugelman, Di Liu, Jia Lu, Christine M Malboeuf, Suzanne Mate, David A Matthews, Christian B Matranga, Luke W Meredith, James Qu, Joshua Quick, Suzan D Pas, My V T Phan, Georgios Pollakis, Chantal B Reusken, Mariano Sanchez-Lockhart, Stephen F Schaffner, John S Schieffelin. Rachel S Sealfon, Etienne Simon- Loriere, Saskia L. Smits, Kilian Stoecker, Lucy Thorne, Ekaete Alice Tobin, Mohamed Vandi, Simon J Watson, Kendra West, Shannon Whitmer, Michael R Wiley, Sarah M Winnicki, Shirlee Wohl, Roman Wölfel, Nathan Yozwiak, Kristian Andersen, Sylvia Blyden, Fatorma Bolay, Miles W Carroll, Bernice Dahn, Boubacar Diallo. Pierre Formenty, Christophe Fraser, George F Gao, Robert Garry, Ian Goodfellow, Stephan Günther, Christian T Happi, Edward C Holmes, Brima Kargbo, Sakoba Keita, Paul Kellam, Marion P G Koopmans, Jens H Kuhn, Nicholas J Loman, N' f'laly Magassouba, Dhamari Naidoo, Stuart T Nichol, Tolbert Nyenswah. Gustavo Palacios, Oliver G Pybus, Pardis C Sabeti, Amadou Sall, Ute Ströher, Isatta Wurie. Marc A Suchard, Philippe Lemey, and Andrew Rambaut. Virus genomes reveal factors that spread and sustained the Ebola epidemic. Nature, April 2017. ISSN 0028-0836, 1476-4687. doi: 10. 10000//auurre22400. URL httt:://xx.ddo orggrrrgg/o 10. 1038/ nat/un2uuur22200
J. Felsenstein. Evolutionary trees from DNA sequences: A maximum likelihood approach. Journal of Molecular Evolution, 17:268-276, 1981.
Joseph Felsenstein. Inferring Phylogenies. Sinauer associates, edition, 2004.
J. R. Garey, T. Near, M. R. Nonnemacher, and S. A. Nadler. Molecular evidence for Acanthocephala as subtaxon of Rotifera. Mol. Evol., 43:287-292, 1996.
Justin Gilmer, Samuel Schoenholz, Patrick F, Riley, Oriol Vinyals, and George E, Dahl. Neural message passing for quantum chemistry, ArXiv, abs/1704.01212, 2017.
S. B. Hedges, K. D. Moberg, and L. Maxson. Tetrapod phylogeny inferred from 18S and 28S ribosomal RNA sequences and review of the evidence for amniote relationships. Mol. Biol. Evol., 7:607-633, 1990.
D. A. Henk, A. Weir, and M. Blackwell. Laboulbeniopsis termitarius, an ectoparasite of termites newly recognized member of the Laboulbeniomycetes. Mycologia, 95:561-564, 2003.
Sebastian Höhna and Alexei Drummond. Guided tree topology proposals for Bayesian phylogenetic inference. Syst. Biol 6 111): January 2012. ISSN 1063-5157. doi: 10.1093/sysbioosyr777. URL httt:://ttt..eee doi org/ 10. /000//////o/////////////
Huelsenbeck, Ronquist, R. Nielsen, and Bollback. Bayesian inference phylogeny and its impact evolutionary biology. Science, 294:23 10-2314, 2001.

---

Wengong Jin, Regina Barzilay, and T. Jaakkola. Junction tree variational autoencoder for molecular graph generation. ArXiv, abs/1802.04364, 2018.
Thomas H Jukes, Charles R Cantor, et al. Evolution of protein molecules. Mammalian protein metabolism, 3:21-132, 1969.
D. Kingma and J. Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
Hazal Koptagel, Oskar Kviman, Harald Melin, Negar Safinianaini, and Jens Lagergren. VaiPhy: variational inference based algorithm for phylogeny. In Advances Neural Information Processing Systems, 2022.
C. Lakner, van der Mark, J. P. Huelsenbeck, B. Larget, and F. Ronquist. Efficiency of Markov chain Monte Carlo tree proposals in Bayesian phylogenetics. Syst. Biol., 57:86-103, 2008.
Bret Larget. The estimation of tree posterior probabilities using conditional clade probability distributions. Syst. Biol., 5014):50 July 2013. ISSN 1063-5157. doi: 00.109//ssbboo/st0014
URL sattp://dx wttp/xxxxxiiii..//ooooooo//..d org/ 10. wtt.://9//yyooooyyo/o0//
Bret R. Larget and L.. Simon. Markov chain Monte Carlo algorithms for the Bayesian analysis of phylogenetic trees. Molecular Biology and Evolution, 16:750-750, 1999.
Y. Li. O. Vinyals, C. Dyer, Pascanu, and Battaglia. Learning deep generative models of graphs. 2018. URL htttps///axxi> org/ abs/ 180108883.0011 03324
Renjie Liao, Yujia Li, Yang Song, Shenlong Wang, Charlie Nash, William L. Hamilton, David Krist- janson Duvenaud, Raquel Urtasun, and Richard S. Zemel. Efficient graph generation with graph recurrent attention networks. ArXiv, abs/19100007a 2019.
B. Mau, M. Newton, and B. Larget. Bayesian phylogenetic inference via Markov chain Monte Carlo methods. Biometrics, 55:1-12, 1999.
Bui Quang Minh, Minh Anh Nguyen, and Amdt von Haeseler. Ultrafast approximation for phyloge- netic bootstrap. Molecular Biology and Evolution, 30:1 188 1195, 2013.
Andriy Mnih and Danilo Jimenez Rezende. Variational inference for monte carlo objectives. In International Conference on Machine Learning, 2016.
Antonio Khalil Moretti, Liyi Zhang, Christian Andersson Naesseth, Hadiah Venner, David M. Blei, and Itsik Pe'er er. Variational combinatorial sequential Monte Carlo methods for Bayesian phylogenetic inference. In Conference Uncertainty in Artificial Intelligence, 2021.
Adam Paszke, Sam Gross, Francisco Massa, Adam Lerer, James Bradbury, Gregory Chanan, Trevor Killeen, Zeming Lin, Natalia Gimelshein, Luca Antiga, Alban Desmaison, Andreas Köpf, Edward Yang, Zach Detito, Martin Raison, Alykhan Tejani, Sasank Chilamkurthy, Benoit Steiner, Lu Fang, Junjie Bai, and Soumith Chintala. PyTorch: An imperative style, nigh-performance deep learning library. In Neural Information Processing Systems, 2019.
Tom Rainforth, Adam R. Kosioreck, Tuan Anh Le, Chris J. Maddison, Maximilian Igl, Frank Wood, and Yee whye Teh. Tighter variational bounds are not necessarily better. In Proceedings of the 36th International Conference on Machine Learning, 2019.
Fredrik Ronquist, Maxim Teslenko, Paul Van Der Mark, Daniel L. Ayres, Aaron Darling, Sebastian Höhna, Bret Larget, Liang Liu, Marc Suchard, and John Huelsenheck. MrBayes 3.2: Efficient Bayesian phylogenetic inference and model choice across a large model space. Systematic Biology. 61(3):539-542, 2012.
A. Y. Rossman, J. M. Mckemy, R. A. Hardo-Schlttheiss and H. J. Schroers. Molecular studies of the Bionectriaceae using large subunit rDNA sequences. Mycologia, 93:100-110, 2001.
Chence Shi, Minkai Xu, Zhaocheng Zhu, Weinan Zhang, Ming Zhang, and Jian Tang. GraphAF: flow-based autoregressive model for molecular graph generation. In International Conference on Learning Representations, 2020. URL htt.tp/se//pppeeeeeiieeeeetootoo
12

---

Martin Simonovsky and Nikos Komodakis. GraphVAE: Towards generation of small graphs using variational autoencoders, In International Conference on Artificial Neural Networks, pp. 412-422, 2018.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Advances in Neural Information Processing Systems, volume 30, 2017.
Liangliang Wang, Alexandre Bouchard-Côté, and A. Doucet. Bayesian phylogenetic inference using a combinatorial sequential Monte Carlo method. Journal of the American Statistical Association, 110:1362- 1374, 2015. httpss//:// www.stieacccoooooooooooo haattpp:s/ppi.seeeattcccccooool
org/ :::::/ooo:::::::::::: 44955399
Yue Wang, Yongbin Sun, Ziwei Liu, Sanjay E. Sarma, Michael M. Bronstein, and Justin M. Solomon. Dynamic graph CNN for learning on point clouds. ACM Transactions on Graphics (TOG), 38:1 12, 2018.
Chris Whidden and Frederick Matsen IV. Quantifying MCMC exploration of phylogenetic tree space. Syst. Biol., 64(3)442--491, May 2015, ISSN 1063-5157, 1076-836X. doi: 10. s0.1093/sysbio/ syv006. URL http::/dx.....
doi wtt.://......oooo//ooo/oo//////e 1wte110///////yyy/eer
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., 52:705-716, 2003.
Ziheng Yang and Bruce Rannala. Bayesian phylogenetic inference using DNA sequences: a Markov chain Monte Carlo method. Molecular Biology and Evolution, 14(7):717-724, 1997.
A. D. Yoder and Yang. Divergence datas for Malagasy lemurs estimated from multiple gene loci: geological and evolutionary context. Mol. Ecol., 13:757-773, 2004.
Jiaxuan You, Bowen Liu, Zhitao Ying, Vijay Pande, and Jure Leskovec. Graph convolutional policy network for goal-directed molecular graph generation. In Advances in Neural Information Processing Systems, pp. 6410-64 (000-621,2220 2018a.
Jiaxuan You, Rex Ying, Xiang Ren, William L. Hamilton, Jure Leskovec. GraphRNN: Generating realistic graphs with deep auto-regressive models. In International Conference on Machine Learning, 2018b.
Cheng Zhang. Improved variational Bayesian phylogenetic inference with normalizing flows. In Neural Information Processing Systems, 2020.
Cheng Zhang. Learnable topological features for phylogenetic inference via graph neural networks. In International Conference on Learning Representations, 2023.
Cheng Zhang and Frederick A Matsen IV. Generalizing tree probability estimation via Bayesian networks. In Neural Information Processing Systems, 2018.
Cheng Zhang and Frederick A Matsen IV. Variational Bayesian phylogenetic inference. In Interna- tional Conference on Learning Representatiinss 2019.
Cheng Zhang and Frederick A Matsen IV, A variational approach to Bayesian phylogenetic inference, 2022. URL https://arxii wwtp:/iiiivrgg/sbbs.eeee wttg/bbb/2204.07////3
N. Zhang and M. Blackwell. Molecular phylogeny of dogwood anthracnose fungus (Discula destruc- tiva) and the Diaporthales. Mycologia, 3655 2001.
13

---

D
Figure 4: Subsplit Bayesian networks and a simple example for a leaf set of 4 taxa (denoted by A, B, C respectively). Left: General subsplit Bayesian networks. The solid full and complete binary tree network is Bx.- The dashed arrows represent the additional dependence for more expressiveness. Middle Left: Examples of (rooted) phylogenetic trees that are hypothesized to model the evolutionary history of the taxa. Middle Right: The corresponding subsplit assignments for the trees. For case of illustration, subsplit , Z) is represented as y in the graph. Right: The SBN for this example, which is By in this case.
Z
A Details of subsplit Bayesian networks
One recent and expressive graphical model that provides flexible family of distributions over tree topologies is the subsplit Bayesian network, as proposed by Zhang & Matsen IV (2018). Let X' be the set of N labeled leaf nodes. A non-empty set C of X is referred as clade and the set of all clades of X, denoted by C(X²), is a totally ordered set with a partial order > (e.g.. lexicographical order) defined on it. An ordered pair of clades (W, Z) is called a subsplit of clade C if it is bipartition of C, i.e. W > Z, W n Z = 0, and W U U = C.
Definition 2 (Subsplit Bayesian Network). A subsplit Bayesian network (SBN) Bx leaf node set X of size N is defined as a Bayesian network whose nodes take on subsplit or singleton clade values of X and has the following properties: (a) The root node of Bx takes subsplits of the entire labeled leaf node set X; (b) Bx contains a full and complete binary tree network By as a subnetwork; (c) The depth of Bx is N - 1, with the root counted depth 1.
Due to the binary structure of Bx. the nodes in SBNs can be indexed by denoting the root node with S and two children of S; with S2; and S2;+1 recursively where S, is an internal node (see the left plot Figure 4.. For any rooted tree topology, by assigning the corresponding subsplits or singleton clades values [S: = 8:):21 to its nodes, one can uniquely map it into an SBN node assignment (see the middle and right plots Figure
As Bayesian networks, the SBN based probability of rooted tree topology takes the following form
Psba (T' = = p(S = s1) II p(S: = ss = =z.):
(14)
1>1
where T. is the index set of the parents of node For unrooted tree topologies, we can also define their SBN based probabilities by viewing them as rooted tree topologies with unobserved roots and integrating out the positions of the root node as follows:
Pebn (T r)
Σ eEF(T)
Pabn (T9)
(15)
where 7" is the resulting rooted tree topology when the rooting position is on edge
In practice, SBNs are parameterized according to the conditional probability sharing principle where the conditional probability for parent-child subsplit pairs are shared across the SBN network, regardless of their locations. The set of all conditional probabilities are called conditional probability tables (CPTs). Parameterizing SBNs, therefore, often requires finding an appropriate support of CPTs. For tree topology density estimation, this can be done using the sample of tree topologies that is given as the data set. For variational Bayesian phylogenetic inference, as no sample of tree topologies is available, one often resorts to fast bootstrap or MCMC methods (Minh et al. [0133] Zhang] 2020). Let S, denotes the root subsplits and SchIpa denotes the child-paren subsplit pairs in the support. The parameters of SBNs are then p = {par $1 E S, U {Pape s E Sch|pa) where
Ps1 p(S 51),
Pajc
p(S, s)Sx: t), vi > 1.
(16)
14

---

As result, the supports of SBN-induced distributions are often limited by the splitting patterns in the observed samples and could not span the entire tree topology space (Zhang & Matsen IV 2022).
The SBN-EM Algorithm For unrooted tree topologies, the SBN based probability 15 can be viewed as hidden variable model where the root subsplit is the hidden variable. In this case, expecation-maxxiizattin maximization (EM) algorithm,
SBNs can be trained using the expectation
proposed by Zhang & Matsen IV (218)). Given training set [Tk we first initialize the parameter estimates as ( (i.c., the simple average estimates as in Zhang & Matsen IV (2018).
In the i-th step, we run
the E-step and M-step as follows
P(TE P(TENNLIEEM.. (EM
E-step: L1 S & S M, compute d (s1) Enjory P (p where 81 ~ ты means
the subsplit 81 can be achieved by placing a "virtual root" on an edge of
M-step: update the parameter estimates by
RM.((41))
Mai + +ams
Σ
= = (s))
K + a E.jes, ma
m
k=
BMC(+11 Paje
ma.t am,.
m
Σ
Σ
dh (si,x)
Σ 1((0.11 (5/ =,
Σ. m + ama.))
1 B(T()
where I is the indicator function, 5i,8 is the node value of Sj in The m, and m,, equivalent counts and is the regularization coefficient that encourages generalization.
When a > 0, this algorithm is called SBN-EM-a.
B Details of variational Bayesian phylogenetic inference
With two variational families (T( and Qu(qlT) over the space of tree topologies and branch lengths, the variational Bayesian phylogenetic inference (VBPI) approach forces @ q) = ((( to approximate the posterior p(7, Y by maximizing the following multi-sample lower bound
1
p(Y
"''r'g'')p''' Qu(
L (t. 4 E
log
Σ
(17)
Q..
K
Gradients of the objective 17 w.r.t. & and ab can be estimated by the VIMCO estimator (Mnih & Rezende (016) and the reparameterization trick respectively. the following, we introduce some common choices of (T) and Q (glT).
Choice (T) Before the proposed ARTree framework in this article, SBNs is the common choice of (T). As introduced in Appendix A SBNs provide a probability distribution over unrooted tree topologies in equation (15.. Given a subsplit support of CPTs, SBNs can be parameterized as follows
exp(o,,)
exp( 1)
Psi Eves exp(0x)) 81 E S,; Pole Evia exp(o.)1)) sit
Schlpa-
(18)
The parameters Φ (0r; $1 e (dsee s/t Sch|pa) are called latent parameters SBNs.
Choice of (glT) The distribution Qu ( is often taken to be diagonal lognormal distribution, which can be parametrized using some heuristic features (Zhang & Matsen IV 2019) or the recently proposed learnable topological features (Zhang) 2023) of as follows. For each edge (u, v) in t, one can first obtain the edge features using bhp = [hh,, hw) where h., is the GNN output at node and f permutation invariant function. Then the mean and standard deviation parameters given by
p(e, -) (he)"(he), a(e,7) MLP" (he)
where MLP and MLP are two multi-layer perceptrons (MLPs).
15

---

Algorithm 2: Tree topology decomposition process
Input: tree topology = with all of the N leaf nodes.
Output: a decision sequence D.
TN (Vx, EN) & the topology T;
for n. = N - 1..... 3 do
Determine the unique neighbor of the leaf node In+1i Determine the two neighbors u and (except Tx+1) of w;
Vn Vntilfur Vnt,\fu, xnt1);
En f (Ent1U) Ufu, ((,,)) Int1), (w, u), (w. u)};
In f (Vn, En);
en t u);
end
D (E3,...., eN -N-
C Details of tree topology decomposition process
The tree topology decomposition process, which maps tree topology to corresponding decision sequence, is indeed the inverse operation of Algorithm [ Also, the decomposition process is implemented in recursive way starting from the tree topology тv of rank N. Intuitively, given ordinal tree topology Ta+1 rank n + 1, one can detach the leaf node In+1 as well as its unique neighbor and reconnect the two neighbors of The remaining graph, denoted Tn. is an ordinal tree topology of rank n; the edge decision en is exactly the reconnected edge. This process continues until the unique ordinal tree topology 7a of rank 3 is reached. We summarize the sketch of tree topology decomposition process in Algorithm 2
Given tree topology with all of the N leaf nodes, by first mapping it to decision sequence D (E3,.... EN-1) following Algorithm 2 and then calculate the probability the product of conditionals
evaluate its ARTree based probability
N-1
Q(T) -
II Q(enleen).
D The proofs of Theorem [nd Lemma
1
Theorem 1
1 be the set of all decision sequences of length N -3 3 and T be the set of all ordinal tree topologies rank N. Let the map
Let D (D D (e3...., CN-1), en E
9: D D
D
D
be the generating process used in ARTree. Then g is a bijection between 'D and
Proof of Theorem It is obvious that is a well-defined map from D to T. To prove it is bijection, it suffices show g is injective and surjective.
We first show that is injective. Assume there are two decision sequences D D E D and
D(1 # D(2), Let & be the first position where (1)) and D begin differ, i.c. e!") e (2) for all i < k and ek
one can take the subtree topology of rank k and l: ++1 of T, тæ and Tk+1- Noting that e₁ refers the edge where to add the new node #n+1s (1)
(1) * (2) If g((()) = ((((2))
the equation E * ck implies they will induce different Ts+1s. This contradicts the uniqueness of
Th+1. Therefore, we conclude that g is injective.
Next, we prove that is surjective. For a tree topology 7 with rank N, we denote its subtree topology of rank & by Tke where & 3. n. In fact, for each k. the topology T corresponds to adding the leaf node Ik+1 to an edge in The and we denote this edge by ex. It is easy to verify that the constructed (e3,.. (ea..... . - 1) is preimage of
16

---

Table 3: The ELBO estimates on DS1 obtained by different combinations of tree topology model Q(T) and branch length model Q(q(T). The results are averaged over 100 independent runs with standard deviation in the brackets.
Model combination
ELBO
SBN + branch length model trained along with SBN SBN + branch length model trained along with ARTree -7110.26(0.03)
-7110.24(0.03)
ARTree + branch length model trained along with ARTree -7110.09(0.04)
Table 4: Runtime comparison in the variational inference setting on DS SBN and ARTree refer to the early stopping of SBN and ARTree that surpass the CSMC baseline in terms of marginal likelihood estimation (-7290.36), respectively. The experiments are run on single core of MacBook Pro 2019.
VCSMC
VaiPhy
@-CSMC
SBN
ARTree
SBN'
ARTree
Total training time (minutes)
248.3
45.1
N/A
659.3
3740.8
10.2
79.5
Evaluation time
2.4
1.6
102.2
0.15
0.41
0.15
0.41
(one estimate of ML, minutes)
Lemma 1 The time complexity of the decomposition process induced by (-) is O(N).
Proof of Lemma I Assume the tree topology T is stored binary tree data structure, where each node other than the root node also has parent node pointer. Before decomposing т, we first build dictionary of the xx) mappings by traversal across that maps n to the leaf node Hn in т, vn S N. The time complexity of building this dictionary is O(N). In the decomposition process, given an ordinal tree topology of rank n + costs O(1) time to locate xn+1 1 and determine the unique parent node w of Tn+1 and the neighbors of w. It is obvious that the time complexity of detaching and reconnecting operations is O(1). One can repeat this procedure for = N ,,...,3. resulting in time complexity of ((N). Therefore, the time complexity of the decomposition process induced by g (-) is O(N).
E
Limitations
Larger variances of ELBO estimates As an autoregressive model for phylogenetic tree topology. ARTree provides reliable approximations for target distributions over tree topologies, as evidenced by our real data experiments. However, when incorporated with the branch length model for VBPI, the ELBOs provided by ARTree tend to have larger variances, which find is caused by the occasional occurrence of "outliers" among samples. In fact, as the support of ARTree spans over the entire tree topology space, it adds to the difficulty of fitting the conditional distribution (g(())) compared to SBNs. When combined with ARTree, the approximation accuracy of Q (g(7) might be related to the cluster structure of peaks in the tree topology posterior. See Figure 3 in Whidden & Matsen IV (2015) for cluster subtree-prune -ruue-andd-egrafff (SPR) graphs of DS1-8. We also examine the approximation accuracies of (9/7) trained along with ARTree for those in the support of SBNs and find significant enhancement in ELBO and reduction of variances. This phenomenon raises an important topic: proper design and optimization of branch length model when the support of tree topology model spans the entire space. We leave this for future work.
"Training: We trained all models following the setting in their original papers: VCSMC was trained for 100 iterations with 2048 particles per iteration; VaiPhy was trained for 200 iterations with 128 particles per iteration; G-CSMC directly estimates ML based on VaiPhy, and therefore does not need extra training; both ARTree and SBN were trained for 400000 iterations with 10 particles per iteration.
Evaluation: We find that the evaluation strategies in their original papers quite different, e.g. VaiPhy, SBN, ARTree used importance sampling to estimate ML with different repetition times; VCSMC and 6CCSMC instead estimated ML with sequential Monte Carlo (SMC), also with different repetition times. To be fair, we report the time for producing one estimate of ML from each of these models (VaiPhy. SBN, and ARTree used importance sampling with 1000 particles; VCSMC and 6-CSMC used SMC with 2048 particles).
17

---

Minor increasement of lower bounds From Table the estimates of lower bounds (as well as marginal likelihoods) of ARTree do not improve much over those of SBNs. We provide two explanations for this: (i) The lower bounds in VBPI are more sensitive to the quality of branch length model Q(gIT) instead the tree topology model Q(T). As SBN and ARTree use the same parametrization of the branch length model, we do not expect a large improvement in lower bounds. (ii) The support of ARTree spans the entire tree topology space, This adds to the difficulty training ((gTT) which is conditioned on tree topology as discussed in the above paragraph. Finally, from T 3
we see that the enhancement on ELBOs of ARTree indeed comes from a better tree topology model as evidenced by the result of the 'SBN branch length model trained along with ARTree' combination.
Time complexity Another limitation of ARTree is its time complexity. We compare the training time and evaluation time in Table Both ARTree and SBNs take more time in training than other methods with unconfined support since they both build machine-learning models with enormous parameters and rely heavily on optimization. However, they take a comparable amount of time to provide good enough approximations for marginal likelihood estimation of similar accuracy, as evidenced by ARTree" and SBN ARTree takes more time than SBNs because it relies on several submodules which, although complicated, accommodate the complex tree space and are widely-used strategies in the literature. The inefficiency of autoregressive generative models also an inherent issue.
designed to promote the expressive power
18


Reviews:
[{"review_id": "JFRdApSJhz", "reviewer": "Reviewer_tXZx", "summary": "This paper presents a new way to construct the variational distribution of tree topologies based on autoregressive generation with GNNs, which is used in the problem of variational phylogenetic inference. The paper mainly compares the new construction of the variational distribution with subsplit Bayesian network approaches. The experiments show that the proposed approach outperforms SBN in terms of learning the ground-truth tree topologies.", "questions": "Can the authors talk about the potential use of the proposed method in computational biology in addition to getting better EBLO or other metrics of modelling the data? such as interpretability.", "limitations": "N/a", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "1. The paper seems to address an interesting sub-problem of the task, which is how to flexibly generate tree topologies for the variational distribution. I'm not a domain expert in variational phylogenetic inference, but previous approaches usually follow the way of SBNs and this paper proposes a novel and better alternative to SBNs.\n\n2. The proposed approach of autoregressive generation with GNNs looks intuitive. The paper provides comprehensive experiments in the comparison with SBN, which support the claims of the proposed method.\n\n3. The paper is well-written and easy to follow.", "weaknesses": "1. As stated in the paper, the main drawback of SBN is that it could not span the entire tree topology space. It seems that there is no analysis on how/why the proposed method is better at doing this in addition to the empirical comparison in the experiments.\n\n2. The main contribution of the paper is in Section 3.1, which is a new parameterisation of $Q(r)$. Most of the techniques in Section 3.2 follows Zhang (2023). It might be better to shorten 3.2, as I think it might not the focus of the paper.\n\n3. Although Table 1 shows that ARTree has better number in terms of revealing the ground-truth trees, it seems that ARTree does not improve much on ELBO and ML in Table 2.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable questions! We address your comments and questions as below.\n\n**Weakness 1**: As stated in the paper, the main drawback of SBN is that it could not span the entire tree topology space. It seems that there is no analysis on how/why the proposed method is better at doing this in addition to the empirical comparison in the experiments.\n\n**Response**:  Thanks for your question! In the generating process of ARTree, we use the softmax function to parameterize the conditional probability of decisions (where to add the new tip node). Therefore, all possible decisions would have nonzero probabilities. This means ARTree can sample any decision sequence with a nonzero probability. As there is a bijection between decision sequences and the entire tree topology space (Theorem 1), this implies ARTree can span the entire tree topology space. We will make it more clear in our revision.\n\n**Weakness 2**: The main contribution of the paper is in Section 3.1, which is a new parameterisation of $Q(\\tau)$. Most of the techniques in Section 3.2 follows Zhang (2023). It might be better to shorten 3.2, as I think it might not the focus of the paper.\n\n**Response**: Thanks for the suggestion! We will modify it accordingly in our revision.\n\n**Weakness 3**: Although Table 1 shows that ARTree has better number in terms of revealing the ground-truth trees, it seems that ARTree does not improve much on ELBO and ML in Table 2.\n\n**Response**: Yes, you are right! In fact, the power of ARTree for VBPI is mainly on tree topology approximation as reflected by the KL results in Table 2. There are two reasons for minor improvements on lower bounds. (i) According to our experience, the lower bounds in VBPI are more sensitive to the quality of branch length model $Q_\\psi(q|\\tau)$ than to the tree topology model $Q_\\psi(\\tau)$, and ARTree and SBN use the same branch length model in VBPI. Also, we want to clarify that significantly improving ELBO and LB-10 are difficult considering they approach the same marginal likelihood. (ii) The support of ARTree spans the entire tree topology space. This adds to the difficulty of training $Q_\\psi(q|\\tau)$ which is conditioned on tree topology $\\tau$, as discussed in Appendix E.\n\n**Question**: Can the authors talk about the potential use of the proposed method in computational biology in addition to getting better ELBO or other metrics of modelling the data? such as interpretability.\n\n**Response**: This is an interesting and open question. Two potential uses: (i) ARTree provides an alternative family of distributions over the entire tree topology with explicit likelihood computation and flexibility. This is itself a useful tool for phylogenetic inference, including tree density estimation and variational posterior approximations, which has a wide range of applications such as genomic epidemiology and conservation genetics. (ii) In ARTree, these learned conditional distributions for the species attaching operations also carry important information about the relationship between the new species to the species on the current tree topologies, and hence can be used to interpret the closeness among these species."}, {"review_id": "e5wdlHZVLQ", "reviewer": "Reviewer_aQBU", "summary": "In this paper, the authors propose a tractable distribution over tree spaces that can be fit for use in density estimation or variational inference.  The key idea is to build a tree by sequentially adding leaves one at a time by adding an additional branch to the tree.  Each step of this process has a reasonable state space, and it is easy to see that such a process generates distributions that span all of tree space.  The authors then parameterize the \"action space\" of this process using (recurrent) graph neural networks, which can be trained using either maximum likelihood in the density estimation case or VIMCO to optimize the ELBO in the VI case.  The authors apply their method to 8 standard phylogenetic benchmarking datasets finding comparable or superior performance to existing density estimation .or VI methods.", "questions": "Have the authors explored how necessary it is to condition the decisions on all previous decisions?  Does a Markov decision process perform substantially worse?  That is, does one need the recurrent GNN, or would a non-recurrent GNN be sufficient?", "limitations": "The authors have adequately addressed the limitations of their study, and I do not foresee any potential negative societal impacts.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 2, "contribution": 3, "strengths": "* The ideas presented in this paper are extremely simple and elegant.\n* The big picture of the approach is easy to describe and conceptually straightforward.\n* The performance of the method seems to be an advance over existing methods, even by metrics that favor existing methods (e.g., inclusive KL is kind to SBNs, as exclusive KL would be infinite for SBNs that do not have support on all of tree space).", "weaknesses": "* A more thorough description of the technical details of the parameterization of the model would be helpful. In particular, it would be useful to have a schematic representing all of the components and how they fit together. Equations 7-11 had a lot of subcomponents (e.g., $P$ and $R$ and $b_n$ and emb, etc... etc...) which were hard to keep track of and see how they all fit together.\n\n\nMinor:\n* I know that it is common in the field, but it is not obvious to me why one would want to take $K$ greater than $1$ in equation (4).  If $K$ is $1$, then (4) is exactly the usual ELBO.  Maximizing the $K=1$ ELBO corresponds to minimizing the KL between the variational and true posteriors, which seems desirable.  Taking $K$ larger than one certainly tightens the lower bound on the evidence, but that doesn't necessarily mean that it will result in a better variational approximation to the posterior.  In fact, as $K \\to \\infty$ equation (4) should become independent of $Q$, which seems undesirable.  See for example https://proceedings.mlr.press/v80/rainforth18b.html\n* Many of the references at the end have minor formatting issues (e.g., lacking capitalization: \"bayesian,  \"markov\", \"monte carlo\", \"Graphrnn\", etc...)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your helpful feedbacks and suggestions! Here are our responses to them.\n\n**Weakness 1**: A more thorough description of the technical details of the parameterization of the model would be helpful. In particular, it would be useful to have a schematic representing all of the components and how they fit together. Equations 7-11 had a lot of subcomponents (e.g., $P$ and $R$ and $b_n$ and emb, etc... etc...) which were hard to keep track of and see how they all fit together.\n\n**Response**: Thanks for the suggestion! We will modify our description and notations accordingly in our revision to make it more clear to the readers.\n\n**Weakness 2**: I know that it is common in the field, but it is not obvious to me why one would want to take $K$ greater than $1$ in equation (4). (...) See for example https://proceedings.mlr.press/v80/rainforth18b.html.\n\n**Response**: Thanks for asking! There are mainly two reasons for taking $K>1$. (i) The gradient of variational bound w.r.t. the discrete component $\\tau$ is generally unstable and suffers from large variance. Taking $K>1$ allows us to use efficient stochastic gradient estimators such as VIMCO (which are designed for multi-sample ELBO) for the tree topology variational parameters. (ii) A sample size $K$ larger than $1$ may encourage exploration over the vast and multimodal tree space to avoid being trapped in local modes.\n\nWe agree that taking $K$ larger does not necessarily lead to a better variational approximation, as the signal-to-noise ratio decreases as $K$ increases (https://proceedings.mlr.press/v80/rainforth18b.html). In practice, a moderate $K$ would be good choice, and we leave a more thorough investigation on the effect of $K$ to future work. Thanks for bringing up this discussion, and we will cite this paper in our revision.\n\n**Weakness 3**: Many of the references at the end have minor formatting issues (e.g., lacking capitalization: \"bayesian, \"markov\", \"monte carlo\", \"Graphrnn\", etc...)\n\n**Response**: We appreciate your careful review of the references. All the formatting issues will be carefully addressed by us in the revised version of the paper.\n\n**Question**: Have the authors explored how necessary it is to condition the decisions on all previous decisions? Does a Markov decision process perform substantially worse? That is, does one need the recurrent GNN, or would a non-recurrent GNN be sufficient?\n\n**Response**: Thanks for your insightful question! We haven't explored the option that uses Markov decision process. However, we expect it to work fairly well given that the current tree topology is a summarization of all previous decisions."}, {"review_id": "smUXLLTw8C", "reviewer": "Reviewer_vFMY", "summary": "The paper introduces ARTree, a deep generative autoregressive model for phylogenetic tree reconstruction. The authors define an autoregressive sequential process for generating tree topologies, $\\tau$, and prove that there is a one-to-one mapping between the resulting topologies and the decision sequence, $D$, instantiated by the process. Utilizing this fact, the authors define a distribution over $D$, letting each decision at time n be drawn from a Categorical distribution given previous decisions 1,..., n-1, and parameterize these Categoricals by calculating and passing learnable topological features to Graph Neural networks (GNN) with a recurrent unit and unit to incorporate time embeddings. \nARTree is then used as variational distribution for tree topologies $Q(\\tau)$, along with a GNN based parameterization of the branch length distribution, within Variational Bayesian phylogenetic inference (VBPI). The representative power of $Q(\\tau)$ is evaluated by comparing it with Subsplit Bayesian networks (SBNs) to a \"ground truth\" posterior distribution (based on the posteriors of long-running MrBayes experiments) and within the context of VBPI on benchmark datasets.", "questions": "Overall, a strong contribution to field of Bayesian phylogenetics and a well-written paper. \n\nThe addressing the following points could make me increase my score further:\n1. Experiment on runtime of ARTree in the context of VBPI\n2. Addressing the issue of minimal increase of raised in Weaknesses\n\nFailing to incorporate the related works mentioned in Weaknesses-section could make me decrease my score.\n\nMisprints:\n56 \"edges of current…\" -> \"edges of the current…\"", "limitations": "Limitations, except for potential runtime concerns (see Weaknesses-section), are properly addressed.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "The vastness of the tree topology space in phylogenetic inference is a well-known obstacle in both classical and tumor phylogenetics. The VI approach to Bayesian phylogenetics is an ongoing area of research and is need of more sophisticated tree topology variational distributions and experiment design to evaluate these varaiational distributions. The tree topology density experiment in 4.1 not only shows strong representative power of ARTree, but the experiment design itself is a contribution to the research field. Furthermore, ARTree is a large improvement when compared to previous VI-methods in phylogeny that use unconfined $Q(\\tau)$ (however this is not highlighted in the paper, see Weaknesses section).", "weaknesses": "The paper fails to mention other works within VI in phylogenetics, e.g., VaiPhy (https://arxiv.org/abs/2203.01121) and VCSMC (https://arxiv.org/abs/2106.00075); these methods do not confine the $Q(\\tau)$ support either and are relevant related work. \n\nThe paper fails to highlight the strong performance of ARTree in VBPI w.r.t. other VI methods with unconfined $Q(\\tau)$ - adding a row to Table 2 with results from, e.g., VCSMC and VaiPhy would greatly accentuate the contribution of the paper.\n\nARTree relies on several subroutines to be able to construct and parameterize the generative decision sequences. This makes the method complicated to grasp and implement, which, given the complex problem at hand, can be regarded as a strength in the authors perseverance rather than a weakness of the paper. However, the amount of steps involved, e.g., L steps of messing passing and calculating the topological node embeddings, naturally invokes questions regarding inference runtime and memory usage. The lack of runtime comparisons between VBPI with ARTree, other methods in VI and MrBayes is a weakness of the paper.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your careful review and valuable feedback! Below are our answers to your comments:\n\n**Weakness 1**: The paper fails to mention other works within VI in phylogenetics, e.g., VaiPhy and VCSMC; these methods do not confine the $Q(\\tau)$ support either and are relevant related work.\n\n**Response**: Thank you for pointing out these missing related works. We will include the comparison between ARTree and the two important related contributions, VaiPhy[1] and VCSMC[2], in our revision. Please see our global response for more discussions.\n\n**Weakness 2**: The paper fails to highlight the strong performance of ARTree in VBPI w.r.t. other VI methods (...) greatly accentuate the contribution of the paper.\n\n**Response**: Thanks for the suggestion! We will add a comparison between ARTree and other VI methods, e.g. VaiPhy and VCSMC (see the following table) in Table 2 in our revised manuscript. We did not include this in our original paper considering SBN as the SOTA model outperforms other methods significantly.\n\nTable:  Marginal likelihood (ML) estimates with one standard deviation in the brankets. $\\phi$-CSMC is proposed along with VaiPhy in [1] and works on bifurcating tree topology space, making it comparable with other methods.\n\n| ---- | ARTree | SBN | VCSMC[2] | $\\phi$-CSMC[1] |\n| ---- | ---- | ---- | ---- | ---- |\n| DS1 | -7108.41(0.19) | **-7108.41(0.15)** | -9180.34(170.27) | -7290.36(7.23) |\n| DS2 | **-26367.71(0.07)** | -26367.71(0.08) | -28700.7(4892.67) | -30568.49(31.34) |\n| DS3 | **-33735.09(0.09)** | **-33735.09(0.09)** | -37211.20(397.97) | -33798.06(6.62) |\n| DS4 | **-13329.94(0.17)** | -13329.94(0.20) | -17106.10(362.74) | -13582.24(35.08) |\n| DS5 | **-8214.59(0.34)** | -8214.62(0.40)| -9449.65(2578.58) | -8367.51(8.87) |\n| DS6 | -6724.37(0.46) | **-6724.37(0.43)** | -9296.66(2046.70) | -7013.83(16.99) |\n| DS7 | **-37331.95(0.27)** | -37331.97(0.28) | N/A | N/A |\n| DS8 | **-8650.61(0.48)** | -8650.64(0.50) | N/A | -9209.18(18.03) |\n\n[1] Koptagel, Hazal, et al. \"VaiPhy: a Variational Inference Based Algorithm for Phylogeny.\" NeurIPS 2022.\\\n[2] Moretti, Antonio Khalil, et al. \"Variational combinatorial sequential Monte Carlo methods for Bayesian phylogenetic inference.\" UAI 2021.\n\n**Weakness 3**: ARTree relies on several subroutines to be able to construct and parameterize the generative decision sequences. (...) The lack of runtime comparisons between VBPI with ARTree, other methods in VI and MrBayes is a weakness of the paper.\n\n**Response**: Thank you for pointing out the lack of runtime comparisons. We will add runtime comparisons in Appendix E in the revised manuscript. The following table is the CPU time and memory of each method in the VI setting on DS1. We do not compare them with the MCMC-based MrBayes because it seems hard to determine a fair time criterion as MrBayes is written in C++ and VBPI is written in Python.\n\nTable: The CPU time and memory usage in the VI setting on DS1. The CPU time is averaged over 100 trials with one standard deviation in the brackets. The experiments are run on a single core of MacBook Pro 2019. N/A: not available due to unresolved memory leak issues. \n| ---- |  SBN   | ARTree  | VCSMC | VaiPhy |\n| ---- |  ----  | ----  | ----  | ----  |\n| CPU time of passing 100 trees (seconds) | 0.99(0.14)  | 5.61(0.22) | 11.54(1.50)  | 34.97(0.63)  |\n| Memory (MB) | 611.74  | 605.78 | N/A  | N/A |\n\nAlthough VCSMC and VaiPhy seem to take longer time if the number of trees is fixed, they generally require hundreds of iterations to converge, since their variational distributions only have a few parameters and are highly structured. In contrast, SBN and ARTree require more than 10w iterations, since they both build machine-learning models with enormous parameters and rely heavily on optimization. ARTree takes more time than SBN because it relies on several submodules which, although complicated, are designed to promote the expressive power to accommodate the complex tree space and are widely-used strategies in the literature. The inefficiency of autoregressive generative models is also an inherent issue.\n\nThe following strategies may help to reduce the computational cost of ARTree. (i) Training on GPUs. As a deep model, ARTree is mainly implemented using vectorized tensor operations in PyTorch. (ii) Early stopping. Although ARTree is trained for 40w iterations in VBPI to get the best numerical results, 10w iterations are enough to reveal the ground truth trees. (iii) More efficient architecture. Several efforts have been made to accelerate autoregressive models, e.g. GraphGEN (https://arxiv.org/abs/2001.08184). Designing efficient architectures for ARTree is an important future direction.\n\n**Question 1**: Experiment on runtime of ARTree in the context of VBPI.\n\n**Response**: Please see our response to weakness 3.\n\n**Question 2**: Addressing the issue of minimal increase of raised in Weaknesses.\n\n**Response**: We are sorry that we could not understand this question. We could not find relevant expressions about 'minimal increase' in Weaknesses.\n\n**Question 3**: Incorporating the related works mentioned in Weaknesses-section.\n\n**Response**: Please see our response to weaknesses 1 and 2.\n\n**Question 4**: Misprints: 56 \"edges of current…\" -> \"edges of the current…\".\n\n**Response**: Thank you for your careful check. We will fix this misprint in our revision."}, {"review_id": "pJT776W9oS", "reviewer": "Reviewer_tL42", "summary": "This work proposes a new approach, ARTree, for obtaining more complex tree-topology approximations by combining deep autoregressive models and GNNs. The existing works in black-box VI for phylogenetic inference predominantly rely on SBNs as approximations of the tree-topology distribution, but here ARTree is shown to be superior in terms of KL divergence to the true distribution (obtained via MCMC). ARTree, in contrast to SBNs, does not rely on presampled tree topologies and thus explores the full tree-topology space, not a subset.", "questions": "**Critique of SBNs**: Could the authors please point to where in Zhang and Matsen, (2022) I may find the discussion regarding \"the limited parent-child subsplit patterns in the observed samples\" (line 37 in the submission)? This is an important argument for ARTree which I have not seen investigated before. To me it could make sense to include this discussion in the Appendix of this submission, as to make the paper more stand alone.\n\n", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "Due to the combinatorial nature of the tree-topology space, designing efficient and powerful density approximation algorithms of the tree-topology posterior is arguably the most complicated aspect of Bayesian phylogenetic inference. In the VI setting, SBNs have been the SOTA algorithms, but, as is mentioned in the paper, SBNs rely on presampled tree topologies. This means SBNs require candidate trees from other tree-sampling algorithms, making them not stand alone. \n\nThe algorithm proposed in the paper is an important contribution to field that is receiving increased interest in the ML community (most VI for phylogenetics papers have been published in NeurIPS, ICLR and UAI). Additionally, there have not been many attempts to improve over SBNs. As such, the proposed work appears to be relevant and could be of significant interest to the NeurIPS community.\n\nThe writing and clarity of the paper is good, although I have some questions below.\n\nThe experiments are performed using appropriate baselines (SBNs and VBPI w. SBNs). However I am missing some important references like VaiPhy (Koptagel et al., selected as oral at NeurIPS 22) and VCSMC (Moretti et al., from UAI 21). More on this below.", "weaknesses": "**Related work**: As this is a paper proposes a new algorithm for doing VI for phylogenetics, VaiPhy by Koptagel et al. (2022) should be referenced, preferably also VCSMC by Moretti et al. (2021). These work are highly related. I do not consider it necessary to experimentally compare ARTree with these algorithms though, as they do not compare favorably with VBPI empirically. \n\nFurthermore, in the VaiPhy paper a sequential algorithm (SLANTIS) is designed for sampling tree topologies. I believe it can be seen as a conditional sampler in the sense that it decided whether to replace (labeled) edges in a topology based on a Bernoulli probability, given the existing edges in the topology. It is does not use GNNs, and the algorithms are clearly distinct as SLANTIS uses precomputed maximum spanning trees. However, I think a conceptual distinction should be included in the paper, nonetheless.\n\n**Evaluating tree topologies**: It is clear that the likelihood of a topology, $\\tau$, simulated by ARTree can be efficiently evaluated. At each decision, the probability of making that decision can be computed and then the probability of proposing $\\tau$ is the product of these intermediate probabilities (as in Eq. 5; is this correct?). \n\nNow, suppose I give you another tree topology, $\\tau'$, that was not simulated by ARTree. Can ARTree compute the likelihood of $\\tau'$, i.e. $Q(\\tau')$? To me it is not clear how this is achieved by reading the paper. This seems to be like an important downstream task which SBNs can handle. If it is possible, I recommend emphasizing this feature, and how a practitioner would achieve it. If it is not possible, I think that this missing feature should be discussed. \n\nNote that I do not regard this to be a crucial feature, ARTree is still an impressive algorithm. However, it would add transparency and promote future work on ARTree.\n\n**\"Domain expertise\"**: It is repeated multiple times as a key flaw of SBNs is the required domain expertise. I do not see, and it is not explained, how ARTree diminishes this requirement? In fact, how do SBNs require more domain expertise? Running MrBAYES or UFBoot to get the presampled tree topologies can be done without understanding these algorithms in depth, as the softwares are very neatly provided. Especially, does not implementing and understanding ARTree require the same domain expertise by the practitioner as implementing and understanding SBNs?\n\nThis should be carefully clarified in the text. Alternatively, I think, the \"domain expertise\" argument should be removed as it does not add information as the paper is written at the moment. There are plenty of compelling arguments for ARTree over SBNs as is.\n\n**Experiments**: In the caption of Table 2: \"The KL results are averaged over 10 independent trainings.\" I was expecting to see standard deviations of these 10 KL numbers. How come they are not included? \n\nAlso in the same caption: \"For ELBO, LB-10, and ML,the results are averaged over 100, 100, and 1000 independent runs\". Do independent runs imply \"independent trainings\" here too? If not, why not use the 10 trained models used for the KL values to get uncertanties wrt the learned model parameters? Finally note that it says \"100, 100 and 1000\", which I figure is a typo.\n\n**Stds of ML results**: I am aware that previous works reward low-variance estimators of the marginal log-likelihood, as is done here in Table 2. My guess is that this is an appropriate way to compare models that provide estimates of lower-bounds of the ML with models that can harshly over-estimate the ML, like the stepping-stone algorithm used with MrBAYES. However, here the comparison is between two VBPI models, using either ARTree or SBNs, which both use estimates of lower bounds. Could the authors please expand on why then the standard deviation is an appropriate measure of the success of the models? \n\nFor instance, I can come up with models with estimators that have zero standard deviation by sacrificing bias. Does this make these models \"better\"?\n\nIf my concerns above are discussed and clarified, I may be willing to raise my score.\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review and helpful feedback! We address your concerns and questions as follows.\n\n**Weakness 1**: Related work\n\n**Response**: Thanks for suggesting these related works! We will reference them and clarify the distinction between SLANTIS and ARTree in our revision. More discussions can be found in the global response.\n\n**Weakness 2**: Evaluating tree topologies\n\n**Response**: Yes, ARTree can compute $Q(\\\\tau')$ even if $\\\\tau'$ is not simulated by ARTree. It is this property that allows us to calculate KL divergence to the ground truth. This fact comes from the decomposition process in Appendix C, where we will add more detailed explanations in our revision.\n\nTo calculate $Q(\\\\tau')$, we can sequentially remove the taxa one by one in a reversed order, starting from the last taxon being added (Algorithm 2). This way, we would get the corresponding decision sequence in linear time (Lemma 1) and use it to compute the likelihood of $\\tau'$.\n\n**Weakness 3**: ''Domain expertise''\n\n**Response**: Thanks for raising this issue! In terms of domain expertise, what we indeed want to emphasize is that SBNs require a pre-selected sample of good candidate trees to provide subsplit supports for parameterization. Although running MrBayes or UFBoot does not require much domain expertise, the choice of using MCMC or bootstrapping indeed demands domain expertise (see Section 4.2 of Zhang and Matsen [2022]). Moreover, those are just some heuristic approaches that are commonly used so far, and designing efficient support estimation methods for SBNs, especially when the posterior is diffuse, remains an unsolved challenge for SBN-based VBPI.\n\nWe apologize for not making this clear enough. In our revision, we will adopt your suggestion to remove the ''domain expertise'' argument for better clarification (e.g., replacing it with more specific descriptions such as pre-selected tree topologies samples).\n\n**Weakness 4**: Experiments\n\n**Response**: We are sorry that our description confused you. Let us illustrate our experiments more clearly. For each dataset, we repeat the experiment 10 times (i.e. ''trainings''). For the $i$-th repetition, we: (i) calculate the KL divergence (a deterministic number) denoted by $KL_i$; (ii) estimate the ELBO for 100 times (i.e. ''runs'') denoted by $ELBO_{i,1},\\ldots, ELBO_{i,100}$, whose sample mean is $mean_{ELBO,i}$ and sample std is $std_{ELBO,i}$. We then report $mean_{KL}=\\sum_i KL_i/10$, $mean_{ELBO}=\\sum_i mean_{ELBO,i}/10$, and $std_{ELBO}=\\sum_i std_{ELBO,i}/10$ in Table 2. The results of LB-10 and ML follow the same way as ELBO.\n\nTherefore, the stds of ELBO, LB-10, and ML across different runs reflect the variance of variational lower bounds, which is a common concern in VI. The std of KL divergence (see the following Table) across different trainings reflects the uncertainties wrt the learned model parameters, and we did not report it due to its different meaning.\n\nTable: KL divergence averaged over 10 independent trainings with standard deviation in brackets. \n|-|DS1|DS2|DS3|DS4|DS5|DS6|DS7|DS8|\n|-|---|---|---|---|---|---|---|----|\n|SBN|0.0707(0.0002)|0.0144(0.0019)|0.0554(0.0082)|0.0739(0.0012)|1.2472(0.0113)|0.3795(0.0015)|0.1531(0.0044)|0.3173(0.0257)|\n|ARTree|0.0097(0.0006)|0.0004(0.0001)|0.0064(0.0003)|0.0219(0.0014)|0.8979(0.0175)|0.2216(0.0014)|0.0123(0.0020)|0.1231(0.0078)|\n\nFinally, we want to clarify that  \"100, 100, and 1000\" is not a typo. In fact, we use 1000 runs for the ML estimation for a more accurate estimation of the variance. \n\n**Weakness 5**: Stds of ML results\n\n**Response**: In our experiments, the ML (in nats) was estimated with importance sampling \n$$\\hat{L}_K = \\log\\left(\\frac{1}{K}\\sum\\_{i=1}^{K}\\frac{P(Y,q_i,\\tau_i)}{Q(q_i,\\tau_i)}\\right)$$ \nusing $K=1000$ samples $(q_i,\\tau_i)\\sim Q(q,\\tau)$, where $Q(q,\\tau)$ is the variational approximation. With that many samples, the ML estimate $\\hat{L}_K$ is more like an exact ML $\\log p(Y)$ instead of a lower bound of it. This strategy is commonly used to access the ML of models (Normalizing Flow: http://proceedings.mlr.press/v37/rezende15.pdf; VIMCO: http://proceedings.mlr.press/v48/mnihb16.pdf). Moreover, the variance of an importance sampling estimator is often used as a measure of the approximation accuracy of the importance distribution to the target (e.g., adaptive importance sampling methods).\n\nThe importance sampling estimator $\\hat{L}_K$ is valid only when $Q(q,\\tau)=0\\Rightarrow P(Y,q,\\tau)=0$. For estimators that have zero standard deviation by sacrificing bias, it seems that this condition is violated because $Q(q,\\tau)$ would collapse to a point. Finally, please note that for ML in VBPI, the comparison of the stds is reasonable only when the means are in their correct range (see the following table).\n\nTable: ML estimates with std in the brackets.\n|-|DS1|DS2|DS3|DS4|DS5|DS6|DS7|DS8|\n|----|----|----|----|----|----|----|----|----|\n|ARTree|-7108.41(0.19)|-26367.71(0.07)|-33735.09(0.09)|-13329.94(0.17)|-8214.59(0.34)|-6724.37(0.46)|-37331.95(0.27)|-8650.61(0.48)|\n|MrBayes stepping-stone|-7108.42(0.18)|-26367.57(0.48)|-33735.44(0.50)|-13330.06(0.54)|-8214.51(0.28)|-6724.07(0.86)|-37332.76(2.42)|-8649.88(1.75)|\n\n**Question**: Critique of SBNs\n\n**Response**: Several relevant expressions can be found in Zhang and Matsen [2022]. For example, in the first two lines of page 9, it reads that ''$P_{\\pi_i}(j\\to i)$ is the conditional probability for the parent-child subsplit pair representing the local splitting pattern of ...''. Our use of the phrase ''parent-child subsplit patterns'' follows this expression. In the last 11 lines of page 11, it reads that ''if we can find a sufficiently large collection of subsplits from these favorable trees and restrict the support of CPDs accordingly ...''. This is just why we say the ''parent-child subsplit patterns'' are ''limited'' in observed samples. We will clarify this argument in Appendix A in our revised manuscript."}]

Author rebuttal:
We thank all reviewers for their constructive feedback.  We have incorporated their suggestions and will revise the paper with the following major changes:
 
**Related works**: We will clarify the distinction between ARTree and two related works - VaiPhy[1] and VCSMC[2] (see a short discussion below). Experimental comparisons will also be added to Table 2. 

 - **A short comparision**. Both VaiPhy[1] and VCSMC[2] have unconfined support over tree topology spaces. VaiPhy employs a novel sequential sampler named SLANTIS, which makes decisions on adding edges in a specific order to sample multifurcating tree topologies. Unlike ARTree which uses parametrized GNNs, SLANTIS derives decisions based on a simply parameterized weight matrix and maximum spanning trees. Secondly, VCSMC samples tree topologies through subtree merging and resampling following CSMC[3], but employs a parametrized proposal distribution. A more powerful variant of it called variational nested SMC (VNCSMC) gives better proposals by incoporating future iterations. In contrast, ARTree takes a different approach by employing GNNs for an autoregressive model that builds up the tree topology sequentially, without requiring a resampling step or a looking-forward step. Empirically, we find ARTree surpasses these two methods significantly in terms of marginal likelihood.

Table:  Marginal likelihood (ML) estimates with one standard deviation in the brankets. $\phi$-CSMC is proposed along with VaiPhy in [1] and works on bifurcating tree topology space, making it comparable with other methods.
| ---- | ARTree | SBN | VCSMC[2] | $\phi$-CSMC[1] |
| ---- | ---- | ---- | ---- | ---- |
| DS1 | -7108.41(0.19) | **-7108.41(0.15)** | -9180.34(170.27) | -7290.36(7.23) |
| DS2 | **-26367.71(0.07)** | -26367.71(0.08) | -28700.7(4892.67) | -30568.49(31.34) |
| DS3 | **-33735.09(0.09)** | **-33735.09(0.09)** | -37211.20(397.97) | -33798.06(6.62) |
| DS4 | **-13329.94(0.17)** | -13329.94(0.20) | -17106.10(362.74) | -13582.24(35.08) |
| DS5 | **-8214.59(0.34)** | -8214.62(0.40)| -9449.65(2578.58) | -8367.51(8.87) |
| DS6 | -6724.37(0.46) | **-6724.37(0.43)** | -9296.66(2046.70) | -7013.83(16.99) |
| DS7 | **-37331.95(0.27)** | -37331.97(0.28) | N/A | N/A |
| DS8 | **-8650.61(0.48)** | -8650.64(0.50) | N/A | -9209.18(18.03) |

[1] Koptagel, Hazal, et al. "VaiPhy: a Variational Inference Based Algorithm for Phylogeny." NeurIPS 2022.\
[2] Moretti, Antonio Khalil, et al. "Variational combinatorial sequential Monte Carlo methods for Bayesian phylogenetic inference." UAI 2021.\
[3] Wang, Liangliang, Alexandre Bouchard-Côté, and Arnaud Doucet. "Bayesian phylogenetic inference using a combinatorial sequential Monte Carlo method." Journal of the American Statistical Association (2015).

**Technical details**: We will provide a more schematic description of the technical details in Section 3.2 and remove redundant statements to present it more clearly to the readers.

**Limitations**: We will add a runtime comparison in Appendix E and give some suggestions for reducing the computational cost.

**Decomposition process**: We will add a description of how to evaluate the tree topology probability using the decomposition process in Appendix C.

**Drawback of SBNs**: We will remove the ambiguous argument about ''domain expertise'' and emphasize the importance of high-quality pre-sampled trees. We will clarify why SBN faces ''the limited parent-child subsplit patterns in the observed samples'' in Appendix A.

We hope our response has adequately addressed the reviewers' questions and concerns, and look forward to reading any other additional comments.

Meta-review: This paper presents a new way to construct the variational distribution of tree topologies based on autoregressive generation with GNNs, which is used in the problem of variational phylogenetic inference. 

The authors have done a pretty good job of rebuttal, and after rebuttal, this paper received scores of 6778. Most of the reviewers are confident that this paper should be accepted. Specifically, reviewers have commented that this submission is a strong contribution to the field of Bayesian phylogenetics and a well-written paper. The ideas presented are simple and elegant, and comprehensive experiments have been done in the comparison with SBN, which support the claims of the proposed method. Overall, given the general positive review feedback, the AC would like to recommend acceptance of the paper.
Decision: Accept (spotlight)