GeoPhy: Differentiable Phylogenetic Inference via Geometric Gradients of Tree Topologies
Authors: Takahiro Mimori
Michiaki Hamada
Venue: NeurIPS 2023 poster (2023)
Keywords: phylogenetic inference
variational inference
control variates
hyperbolic space
TL;DR: Novel differentiable phylogenetic inference approach with continuous topological distributions

Abstract: Phylogenetic inference, grounded in molecular evolution models, is essential for understanding the evolutionary relationships in biological data. Accounting for the uncertainty of phylogenetic tree variables, which include tree topologies and evolutionary distances on branches, is crucial for accurately inferring species relationships from molecular data and tasks requiring variable marginalization. Variational Bayesian methods are key to developing scalable, practical models; however, it remains challenging to conduct phylogenetic inference without restricting the combinatorially vast number of possible tree topologies. In this work, we introduce a novel, fully differentiable formulation of phylogenetic inference that leverages a unique representation of topological distributions in continuous geometric spaces. Through practical considerations on design spaces and control variates for gradient estimations, our approach, GeoPhy, enables variational inference without limiting the topological candidates. In experiments using real benchmark datasets, GeoPhy significantly outperformed other approximate Bayesian methods that considered whole topologies.

GeoPhy: Differentiable Phylogenetic Inference via Geometric Gradients of Tree Topologies
Takahiro Mimori
Michiaki Hamada
Waseda University
Waseda University
RIKEN AIP
AIST Waseda CBBD-OIL
takahiro.mimorf @ aoni. waseda.jp
Nippon Medical School
mhamada @ waseda.jp
Abstract
Phylogenetic inference, grounded in molecular evolution models, is essential for understanding the evolutionary relationships in biological data. Accounting for the uncertainty of phylogenetic tree variables, which include tree topologies and evolutionary distances on branches, is crucial for accurately inferring species rela- tionships from molecular data and tasks requiring variable marginalization, Varia- tional Bayesian methods are key to developing scalable, practical models; however, it remains challenging conduct phylogenetic inference without restricting the combinatorially vast number of possible tree topologies. In this work, we introduce novel, fully differentiable formulation of phylogenetic inference that leverages a unique representation of topological distributions in continuous geometric spaces. Through practical considerations on design spaces and control variates for gradient estimations, our approach, GeoPhy, enables variational inference without limiting the topological candidates. In experiments using real benchmark datasets, GeoPhy significantly outperformed other approximate Bayesian methods that considered whole topologies.
Introduction
Phylogenetic inference, the reconstruction of tree-structured evolutionary relationships between bio- logical units, such genes, cells, individuals, and species ranging from viruses to macro-organisms, a fundamental problem in biology. As the phylogenetic relationships are often indirectly inferred from molecular observations, including DNA, RNA and protein sequences, Bayesian inference has been an essential tool to quantify the uncertainty of phylogeny. However, due to the complex nature of the phylogenetic tree object, which involve both a discrete topology and dependent continuous variables for branch lengths, the default approach for the phylogenetic inference has typically been Markov-chain Monte Carlo (MCMC) method (29), enhanced with domain-specific techniques such as a mixed strategy for efficient exploration of tree topologies.
As an alternative approach to the conventional MCMCs, variational Bayesian approach for phyloge- netic inference was proposed in (42) (VBPI), which has subsequently been improved in the expressive powers of topology-depeneee presented accurate joint posterior distributions of topology and branch lengths for real datasets, they required reasonable preselection of candidate tree topologies avoid combinatorial explosion in the number of weight parameters beforehand. There have also been proposed variational approaches
branch length distributions (39 40). Although these methods have
24 [6] on top of the combinatorial sequential Monte Carlo method (CSMC(34)) ). which iteratively updated weighted topologies without the need for the preselection steps. However, the fidelity of the joint posterior distributions still largely behind that of MCMC and VBPI reported in [16
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

In this work, we propose a simple yet effective scheme for parameterizing binary tree topological distribution with a transformation of continuous distributions. We further formulate novel differen- tiable variational Bayesian approach named GeoPhy to approximate the posterior distribution of the tree topology and branch lengths without preselection of candidate topologies. By constructing practical topological models based on Euclidean and hyperbolic spaces, and employing variance reduction techniques for the stochastic gradient estimates of the variational lower bound, we have developed robust algorithm that enables stable, gradient-based phylogenetic inference. In our exper- iments using real biological datasets, we demonstrate that GeoPhy significantly outperforms other approaches, all without topologiaal restrictions, in terms of the fidelity of the marginal oog-likelihood (MLL) estimates to gold-standard provided with long-run MCMCs.
Our contributions are summarized as follows:
We propose representation of tree topological distributions based on continuous distri- butions, thereby we construct a fully differentiable method named GeoPhy for variational Bayesian phylogenetic inference without restricting the support of the topologies.
By considering variance reduction in stochastic gradient estimates and simple construction of topological distributions on Euclidean and hyperbolic spaces, GeoPhy offers unprece- dently close marginal log-likoe ligliikei estimates to gold standard MCMC runs, outperforming comparable approaches without topological preselections.
2 Background
2.1 Phylogenetic models
Let r be represent an unrooted binary tree topology with N leaf nodes (tips), and let B, denote set of evolutionary distances defined each of the branches of r. A phylogenetic tree (r, B,) represents an evolutionary relationship between N species, which is inferred from molecular data, such DNA, RNA protein sequences obtained for the species. Let Y {Yij E D N,1S/ be a set of aligned sequences with length M from the species, where Y, denote a character (base) of the i-th sequence at j-th site, and is contained in a of possible bases Ω. For DNA sequences, Ω represents a set of 4-bit vectors, where each bit represents 'A', 'T', 'G', or 'C". A likelihood model of the sequences P(Y rr B,) is determined based evolutionary assumptions. In this study, follow a common practice for method evaluations [42] 40) as follows: Y is assumed to be generated from a Markov process along the branches of in a site-independent manner; The base mutations are assumed to follow the Jukes-Cantor model [12]. The log-likelihood In P(Y B,) can be calculated using Felsenstein's pruning algorithm which is also known as the sum- cumprroduct algorithm, and differentiable with respect to B,.
2.2 Variational Bayesian phylogenetic inference
The variational inference problem for phylogenetic trees, which seeks to approximate the posterior probability PY), formulated as follows:
min DKL (2(r)Q(())))))( B,()),
where DKL. Q(T) and ((B+lT) denote the Kullback-Leibler divergence, a variational tree topology distribution, and a variational branch length distribution, respectively. The first variational Bayesian phylogenetic inference method (VBPI) was proposed by Zhang and Matsen IV [2), which has been successively improved for the expressiveness of Q (T) [3 40). For the expression of variational topology mass function Q(T). they all rely on a subsplit Bayesian network (SBN) [1]. which represents tree topology mass function (()) as product of conditional probabilities of splits (i.e.. bipartition of the tip nodes) given their parent splits. However, SBN necessitates preselection of the likely set of tree topologies, which hinders an end-to-end optimization strategy of the distribution over all the topologies and branch lengths.
OOu implementation is found at htttpps///gttuuu..oo///ta///////

---

((aa)
23
23
0.4
Q(ap)
@(za)
21
r(z)
24
21
1.5
0.4
24
1.6
1.6
2 ~ Q(z ...9 2a)
B ~ Q(B.lz)
Q(za)
22
Z2
The proposed scheme for constructing a variational distribution Q(T, Br) of tree topology and branch lengths by using distribution (=)) defined the continuous space. Left: The marginal distributions of four tip nodes (2)))..... (zz). Middle: Coordinates of the tip nodes
Figure 1:
(21..... 21} sampled from the distribution Q(2), and tree topology =(2) determined from
Right: A set of branch lengths B, (red figures) of the tree 7 is sampled from topology dependent distribution Q (B.lT))
2.3 Learnable topological features (L.TFs)
Zhang 40
introduced deterministic embedding method for tree topologies, termed learnable topological features (LTFs), which provide unique signatures for each tree topology. Specifically, given set of nodes V and edges E of tree topology node embedding function femb V - Rd is determined to minimize the Dirichlet energy defined as follows:
Σ (3,))E
- femb(v)) (v)12.
(2)
When feature vectors for the tip nodes determined, Zhang 40 showed that the optimal femb for the interior nodes can be efficiently computed with two traversals of the tree topology. They also utilized femb parameterize the conditional distribution of branch length given tree topology ((B(Tr) by using graph neural networks with the input node features (J) E V}.
2.4 Hyperbolic spaces and probability distributions
Hyperbolic spaces are known to be able to embed hierarchical data with less distortion in considerably fewer dimensions than Euclidean spaces [32]. This property has been exploited to develop various data analysis applications such as link predictions for relational data [26]. as well phylogenetic analyses [21 [9]. For representing the uncertainty of embedded data on hyperbolic spaces, the wrapped normal distribution proposed in (25] has appealing characteristics: it is easy sample from and evaluate the density the distribution.
The Lorentz model Hd is representation of the d-dimensional hyperbolic space, which is sub- manifold d d+1 dimensional Euclidean space. Denote y4 E Hd and Σ E Rgaxd location and scale parameters, respectively, random variable e H that follows wrapped normal distribution WA Σ) is defined as follows:
N.O, E), = =exx,, oPT, ur -
where nt, exp,, TH - d and PT p Thel T2 HH denote the origin, the exponential map, and the parallel transport defined on the Lorentz model. Also, the probability density WA( Σ) has a closed form and differentiable with respect to the parameters μu, More details and properties summarized in Appendix A

---

3 Proposed methods
3.1 Geometric representations of tree topology ensembles
Considering the typically infeasible task of parameterizing the probability mass function of unrooted tree topologies T. which requires (2.V 5)!! degrees freedom, we propose an alternative approach. We suggest constructing the mass function Q((T) through transformation of a certain
probability density (2)) over a continuous domain 2, follows:
Q(T) == [[[ T(2)]l.
where T 2 + T denotes a deterministic link function that maps N coordinates to the corresponding tree topology (Fig. D. Note that we have overloaded 7 to represent both variable and function for notational simplicity. An example of the representation space Z is product of the Euclidean space R.xd or hyperbolic space H. xd where denotes the dimension of each tip' representation coordinate. For the link function, we can use T(2) = TNu o Dz2), where D Z + R X denotes a function that takes N coordinates and provides a distance matrix between those based on a geometric measure such as the Fuclidean or hyperbolic distance. TNJ TRN + 'T denotes map that takes this distance matrix and generates an unrooted binary tree topology of their phylogeny, determined using the Neighbor-Joining (NJ) algorithm [3I] While the NJ algorithm offers rooted binary tree topology accompanied by estimated branch lengths, we only use the topology information and remove the root node from to obtain the unrooted tree topology T E T.
3.2 Derivation of variational lower bound
Given distribution tip coordinates Q(2) and an induced tree topology distribution Q(T) according to equation 4 the variational lower bound (I is evaluated as follows:
P(Y, (,(()) + P(T(2))
[[Q] Eq(=) (Q(R-I=(2))
In Q(B.IT(=))
(5)
Q(T(2))
Thanks to the deterministic mapping T(z), we can obtain an unbiased estimator of ([Q] by sampling from Q(2) without summing over the combinatorial many topologies However, even when
the density (2)) is computable, the evaluation of In Q(T) remains still infeasible according to
the definition 44. We resolve this issue by introducing the second lower bound with respect
conditional variational distribution R(zlT) as follows:
P(Y, B+(r(=)) In Q(B.IT(2))
+ +In P(T(3)) Q(2)
L,Q, R) = EQ(=) (Q())
(6)
This formulation is similar in structure to the variational auto-encoders [1] where Q(r[z) and R(2(T) are viewed as a probabilistic decoder and encoder of data т, respectively. However, unlike typical auto-encoders, we design Q(T[z) to be deterministic, and instead adjust Q(z) approximate Q(7).
Proposition 1. Assuming that supp R(=l7) supp Q(=l), the inequality In P(Y) LIQ) > LIQ, R| holds, where supp denotes the support of distribution. The first and the second equality holds when Q(T, B+) P(r, (Y) and R(2(T) Q=|rT), respectively.
Proof. The first variational lower bound of the marginal log-likelihhe likelihoodd is given as follows:
In P(Y) 2 In P(Y) - DKL [Q(T, B,) P((,, BY] = RQ(..B.) In P(Y, Br)
= L(Q], (7)
Q(r,
where the equality condition of the first inequality holds when Q(T, B,) = P(r, Br YY) Since we have defined Q(T) in equation (4)
further transform the lower bound 4())
In P(Y,BIIT) Q(B,lT)
+ In P(T) Q(T)
4[Q] Eq(+)
Σ [/r = -()) Eq(B, LEET
(8)
from which equation (5) immediately follows. Hence the inequality In P(Y) 2 L(Q) and its equality condition Q(T, B,) P(T, B, have been proven.

---

Next, the entropy term ( I In (((())) = BQ(r) In Q(T) in equation can be transformed
derive further lower bound as follows:
QQ()) [In Q(T)) 2 [+)))) [In Q(T) Dri (Q(2)T))
In ((1)Q(=]T)
Q(2)
(2())())))
Eo(=) In R(z[T(=))]
R(z|T)
where the last equality derived by using the relation ( = (Z)(z)) [2, 1|T [
and ((T)Q(=1T) = Q(2)I[T =(=)]. The equality condition of the first inequality holds when R(alT) = (=[TT). Hence, the inequality L[Q) 2 L(Q, R) and the equality condition is proven. n
Similar to Burda et al. [I, we can also derive a tractable eeportance-weigheee lower-bound of the model evidence (IW-ELBO), which is used for estimating the marginal log- likelihood (MLL), In P(Y), or an alternative lower-bound objective for maximization. The details and derivations are described in Appendix B
3.3
Differentiable phylogenetic inference with GeoPhy
Algorithm 1 GeoPhy algorithm
8, Ø. < Initialize variational parameters
2: while not converged do
(1:K) (1:K) Random samples from distributions Pz. PB
3:
.EB
((:KK +- ho(es (1:K)
4:
B(l:K)
holen (1:K) ((2(1))))
(1:K)
5:
6:
Go, Do: Fv t- Estimate the gradients VoL. VoL, V.L.
7:
0, tb, 1/ 4- Update parameters using an SGA algorithm, given gradients go. To- Dv
8: end while
9: return 0, o.
Here, we introduce parameterized distributions, Qo(=), Qo( and R. (2|7). aiming to optimize the variational objective L(Qu,a. R,y using stochastic gradient ascent (SGA). The framework, which we term GeoPhy, is summarized Algorithm [ While aligning with black-box variational inference 27 and its extension to phylogenetic inference 42), our unique lower bound objective enables us to optimize the distribution over entire phylogenetic trees within continuous geometric spaces.
Gradients of lower bound We derive the gradient terms for stochastic optimization of L, omitting parameters 8, c. &i for notational simplicity. We assume that Qo(2) and Qs (B)((((()) are both reparameterizable; namely, we sample from these distributions follows:
= ho(()) B, = (G T(=)),
where the terms (ex) and eg PB({B) represent sampling from parameter-free distributions Ps and PB. respectively, and the functions he and he are differentiable with respect and ø. respectively. Although we cannot take the derivative of T(2) with respect to the gradient of L with respect is evaluated as follows:
P(Y, VoL Eqo(=) (vo2)Q)()) (EQ.(2. In
+nP(((())R(((((()) Qo(B,IT(2))
(10)
where Hi denotes the differential entropy. The gradient of with respect to & and th are evaluated as follows:
VoL = Eqo(s) Epa(t) V In P(Y = [(((e)))(((2)
VeL = EQo(s) Vo In (x(
Qo ho(c,T))
(11)
where assume tractable density model for Ryg (2/7).
5

---

Gradient estimators and variance reduction From equation (TO). an unbiased estimator of VoL. using K K-saplle Monte Carlo samples, can be derived as follows:
98
KK
1 Σ K
In Qo( f(2(k) B Vo (olno((ho(ee)))
(12)
where we denote
e Ps p=2211) ho(e: EBB ~ PBB B ho(e , ((2(k))), and
(k)
f(z (,,,B))::ll P (())
P(T())),(()(((()
(13)
(B(
Note that we explicitly distinguish and he(es) to indicate the target of differentiation with respect to 0.
In practice, it is known to be crucial to reducing the variance of the gradient estimators proportional the score function V Qo to make optimization feasible. This issue addressed by introducing control variate c(z), which has a zero expectation [x[c[=)] = 0 and works to reduce the variance of the term by subtracting it from the gradient estimator. For the case of > is known to be effective to use simple Leave-one-out (LOO) control variates [15] 28 as follows:
K
(K)
VolnQu(2(1)) (I(2(t), BB!)) - [m (I Vo In Qa
I0,100
K
(14)
where we defined [x(2(1)) B(W)) respect to the k-the sample.
Ex - P (z( Bla), which is constant with K-T
We can also adopt gradient estimator called LAX 77. which uses a learnable surrogate function 8x (2) differentiable in = to form control variates as follows:
ēeLLAX (Vo In Qo (2)) (f(2 B,) - (x)) + (osy(ho(ez)) - Ve In Qo(ho(ez)), (15)
where we write the case with K = for simplicity. To optimize the surrogate function 8x. we use the estimator Fx = V. x (g8) where (.) denotes the average of the vector elements. We summarize the procedure with LAX estimator in Algorithm??
We also consider other options for the gradient estimators, such as the combinations of the LOO and LAX estimators, the gradient of IW-ELBO, and its variance-reduced estimator called VIMCO 22 in our experiments. For the remaining gradients terms v and v.,.., we can simply employ reparameterized estimators as follows:
To = V P(Y, B, = (e ()))
Te = vn Rp(27(23))
(16)
Q = ho(CB, (ao(CB.T)le T²))
where we denote cB ~ PB and z ~ Qo. More details of the gradient estimators are summarized in Appendix B
Variational distributions To investigate the basic effectiveness of GeoPhy algorithm, we employ
simple constructions for the variational distributions (z).,, (B+lT,, and (2/7)) We use an independent distribution for each tip node coordinate, i.e. Qo(2) = - Qo, (zx). where we use d-dimensional normal or wrapped normal distribution for the coordinates of each tip node 2, For the conditional distribution of branch lengths given tree topology, Q Q, (B,ll))
the diagonal lognormal distribution whose location and scale parameters are given as a function of the LTFs of the topology т, which is comprised of GNNs as proposed in (40 For the model of Rg (z(7), also
employ an independent distribution: R (2|T) = II Re, (21lT), where, we use the same type of distribution Qo, (24), independent of
"x is unbiased estimator of V (Vo()) ((oo)) as shown in 7

---

Related work
Differentiability for discrete optimization Discrete optimization problems often suffer from the lack of informative gradients of the objective functions. To address this issue, continuous relaxation for discrete optimization has been actively studied, such as a widely-used reparameterization trick with the Gumbel-softmal distribution [10] 20]. Beyond categorical variables, recent approaches have further advanced the continuous relaxation techniques to more complex discrete objects, including spanning trees [33] However, it is still nontrivial to extend such techniques to the case of binary tree topologies. As outlined in equation 4 we have introduced distribution over binary tree topologies T derived from continuous distributions (2²). This method facilitates gradient-based optimization further aided by variance reduction techniques.
Gradient-based algorithms for tree optimization For the hierarchical clustering (HC), which reconstructs a tree relationship based on the distance measures between samples, gradient- based algorithms 23 2 have been proposed based on Dasgupta's cost function [4). In particular, Chami et al. (2 proposed to decode tree topology from hyperbolic coordinates while the optimization is performed for relaxed cost function, which is differentiable with respect to the coordinates. However, these approaches are not readily applicable to more general problems, including phylogenetic inference, their formulations depend on the specific form of the cost functions.
Phylogenetic analysis in hyperbolic space The approach of embedding phylogenetic trees into hyperbolic spaces has been explored for visualization and interpretation of novel samples with existing phylogeny
For the inference task, maximum-likelihow likelihood approach was proposed in (35), which however assumed a simplified likelihood function of pairwise distances. A recent study [9 proposed an MCMC-based algorithm for sampling (T, B,), which were linked from coordinates using the NJ algorithm [3I. However, there remained the issue of an unevaluated Jacobian determinant, which posed challenge in evaluating inference objectives. Given that we only use topology т as described in equation (4. the variational lower bound for the inference can be
unbiasedly evaluated through sampling, as shown in Proposition] I
5
Experiments
We applied GeoPhy to perform phylogenetic inference on biological sequence datasets of 27 to 64 species compiled in Lakner et al. [17].
Models and training As a prior distribution of P(T) and P(B-IT), we assumed a uniform dis- tribution over all topologies, and an exponential distribution Exp(10) 10) independent for all branches, respectively, as commonly used in the literature [42 [6] For the neural network used in the parameter-
ization of Q (B= (T). we employed edge convolutional operation (EDGE), which was well-performed architecture in 40). For the stochastic gradient optimizations, we used the Adam optimizer [33 with a learning rate of 0.001. We trained for GeoPhy until one million Monte Carlo tree samples were consumed for the gradient estimation of our loss function. This number equals the number of likelihood evaluations (NLEs) and used for a standardized comparison of experimental runs 34 42). The marginal ligeliiel (MLL) value is estimated with 1000 MC samples. More details of the experimental setup are found in Appendix C
Initialization of coordinates We initialized the mean parameters of the tip coordinate distribution Qo(z) with the multi-dimensional scaling (MDS) algorithm when Qo was given as normal distribu- tions. For Qo comprised of wrapped normal distributions, we used the hyperbolic MDS algorithm (hMDS) proposed in 33 for the initialization. For distance matrix used for MDS and hMDS, we used the Hamming distance between each pair of the input sequences Y similar to [19]. For the scale parameters, we used 0.1 for all experiments. For Rig. (z), we used the same mean parameters as Qz(z) and 1.0 for the scale parameters.
5.1 Exploration of stable learning conditions
To achieve stable and fast convergence the posterior approximations, we compared several control variates (CVs) for the variance reduction of the gradients by using DSI dataset [8]. We
7

---

(diag d-2,
(diag K-3)
d=2, K-3)
-7500
-7500
estimate
eetirrate
MI. eelerratte
-8000-
-8000
LAX
LAX
LGO
-8500
LOO
LOO
LGD KLAX
NCCV LAX
-9000
VINCO
-9000
VINO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Iterations
1e5
Iterations
1e6
Iterations
166
Figure 2: Comparison of marginal log- likeliioo (MLL) estimates for DS1 dataset with different control variates. For variational distribution, Q(z), an independent two-dimensional Normal distribution with a diagonal (diag) or full covariance matrix was used for each tip node. K stands for the number of Monte Carlo samples used for gradient estimation. For reference, show the mean MIL value (gray dashed lines; -7108.42 + 0.18) estimated with the MrBayes stepping-ston (SS) method in 42. We also show MLL estimates obtained with VBPI-GNN (red dashed lines) using VIMCO estimator (K = 10). The iterations are counted as the number of likelihood evaluations (NLEs). Legend: NoCV stands for control variates, LOO for leave-one-out CV, and IW for the case of using weportence-weiighed ELBO as the learning objective.
demonstrate that the choice of control variates is crucial for optimization (Fig. 2 In particular, adaptive control variates (LAX) for individual Monte Carlo samples and the leave-one- (LOO) CVs for multiple Monte Carlo samples (with K 3) significantly accelerate the convergence of the marginal log- ligeellld (MLL) estimates, yielding promising results comparable to VBPI-GNN. Although IW-ELBO was effective when ((2) was comprised of diagonal Normal distributions, no advantages were observed for the case with full covariance matrices. Similar tendencies in the MLL estimates were observed for the other dimension (d = 3, 4) and wrapped normal distributions (Appendix D.I..
5.2
Comparison of different topological distributions
We investigated effective choices of tip coordinate distributions Q(=), which yielded tree topologies, by comparing different combinations of distribution types (normal N wrapped normal W.N), space dimensions (d = 2,3, 4), and covariance matrix types (diagonal or full), across selected CVs (LAX, LOO, and LOO+LAX). While overall performance was stable and comparable for the range of configurations, the most flexible ones: the normal and wrapped normal distributions with d = 4 and full covariance matrix, indicated relatively high MLL estimates within their respective groups (Table 2 in Appendix D). implying the importance of flexibility in the variational distributions.
5.3 Performance evaluation across eight benchmark datasets
To demonstrate the inference performance of GeoPhy, we compared the marginal lieliiol
(MLL) estimates for the eight real datasets (DS1-9 1081 8 6
37
g 17
4
38
30
The gold-standard
values were obtained using the stepping-stone (SS) algorithm 36 in MrBayes 29 wherein each evaluation involved ten independent runs, each with four chains of 10,000,00 iterations, as reported in [42]. In Table [ we summarize the MLL estimates of GeoPhy and other approximate Bayesian inference approaches for the eight datasets. While VBPI-GNN [4) employs a preselected set of tree topologies as its support set before exccution, is known provide reasonable MLL estimates near the reference values. The other approaches including CSMC 34,, VCSMC 24], CCCSMC [16) and GeoPhy (ours) tackles the more challenging general problem of model optimization by considering all candidate topologies without preselection. We compared the GeoPhy for two configurations of Q(2): a wrapped normal distribution WW with 4-dimensional full covariance matrix, and a 2-dimensional diagonal matrix, with three choices of CVs for optimization. Results for other settings, including those of the normal distributions in Fuclidean spaces, are provided in Table in Appendix [ There, noted slight advantage of wrapped normal distributions over their Euclidean counterparts. In

---

we observe that the superior performance of the larger model (fullll)4)) over (dNadiag,)) observed in T is consistently reproduced across all datasets. Moreover, although (Fflf)L,22 is relatively less performant in our comparison, almost all the estimates outperformed the other CSMC-based methods, implying the stability and efficiency of our approach.
Table 1: Comparison of the marginal log- ligelikeli (MLL) estimates with different approaches in eight benchmark datasets. The MLL values for MrBayes SS and VBPI-GN GN-G which employs the VIMCO estimator of 10 samples and EDGE GNN, are sourced from 40), The values of CSMC, VCSMC, and Ø-CSMSC are referenced from [6]. We adopt the nomenclature from the literature 42] and refer to the dataset named DS7 in [16 as DS8. The MLL values for our approach (GeoPhy) are shown for two different (=)) configurations: wrapped normal distribution WWV with 4-dimensional full covariance matrix, and 2-dimensional diagonal matrix. In each training, we used three different sets of CVs: LAX with K 1, LOO with K = 3 which we labeled LOO(3), and combination of LOO and LAX. denoted as LOO(3)+. The bold figures are the best (highest) values obtained with GeoPhy and the tree CSMC-based methods, all of which perform an approximate Bayesian inference without the preselection of topologies. We underlined GeoPhy's ' MLL estimates that outperformed the other CSMC-based methods, demonstrating superior performance across various configurations.
Dataset
DSI 27
DS3
DS4
DS5
*Taxa (N)
29
36
41
50
50
59
64
#Sites (M)
1812
1137
MrBayes SS
7108.42
26336.57
33735.44
13330.06
8214.51
6724.07
33333.77
8649.88
22 42
(0.18)
(0.48)
(0.5)
(0.54)
(0.29)
(0.86)
(2:42)
(1.75)
VBPL-GNN
7108.41
26366.77
337735.2
13329.91
8214.64
6721.3
-37332.01
8650.65
40
(0.14)
(0.07)
(0.09)
(0.19)
(0.38)
(0.4)
(0.26)
(0.45)
CSMC
8306.76
27884.37
35381.01
15019.21
8940.62
8029.55
11013.57
34
(166.27)
(226.0)
(218.15 18)
(100.61)
(46.10)
(83.67)
(113.19)
VCSMC
-9180.34
28700077
-32211.2
-17106.1
-9449.65
12296 66
24
(170.27)
(1892.67)
(397797)
(362.74)
(2578.58)
(20-6.77
+ CSMC
-7290.36
30366//0
337788006 00
13582-24
8367.51
7013.83
-9209 18
16
(7.23)
(31.34)
(6.62)
(35.08)
(8.87)
(16.99)
(18.03)
(f
7111.55
26377.99 18
-33757.79 79
13312.71
-8210.87
66335.11
33377788 80
8663 51
(0.07)
(8.07)
(1,61)
(2.64)
(29.48)
(6,85)
L0O(3)
77199.77
26368.44
41
33736.01
13339.992 26
-8231.06
-6733.91
377550777
8671. 32
(0.13)
(0.03)
(3.19)
(0.37)
(11.74)
(5.98)
L00(3)
7116.09
26368 51
33735.85
- 13337777 42
8233.89
6733590
-37358.0 96
8660 48
(10.67)
(0.12)
(1.32)
(13.06)
(0.78)
(diddigg2))
7126. 89
26-111.84
-3823.774 74
-33358 16
-8251 .5
-6745.60
-37516.88
-8719 41
(10.06)
(27.91)
(15.62)
(9.79)
(69.85
(60.54)
L0O(3)
-1300.67
26380-11
-33737.75
133316.91
82239 36
-6711.63
-37382,28
8690 11
(14.40)
(2.48)
(4,25)
(3:23)
(31.96)
L00(3)+
-1228 10
26375 28
-33336991
133317.32
8223 41
66712.40
-37411.28
8683 22
(9.78)
(11.78)
(1.91)
(1.12)
(5.70)
(1.94)
(56.74)
(13.13)
Figure 3: Superimposed probability densities of topological distributions ΣN ((z;) up to 100,000, 200,000, and 500,000 steps (MC samples) from the left. For Q(z). we employed a wrapped normal distribution with two-dimensional full covariance matrix. The experiments used the DS3 dataset (N 36). The majority-rule consensus phylogenetic tree obtained with MrBayes and each step of GeoPhy are shown in blue and red lines, respectively. The center area is magnified by transforming the radius of the Poincaré coordinates into tanh 2.1r.
We further analyze the accuracy and expressivity of tree topology distribution Q(T) in comparison with MCMC tree samples. As the training progresses, the majority consensus tree from GeoPhy

---

D51
DS3
DS7
DS1
30
1.0
ayaas
20
VBPI-NN
25
0.8
15
20
Cons RF
0.6
Freq
15
dim
0.4
10
0.2
2
0.0
-7110 -7100-711110 -7130
-33750 -33800
-37400
-37600
100
MLL estimate
estimate
estimate
Bipartitions
Figure 4: First to Third: Marginal log- lillliieehoodd (MLL) estimates and Robinson- Foulds (RF) dis- tance between the consensus trees obtained from topology distribution Q(T) and MrBayes (MCMC) for datasets DSI, DS3, and DS7, respectively. Fourth: Comparison of bipartition frequencies derived from the posterior distribution of tree topologies for MrBayes, VBPI-GNN, and GeoPhy (ours).
progressively aligns with those obtained via MCMC (Fig. 3). Furthermore, we also observed that the Robinson-Foulds (RF) distance between the consensus trees derived from the topological distribution Q(T) and MCMC (MrBayes) are highly aligned with the MLL estimates (Fig. First to Third), underscoring the validity of MLL-based evaluations. While GeoPhy provides tree topology mode akin to MCMCs with close MLL value, it may require more expressivity in Q(=) to fully represent tree topology distributions. In Fig. Fourth, we illustrate the room for improvement in the expressivity of (T)) in representing the bipartition frequency of tree samples, where VBPI-GNN aligns more closely with MrBayes. Further analysis is found in Appendix D.3 terms of execution time, we compared GeoPhy with VBPI-GNN and MrBayes, obtaining results that are promising and comparable to VBPI-GNN (Fig. TO Left). Though we have not included CSMC-based methods in the comparison due to deviations in MIL values, they tend to capitalize on parallelism for faster execution demonstrated in [16].
6 Conclusion
We developed novel differential phylogenetic inference framework named GeoPhy, which optimized variational distribution of tree topology and branch lengths without the preselection of candidate topologies. We also proposed a practical implementation for stable model optimization through choices of distribution models and control variates. In experiments conducted with real sequence datasets, GeoPhy consistently outperformed other approximate Bayesian methods that considered whole topologies.
Limitations and Future work
Although GeoPhy exhibits remarkable performance standard benchmarks without preselecting topologies, it may be necessary to use more expressive distributions for Q(z) than independent parametric distributions reach level of performance comparable to state-of-the-art VBPI or gold-standard MCMC evaluations. Another limitation of our study lies in its exclusive focus on efficiency, measured in terms of the number of likelihood evaluations (NLEs), without paying significant attention to optimizing the computational cost per iteration. The general design of our variational distribution allows us to replace the tree link function such as UPGMA as presented in Fig. 9 Right. Future work should explore alternative design choices in continuous tree topology representations and the link functions to address these limitations. Extending our framework to include complex models remains crucial. Given the potential of scalable phylogenetic methods to enhance our understanding of viral and bacterial evolution, further exploration could substantially impact public health and disease control.
10

---

Acknowledgments and Disclosure of Funding
We appreciate the valuable comments on our study provided by Dr. Tsukasa Fukunaga of Waseda Uni- versity and the anonymous reviewers. TM was partially supported by JSPS KAKENHI JP20H04239 and JST ACT-X JPMJAX22AL MH was partially supported by JST CREST JPMJCR2IFI and AMED AP22amal210SS, ae011 and JP21 JP2Igm0010008.
References
[1] Yuri Burda, Roger B. Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. In Yoshua Bengio and Yann LeCun, editors, 4th International Conference on Learning Repre- sentations, ICLR 2016, San Juan, Puerto Rico, May 2-4, 2016, Conference Track Proceedings, 2016. URL http: http://aaxxvvvv.. org/ aba/ http///aa//1...... 00519
[2] Ines Chami, Albert Gu, Vaggos Chatziafratis, and Christopher Ré. From trees to continuous embeddings and back: Hyperbolic hierarchical clustering. Advances Neural Information Processing Systems, 5065-15555 2020.
[3] Eli Chien, Puoya Tabaghi, and Olgica Milenkovic. for tree-fitting and hierarchical clustering. arXiv preprint arXiv:2205.09771, 2022.
Hyperaid: Denoising in hyperbolic spaces
[4] Sanjoy Dasgupta. cost function for similarity- hierarchical clustering. In Proceedings of the forty-eighte anmual ACM symposium on Theory of Computing, pages 118-127, 2016.
[5] Joseph Felsenstein. Maximum likelihood and minimum-steps methods for estimating evolution- ary trees from data on discrete characters. Systematic Biology, 22(3):240-249 1973.
[6] James R Garey, Thomas J Near, Michael R Nonnemacher, and Steven A Nadler. Molecular evidence for acanthocephala subtaxon of rotifera. Journal of Molecular Evolution, 43: 287-292, 1996.
[7] Will Grathwohl, Dami Choi, Yuhuai Wu, Geoffrey Roeder, and David Duvenaud. Backprop- agation through the void: Optimizing control variates for black-box gradient estimation. In International Conference on Learning Representations, 2018.
[8] Blair Hedges, Kirk D Moberg, and Linda R Maxson. Tetrapod phylogeny inferred from 18s and 28s ribosomal rna sequences and review of the evidence for amniote relationships. Molecular Biology and Evolution, ((6):60-633,, 1990.
[9] Daniel A Henk, Alex Weir, and Meredith Blackwell. Laboulbeniopsis termitarius, an ectopara- site of termites newly recognized as a member of the laboulbeniomyceets. Mycologia, 95(4): 561-564, 2003.
[10] Eric Jang, Shixiang Gu, and Ben Poole. Categorical reparameterization with gumbel-softmax. arXiv preprint arXiv: 1.N1:1111..0144. 2016.
[11] Yueyu Jiang, Puoya Tabaghi, and Siavash Mirarab. Phylogenetic placement problem: A hyperbolic embedding approach. In Comparative Genomics: 19th International Conference, RECOMB-CE 2022, La Jolla, CA, USA, May 20-21, 2022, Proceedings, pages 68-85. Springer, 2022.
[12] Thomas H Jukes, Charles R Cantor, et al. Evolution of protein molecules. Mammalian protein metabolism, 3:21-132, 1969.
[13] Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980. 2014.
[14] Diederik P Kingma and Max Welling. Auto-encoding variational bayes. arXiv preprint 131261 14, 2013.
[15] Wouter Kool, Herke van Hoof, and Max Welling. Buy 4 reinforce samples, get baseline for free! 2019.

---

[16] Hazal Koptagel, Oskar Kviman, Harald Melin, Negar Safinianaini, and Jens Lagergren. Vaiphy: a variational inference based algorithm for phylogeny. In Advances in Neural Information Processing Systems, 2022.
[17] Clemens Lakner, Paul Van Der Mark, John Huelsenbeck, Bret Larget, and Fredrik Ronquist. Efficcency of markov chain monte carlo tree proposals in bayesian phylogenetics. Systemalic biology, 57(1) (7(1):86-100) 1006-103,000 2008.
[18] Marc Law, Renjie Liao, Jake Snell, and Richard Zemel. Lorentzian distance learning for hyperbolic representations. In International Conference on Machine Learning, pages 3672- 3681. PMLR, 2019.
[19] Matthew Macaulay, Aaron Darling, and Mathicu Fourment. Fidelity of hyperbolic space for bayesian phylogenetic inference. PLoS computational biology, 2023.
[20] Chris Maddison, Andriy Mnih, and Yee Whye Teh. The concrete distribution: A continuous relaxation of discrete random variables. arXiv preprint arXiv:1611.00712 2016.
[21] Hirotaka Matsumoto, Takahiro Mimori, and Tsukasa Fukunaga. Novel metric for hyperbolic phylogenetic tree embeddings. Biology Methods and Protocols, 6(1):bpab006, 2021.
[22] Andriy Mnih and Danilo Rezende. Variational inference for monte carlo objectives. In International Conference Machine Learning, pages 2188-2196. PMLR, 2016.
[23] Nicholas Monath, Manzil Zaheer, Daniel Silva, Andrew McCallum, and Amr Ahmed. Gradient- based hierarchical clustering using continuous representations of trees in hyperbolic space. In Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 714-722, 2019.
[24] Antonio Khalil Moretti, Liyi Zhang, Christian Naesseth, Hadiah Venner, David Blei, and Itsik Pe' er. Variational combinatorial sequential monte carlo methods for bayesian phylogenetic inference. In Uncertainty in Artificial Intelligence, pages 971-981. PMLR, 2021.
[25] Yoshihiro Nagano, Shoichiro Yamaguchi, Yasuhiro Fujita, and Masanori Koyama. A wrapped normal distribution on hyperbolic space for gradient- based learning. In International Conference on Machine Learning, pages 4693-4702. PMLR, 2019.
[26] Maximillian Nickel and Douwe Kiela. Poincaré embeddings for learning hierarchical represen- tations. Advances in neural information processing systems, 30, 2017.
[27] Rajesh Ranganath, Sean Gerrish, and David Blei. Black box variational inference. In Artificial intelligence and statistics, pages 814-822. PMLR, 2014.
[28] Lorenz Richter, Ayman Boustati, Nikolas Nüsken, Francisco Ruiz, and Omer Deniz Akyildiz. Vargrad: a low-variance gradient estimator for variational inference. Advances in Neural Information Processing Systems, 33:13481 11111388113444 13492, 2020.
[29] Fredrik Ronquist, Maxim Teslenko, Paul Van Der Mark, Daniel Ayres, Aaron Darling, Sebastian Höhna, Bret Larget, Liang Liu, Mare A Suchard, and John P Huelsenbeck. Mrbayes 3.2: efficient bayesian phylogenetic inference and model choice across large model space. Systematic biology, 61)::3339554 2012.
[30] Amy Y Rossman, John M McKemy, Rebecca A Serdo--chultheiss, and Hans-Josef Schroers. Molecular studies of the bionectriaceae using large subunit rdna sequences. Mycologia, 93(1): 100-110, 2001.
[31] Naruya Saitou and Masatoshi Nei. The neighbor-joining method: new method for reconstruct- ing phylogenetic trees. Molecular biology and evolution, 4(4):406-425, 1987.
[32] Frederic Sala, Chris De Sa, Albert Gu, and Christopher Ré. Representation tradeoffs for hyperbolic embeddings. In International conference on machine learning, pages 4460-4469. PMLR, 2018
12

---

[33] Kirill Struminsky, Artyom Gadetsky, Denis Rakitin, Danil Karpushkin, and Dmitry P Vetrow. Leveraging recursive gumbel-max trick for approximate inference in combinatorial spaces. Advances in Neural Information Processing Systems, 1491999--11011,, 2021.
[34] Liangliang Wang. Alexandre Bouchard-CCöte and Arnaud Doucet. Bayesian phylogenetic inference using a combinatorial sequential monte carlo method. Journal of the American Statistical Association, 1106512):1362-3333 1374, 2015.
[35] Benjamin Wilson. Learning phylogenetic trees hyperbolic point configurations. arXiv preprint arXiv:2104. 1430, 2021.
[36] Wangang Xie, Paul O Lewis, Yu Fan, Lynn Kuo, and Ming-Hui Chen. Improving marginal likelihood estimation for bayesian phylogenetic model selection. Systematic biology, 60(2): 150-166, 2011.
[37] Ziheng Yang and Anne 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. Systematic biology, 52(5):705-716, 2003.
[38] Anne D Yoder and Ziheng Yang. Divergence dates for malagasy lemurs estimated from multiple gene loci: geological and evolutionary context. Molecular Ecology. 13(4):757-773, 2004.
[39] Cheng Zhang. Improved variational bayesian phylogenetic inference with normalizing flows. In H. Larochelle, M. Ranzato, R. Hadsell, M. Balcan, and H. Lin, editors, Ad- vances in Neural Information Processing Systems, volume 33, pages 18760-18771. Cur- ran Associates, Inc., 2020. 3B649999441166611222766666666898 .pff
nww.ine.rip.... http://eerrrr...ppppp/////////// cc/ppaper/200/ffiee//
[40] Cheng Zhang. Learnable topological features for phylogenetic inference via graph neural networks. The Eleventh International Conference on Learning Representations, 2023. URL
[41] Cheng Zhang and Frederick A Matsen IV. Generalizing tree probability estimation via bayesian networks. In Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa- Bianchi, and R. Gareett, editors, Advances in Neural Information Processing Systems 31, pages 1449-1458. Curran Associates, Inc., 2018. URL http::/papers wiw.. cc/paper/
7418 hatp://ai..llin......oooaaii.ita pdf
[42] Cheng Zhang and Frederick A. Matsen IV. Variational bayesian phylogenetic inference. In International Conference on Learning Representations, 2019. URL |Rtp:sttt//ooeeeeeiiie
net/forum?id=$JWa JR9FX
[43] Ning Zhang and Meredith Blackwell. Molecular phylogeny of dogwood anthracnose fungus (discula destructiva) and the diaporthales. Mycologia, 93(2):355-365, 2001.
13

---

A Details of Background
A.1 Summary of calculations for wrapped normal distributions
The wrapped normal distribution, proposed in [25], is a probability distribution defined on hyperbolic spaces, which is easy to sample from and evaluate its probability density at arbitrary coordinates in the hyperbolic spaces. In the following, provide detailed summary of the calculation involved in applying the wrapped normal distributions within the context of the Lorentz model of hyperbolic spaces.
The Lorentz model, denoted as Hi, represents the d-dimennstinal hyperbolic space as a submanifold of d dimensional Fuclidean space. Given u, u E Rd+1, can define the pseudo-inner product and pseudo-norm as follows:
(u, 1) == -0 + Σ
uyus
|ull
V(u,u), u)
(17)
/=1
The distance between hyperbolic coordinates 14, E Hd is defined as follows:
d(u, (4) cosh
cosh (oshh'(-0 ()),
(18)
where cosh -1 denotes the inverse function of the hyperbolic cosine function. Consider hyperbolic coordinates v, μa E HH and tangent vectors u E Ty, and v E T.H An exponential map cxp,.u) E H", a logarithm map log E T,H and a parallel transport map (T,-a(()) E T,H can be calculated as follows:
(xp = cosh([lulle.)e sinh((]u/le)
2 Tull
(19)
log, (v) cosh (a) (v aq4) Va² 11
(20)
PT PTu-pp (v) + (n-ov, u) (v 14), a a+1
(21)
where we denote a = (v,,), and cosh -1 represents the inverse function of cosh.
Given location and scale parameters denoted as μ Hd and Σ E 1dd respectively, the procedure for sampling from a wrapped normal distribution ~ WN(U, 2.) defined over H is given as follows:
exp, oPT nr (p)
Wiad NO E))
(22)
where n1º = 01,0,..., 0) denotes the origin of the H.. Note that we set no = 0, and u u= Wo:d
ThoH represents a tangent vector at the origin TzelHd.
From the sampling definition in equation 22 the probability density function WNN (4, Σ) can be derived as follows:
log WA p. E) log (w; 0, (d 1) In sinh (ull//
(23)
||
where u is defined as u = PT pT log (2). For detailed derivation, we refer to Appendix A of 25
A.2
GNN-based parameterization for variational branch length distributions
In this work, we employ variational branch length distribution (B,lTT parameterized with graph neural network (GNN) as described in [40]. In concrete, each of the branch lengths follows an independent lognormal distribution, where location and scale parameters are predicted with
GNN that takes the tree topology and the learnable topological features (LTFs) of the topology which are computed with method described in 40 Below, we summarize an architecture that we in this study.
14

---

Branch length parameterizations Let V and E, respectively represent the sets of nodes and branch edges for a given unrooted binary tree topology т. The imput to the GNN consists of node features represented by LTFs denoted as {/1e ( JoV,. These features undergo transformation L times follows:
( hev,
hev, = C Joev,) - -9 gbb) o.... o.... g. Jevv,)
(24)
where we set L = 2. The function (()) represents a GNN layer. For this function, we utilize edge convolutional layers, which will be described in more detail the following paragraph.
Next, the last node features (L) are transformed to output parameters of edge length as follows:
h MLPv (h!"))
(Vo E V+)
(25)
(o.u) MAX (Thy. Rw),
(v(v, E Ex) (V(v, Ex)
(26)
(()a):: log (v.o) = MTP A.LPP(M(u.. (M)()))))
(27)
where MLP N. MAX, and MLPg denotes multi-layer perceptron for node features with two hidden layers, the element-wise max operation, and multi-layer perceptron with hidden layer that outputs the location and scale parameter (Me: o.) of the lognormal distributions for each edge E E,, For each of the hidden layers employed in MLP, and MLP g. we set its width to 100 and apply the ELU activation function after the linear transformation of input values.
Edge convolutional layers In a previous study 44), GNN with edge convolutional layers, referred to as EDGE, demonstrated strong performance when predicting the posterior tree distributions. In EDGE, the function (()) transforms node features {h! JUEv, according the following scheme:
= g ((h))ovv.)
(28)
where ((9) is comprised of the edge convolutional operation with the exponential línear unit (ELU) activation function. Specifically, the transformation with the layer g(() is computed follows:
e2 = MLP(() (NO9 -100) -h())) Vu (T
(29)
A (((+1) AGG(O (9) wEN,(v) uur
(30)
h(++1) ELU (N(eet))
(31)
where N₂ (v) represents set of neighboring nodes connected to node in the tree topology refers to the concatenation operation of elements, MLP(() denotes a full connection layer and the exponential linear unit (ELU) activation unit, and (G represents an aggregation operation that takes the maximum value of neighboring edge features ex Vu E N (v) for each element.
B Variational Lower Bounds and Gradient Estimators
B.1 Variational lower bound
In Proposition we present that the following functional is lower bound of the marginal log- likelihood In P(Y).
fiQ, R| == B.). [In F2, B,)] = EQ(a) [Q(( F(2, B,) Q(=)]
(32)
< In P(Y),
(33)
where F and F² are respectively defined as follows:
P(Y. (+I(()) P(())) B,),, F(z, B,)
F(z, BB) Q(B.IT(=))
(34)
Q(z)
15

---

B.2
Gradient estimators for variational lower bound
The gradient of L Ryl with respect is given by [Bqn(B. P(,, B,)] In Qu(z)) In P(Y, P(Y,B,lT(E)) Qa(BBlT(2))
VoL (os))
(35)
(x)(s)) (VolnQo(z)) (F2L(B-1(())
+In (ore))
[oH]Qo
(36)
where HI denotes the differential entropy. We assume that (e(B,lT) is reparameterizable as in 40]: namely, B, can be sampled through B, ha(eB, T). where ep PB((B), where (B(
and hy denote parameter-free base distribution and a differentiable function with Φ, respectively. Consequently, the gradient of with respect to Φ is evaluated as follows:
V VoL
V Bqu(=) BQo(B, (z) [In F(z, By)]
(37)
EQu(+)) (ps(cs) V In P(Y, B, (o(en,T)(T(s)) Qo(B, ho(eB.T)) (2))
(38)
Lastly, the gradient of L with respect to can be evaluated with a tractable density model (=(())) as follows:
(x((z) (Qo(B.)s) F(z, B,)] = (Q((z) VV (o
(39)
Given samples c ps and é PB, we can compute - holez), T22). and B, = ( 7 r())).
Then, the below equations are estimators of gradients VoL, Vaf. and V S., respectively:
Je V In Qo(z) In F(z (2,,,, Br) Ve In (ha(()))),
(40)
V In (2,,F(,,, (()))) = Vn ((()(() P(Y, (,(Ee,,
(41)
Io V In F F(z, he (()))) = v In (2)T(2))
(42)
The gradients can be computed through the auto- auto-gradient of the following target:
In (=(=)) detach[/(z, B+)] + (1,,,, ho r(2)) - In Qo(ho(ez)),
(43)
where we denote f(z, BB) In F(, B,), and detach[-] refers to operation that blocks backprop- agation through its argument. For clarity in terms of differentiability with respect to the parameters, we distinguish between expressions (z, B,) and (ho(ez), he(en, r())))
B.3
Multi-sample gradient estimators
(k) ,(00 = ho(e)) and Given K set of Monte Carlo (MC) samples from Q0,((z, 2, Br), i.e. (ez
(k) _b(b) = (e ((2(k))) K we can simply estimate LoloQo,s R&l as follows:
(k)
1 K
9o K (VolnQ(=(()) fk2(1), B(()) - Vo In (n
Σ
(44)
As simple extension of equation (43). the gradients are obtained through an auto- aato-gradient computa- tion of the following target:
K
1
Σ
(n (o((t())))) (xtach(s(x(*), B( + f(x(M) ho (k) ((x())) In 2o (ho(e
K
(45)
B.4
Leave-one-out (LOO) control variates for variance reduction
For the term of K K-sample gradient estimator To proportional to the score function Vo In Qa. a leave-one-out (L.O0) variance reduction is known to be effective [15 28 which is denoted as
[()
follows:
K
((K)
Σ
Vo Qo(2(k)) (f(=(k),B(4 - Fz(2( B(Wi)) Vo In Qo
F1.00.0
K
k=1
(46)
16

---

where fk denotes:
Ju(=(Vk), BB()))
Σ
s(> B((r)
(47)
K
To employ the LOO gradient estimator for 6, the target of auto-auadaan auto-gradient computation equation 45 needs to be adjusted as follows:
1 K
Qo(((()) detach(f(2()) B [x( B(W))]
+f(2(0), hole (K) r(=())))) , (n
(48)
B.5 LAX estimators for adaptive variance reduction
Algorithm 2 GeoPhy algorithm with LAX
0, Φ. v 2.X*= <- Initialize variational parameters
2: while not converged do
(1:K) (1:K) Random samples from distributions Ps. PB
CB
2(1K) + ho(c (1:K)
4:
B(I.K) -ho( (1:K) (()(L))))
5:
6:
Do. Do: Iv + Esttmatee the gradients VoL, V VoL. VoL
Ox + v= (x3
7:
8:
0, v. x,
x + Update parameters using an SGA algorithm, given gradients go Jos Fo. Ix
9:
end while
10: return 0, b.
The LAX estimator 77 a stochastic gradient estimator based on a surrogate function, which can be adaptively learned to reduce the variance regarding the term Vo In Qa(z). In our case, the LAX estimator given as follows:
0LAX, V In Q=(z) (f(z, B+) - sx)) Vasy (ho(E))
(49)
We summarize the modified procedure of GeoPhy with LAX in Algorithm 2 As we assume Qa(z) is differentiable with respect to 2, we can also use a modified estimator as follows:
0LAX. V In Q=(z) (f(z, B+) sx())) Vosx(ho(e.)) Ve In (hoo(e))))
(50)
Since it is favorable to reduce the variance of ILAX.O, we optimize to minimize the following objective as proposed in 77:
(Voo =
V Qaiz)
0o.1
(Qo(s) D3. 1----qgtee Boo(2) 1).
(51)
ng 1=1
ng
where no denotes the dimension of 0. As the gradient in equation 49 estimator of V,L, which is dependent on X, we can use the relation E E
is given as an unbiased
Therefore, the unbiased estimator of the gradient Vx (V)6())[LAXx,)) is given as follows:
Tx
(52)
As we require the gradient of V with respect to x for the optimization, we use different objectives for auto-grate anto-gradient computation with respect to 0 and the other parameters and & as follows:
L'xx = Qo())) (detach[](z, B,)l 8x xx)) + + (ho(es)) In Qo(ho(ez)),
(53)
E. E'v = f(e,
(54)
17

---

B.6 LAX estimators with multiple MC-samples
For the cases with K MC-samples, we use LAX estimators by differentiating the following objectives:
(K)
1
Σ
In (=2(0))
(detach(S (ddoo=t))((deeaah((((())) B sx(2())) + sx(hole!))) In Qo (h)(e
K
k=1
(55)
K
(,(
f(z(h) Σ k=1
f(eth) =(2(()))). )))))))(()))))))
(56)
0,0
When combine LAX estimators with LOO control variates. the target for auto-gradient computa- tion changes to the following:
(K)
K (m]Q)((()) (dotachlf(2 (B))) FT( B( 8x (2(((())
1
Σ
LOOHLAX.0
+sx(ho(e!")) (n
(57)
We note that (i) 0,6 is not affected by the introduction of L.O control variates.
B.7 Derivation of importance-weighted evidence lower bound (IW-ELBO)
An importance-weighted evidence lower bound (IW-ELBO) 1. tighter lower bound of the log-likelihood In P(Y) than ELBO. For our model, a conventional K-sample IW-ELBO is given as follows:
P(Y, B , Q(B ,
C [Q] == E Q(e (.).... ((20))0))00BB""" B(K))
ln K
(58)
The fact that L(K) [Q] is the lower bound of In P(Y) is directly followed from Theorem 1 in M. However, our model cannot directly evaluate the mass function Q(T). we must resort to considering the second lower bound, similar the case of K 1 as depicted in Proposition [ K K-sample tractable IW-ELBO as follows:
We define the
K
(K)
Q, R = E Q2(20) Q(200,
In
F' (2() B(())
(59)
K
where F' is defined in equation (34). We will prove Theorem that C/W [Q, Rt serves as a lower
bound of the In P(Y). Although this inequality holds when K I. shown by In P(Y) > L(Q) > CIQ, R] in Proposition [ the relationship less obvious when > 1. Before delving into that, we prepare the following proposition.
Proposition 2. Given Q(z,T) as defined in equation 4 and an arbitrary conditional distribution R(2[T), it follows that
ERCH+ = T(=)]] <
(60)
where setting R(=/r) = Q(=(T) is a sufficient condition for the equality to hold.
Proof. The inequality immediately follows from the definition as follows:
En [)[ = T(=)]] < En =
(61)
Next, when we set R(2/r) (2)1T), the condition for equality is satisfied follows:
BQ [[[ = r(=)]] [l(=) r(=))=]
Q(7) =1
(62)
Q(T)
Q(T)
where have used the definition of Q(T) EQ(s) 7 (2)]] from equation 4 and the resulting relation Q(=lT)Q(7) 1[[T r(=)]Q(z).
18

---

Theorem 1. Given T2, as defined in equation and an arbitrary conditional distribution R(2 that satisfies suppR(2/1 (2|T) 2 (uppQ(=(T)) for any natural number K the following relation holds:
In P(Y) 2 S/w QQ. R 2 f|w - 1 Rl.
(63)
(=)]]] Additionally, if Fz, B,) is bounded and VT, Bn(=) = = In PYY) as K + oc.
then C [Q. R] approaches
Proof. We first show that for any natural number K > M,
c/W [Q, R| 2 (M) [Q, R..
(64)
For simplicity, we denote Q(2(k), B and F' ) as Qs and FY, respectively, in the following discussion. Let UK represent a uniform distribution over a subset with M distinct indices chosen from the K indices {I...., K}. Similar to the approach used in [I], we will utilize the following relationship:
I M
F E
(65)
Fm
m==1
Now, the inequality (64) derived as follows:
M
EQ....Q
E
In
EQi...x
In EE f4 1m
F
(66)
K
M
2 Qr
(m Jub M
In
M
Fy
(67)
m=1
M
In
F
(68)
M
172
where we have also used Jensen's inequality.
Next, we show that In P(Y) 2 F(W) [Q, R).
We again use Jensen's 's inequality as follows:
(K) [Q, R EQ
In
1
Fk
(69)
K
k=1
K
< In SIn Qu
F
In EQ [F' B,)]
(70)
K
The last term is further transformed as follows:
P(Y, ((Y)))))
In E.). B,)] = In EQ(,,))))
(71)
Q(B,(T()))Q(z)
(72)
Q(z,
P(Y, ((()(())
In())()
(73)
Q(2)
(((E7')
In
Σ EQ(2)
I[7' = r())]
(74)
Q(z)
TET
= In
Σ
[T = r(})]
(75)
NET
< In
Σ
(76)
MET
19

---

where, in the transition from the first to the second row, we employed the following relation:
Q(B-IT(2)) Σ Q(B.IT)[[ r(=) = Σ TET
Q(B.))) Q(B,l²),
(77)
FET
and we have used Proposition
2 for the last inequality.
Finally, we will show that the following convergence property assuming that B,) is bounded:
TW 1Q, R8] + In Σ P(Y, U1 = -T())]
(K o0).
(78)
F'ET
From the strong law of large numbers, it follows that almost surely as + oc:
K Ex= F% converges to the following term
EQ B2) (Ek, B&]]
Σ P(Y, +') ERisTr) 1/7 [/)=+]()]]] F'ET
(79)
where we have employed the same transformations as used from equation 71 to 75 Observe that the r.h.s term of the equation 78) equals to In P(Y) when Vr E T, E
-[))]]
which completes the proof.
Estimation of marginal log-likelihood For the estimation of In P(Y), we employ [() R] with
K 1, 000 similar to 40). From Theorem IW-ELBO ((K) Rt is at least a better lower bound of In P(Y) than ELBO L RJ. and converges to In P(Y) when Vr, R<)))) [T[ = r(=)]
According Proposition 2 this equality condition satisfied when we set R(=/T) (2/7), which
is approached by maximizing L[Q, R) with respect to R as indicated in Proposition
B.8 Gradient estimators for IW-ELBO
The gradient of IW-ELBO (K) QQ. Rsj] with respect 0 is given by
(K) V 1W (20,) (0. ((00)
wh (1:KK) B(IK)) In In (x)
B(ho))
V In Q0(2 (1:(1 (
+ E Qo.0( (2..)
(80)
where we have defined
(2(2((), BB(b))
w/ (2 (1:K) B(I:K))
(81)
Σ F'(2()), 13
K
(1:(1 ((li)K)(((Kk)) B(ia)) =
(2) B((())
(82)
K
Similarly, as FL is differentiable with respect to B,, and B, = ho(EB, T) is differentiable with respect 6, the gradient of c(W with respect to & can be evaluated follows:
K
(K) V VoL IW EQ E papc
(1-2( (,BK)),((B))) B.l:K) In F2 ha 7))
(K)
(83)
Since Vo In (()
= = (o from equation 34 an unbiased estimator of the
gradient (k)( is given as follows:
K
91w
(K) =
la (0(2(1))
(-wx(2(l:k)) ( (((l:)),((BKK))
B(Iik) + ((2(1.8), B(:KI]]
(84)
20

---

The remaining gradient estimators are given as follows:
K
91W.6 .
(x(2(1 B(ikk)y In F(2 (a ,7))
(85)
k=1
K
(UK
(x(2(1:K), In F(2()) B(I.k)J Bllik) ho(e()),T)) T)).
(86)
In total, the target for computing auto aato-gradient for these gradients is given as follows:
K
In Q0(2 detach [- -=x(2(2kk), B(ikk) + (1:Kk B(Iki)]
IW
K
+
detach(wa /ux ( (2 B(I:K))] In (( ho(en (()
(87)
B.9 VIMCO estimators
The VIMCO estimator 22 aims reducing the variance in each of the terms proportional
to ((2(:)), B(:KK) in equation (84) This is accomprished by forming for control variate
(a(2(14), B((B)) that consists of only variables other than 2(x) and ((a) for every k as follows:
K ( (3(Kk) B(L(K) + ((2(1:K) - (=((1)) ((l.K),B(..))) B(I:k))
(K)
Σ Voln (o((A)).
B(iN))]
9o, VIMCO =
(88)
where we denote that
K Fl((20))) B(A))+ ((A0) +
Σ F((e()), Bles) )
(z(2(1)) B((N) = In
(89)
K
W=
K
(BC))))B(BB=====)))))))))))))))) exp
1
Σ
In F' .B(())
(90)
K
Note that the unbiasedness of go, VIMCO can confirmed through the following observation:
E Qu. In (2(11))) B(W))
=Z1(2(14), B(W)) E Qu B
Vo Qo(
= =Z4(2(10), B(W))) (
V (o =0.
(91)
To employ the VIMCO estimator for 0, the target for auto-gradient computation as stated equation 87 needs to be adjusted as follows:
K( VIMCO
In (0()(1)) detach [-w((((:k) B((:k)) + (12(1:1) (((0(1,,,(()))) B(Ik))) - (2(11)) (0(()1),(()))
B((1))
K
detach (aachh[w((e [w (2:K) B(2:K0)) ((2(k), ho(e
(N) ((2())))
(92)
C Experimental Details
C.1 Training process
For the training of GeoPhy, we continued the stochastic gradient descent process until a total of 1,000,000 Monte Carlo (MC) tree samples were consumed. Specifically, if K MC-samples were used
21

---

per step, we performed up to 1,000,000/ / K steps. is noteworthy that the number of MC samples equaled the number of likelihood evaluations (NLEs), which provided us with basis for comparing convergence speed between different runs, shown in Fig. In all of our experiments, we used Adam optimizer with an initial learning rate of 0.0001. The learning rate was then multiplied by 0.75 after every 200,000 steps. Similar to approaches taken by Zhang and Matsen IV 42), we incorporated an annealing procedure during the initial consumption of 100,000 MC samples. Specifically, we replaced the likelihood function in the lower bound with Pr 7) and linearly increased the inverse temperature 3 from 0.001 to throughout the iterations. Note that all the estimations of marginal log-likelihood (MLL) were performed with 3 set to 1.
For each estimation of MLL values, GeoPhy and VBPI-GNN used IW-ELBO with 1000 samples, whereas all CSMC- based methods used 2048 particle weights [16].
C.2 Variational branch length distribuions
For the variational branch length distribution (+.(Bl we followed an architecture of [40]: namely, each branch length independently followed lognormal distribution which was parameterized with graph neural network (GNN). Details are described Appendix A.2
C.3 LAX estimators
As input features of a surrogate function 8x (2) used in the LAX estimators, we employed a flattened vector of coordinates 2 E IRNx xd when resides in Euclidean space. In cases where the coordinates were E HH! we first transformed with a logarithm map log E Toold, then omitted their constant value O-th elements and subsequently flattened the result. We implemented a simple multi- layer perceptron (MLP) network with single hidden layer of width 10Nd and a subsequent sigmoid linear unit (SiLU) activation function as the neural network to output s (z).
C.4 Replication of MLL estimates with MrBayes SS
Given the observed discrepancies in marginal log-likelihood (MLL) estimates obtained with the MrBayes stepping-stone (SS) method between references [42 and [16], we replicated the MrBayes SS runs using MrBayes version 3.2.7a. The script we used provided below.
BEGIN MRBAYES;
set autoclose=yes nowarn-yes Seed=123 Swapseed-123;
1set nst=1; prset statetrropr--ixxeeeppaa)) prset ertennppr ss ngen=10000000 nruns=10 nchains=4 printfreq=1000 samplefreq=100 \ savebrlens-yes silename-mrbaess.ce END;
We incorporated the results in the row named Replication in Table 3 where the values aligned more closely with those found in [42). We deduced that the prior distribution used in 16 might have been set differently as the current default values of brlenspr are Unconstrained : GammaDir((10,, 0.100, 1.0, 1.0) [] which deviates from the model assumption used for the benchmarks. We observed that the line brlenspr was not included in the code provided in Appendix F of [16). Having been able to replicate the results found in (42), we opted use their values as reference in [
C.5 Visualization of tree topologies
In Fig. we visualized the sum of the probability densities for tip node distribution Σn' Q(x;) by projecting each hyperbolic coordinate z; E HHd onto the Poincaré coordinates 3, zik/ (1 + 250) (k- = I,...,d), d), then applying transformation E₁ + tanh (a 22) / 11
pdf
22

---

with a = 2.1 to emphasize the central region. To display the density Q in the new coordinates, the Jacobian term was also considered to evaluate the density Qzz)).
For the comparison of consensus tree topologies, we plotted the edges of the tree by connecting each of their end node coordinate pairs with geodesic line. The coordinate in Hd of the i-th tip node was determined as the location parameter 14 E HH of the wrapped normal distribution Q(z;) (z; Ma; 21). Let É., E HI denotes the coordinate of an interior node u, we defined & by using the Lorentzian centroid operation C [8] as follows:
Ex
Eu = C(ccslesss fussssssssssse
(93)
V- (EnrEun
where E. = Eses, (u) M2Ca denote an unnormalized sum of weighted coordinates, 8 E (,( denote subset of tip node indices partitioned by the interior node in the tree topology T, Cs C(fmiiiee {1 hee.) denote the Lorentzian centroid of the tip nodes contained in the subset and V/ = N |s] denote the number of the tip nodes in the complement set of where |s] represents the number of tip nodes in the subset s. As an unrooted tree topology T can be identified by the set of tip node partitions introduced by the interior nodes of т, the same unrooted tree topologies give the same set of interior coordinates fEu BEv according to equation (93).
D Additional Results
D.1 Marginal log-likelihood estimates throughout iterations for DS1 dataset
In Fig. 2 we highlight several key findings of training processes. Specifically, (1) control variates
were crucial for optimizations; (2) The LAX estimator for individual Monte Carlo (MC) samples
(K = 1) and leave-one-out (LOO) CVs for multiple MC samples were effective; (3) IW-ELBO and
VIMCO were not particularly effective when a full diagonal matrix is used for Q(z). In Fig. 5 and 6
we show that these findings hold consistently across various control variates and configuration of (2),, including, different types of distributions (normal A and wrapped normal W.)), embedding
dimensions (d = 2, 3, and 4), and varying numbers of MC samples (K~ = 1, 2, and 3).
D.2 Marginal og-likelihood estimates for eight datasets
We present more comprehensive results in Table 3 extending upon the data from Table I This table showcases the marginal log- ligellkiioood (MLL) estimates obtained with various GeoPhy configurations and other conventional methods for the datusets DS1-DS8. Once again, GeoPhy demonstrates its superior performance, consistently outperforming CSMC-based approaches that do not require preselection of tree topologies across the majority of configurations and datasets. This reaffirms the stability and excellence of our approach. Additionally, we found that a ((z) configuration using 4-dimensional wrapped normal distribution with full covariance matrix was the most effective among the tested configurations.
D.3 Analysis of tree topology distributions
To assess the accuracy and expressivity of GeoPhy's tree topology distribution, we contrasted Q(T) with MCMC tree samples. For our analysis, we considered experimental runs using LAX estimators with either K = or 3. For configuration of Q().. We employed diagonal or multivariate wrapped normal distribution with 2, 4, 5, and dimensions. We sampled 1000 MC trees from Q(T) for each experiment and then deduced the majority consensus tree. As depicted in Fig, 8 across the eight datasets, the Robinson-Fould (RF) distances between the majority consensus tree of Q(T) and the one obtained with MrBayes were well aligned with the MLL values.
To highlight the expressivity limitation of Q(7) obtained from a simple independent distribution of ((2), we compiled three statistics related to the tree topology samples for DS1, DS4, and DS7 datasets (Table 4). wherein we employed the model with the least RF distance for each dataset. For DSI and DS4, GeoPhy's 's tree topologies were less diverse than those of MrBayes and VBPI-GNN and had more concentration on the most frequent topology. For more diffused tree DS7, GeoPhy also showed diverse tree samples with a close diversity index to the other two.
23

---

Table 2: Comparative results of the mean MLL estimates for DS1 dataset obtained with GeoPhy for different combinations the distribution (=)), per-step Monte Carlo samples (K), and control variates (CVs). Each of the mean MLL values is obtained from five independent runs with 1000 MC samples at the last steps. Standard deviations are shown parentheses. The bold and underlined numbers represent the best (highest) values in row-wise and column wise comparisons, respectively.
N (diag)
N (full)
CV
3
d=2
d
K=I (LAX)
-7126.79
-7130.11
-7123.95
-7129.70
71199.51
-7120.03
(10.06)
(10.63)
(12.03)
(6.14)
(10.90)
(11.92)
K=2 (LOO)
-7122.84
-7125.72
-7121.99
-7121.23
-7130.43
-7114.92
(11.96)
(13.08)
(12.52)
(14.62)
(10.82)
(8.32)
K=2 (L00OLA)))
-6129.26
-7128.18
-7116.96
-7121.74
-7129.45
-7115.14
(8.35)
(9.83)
(10.39)
(10.73)
(10:23)
(8:26)
K=3 (LO0)
-7130.60
-7124.61
-7120.88
-7132335
-1128.26
-7124.62
(10.66)
(12.21)
(13.15)
(6.89)
(9.80)
(12.33)
K=3 (L00+LAX)
7128.36
7125.66
7119.81
7122.76
7123.67
7123.37
(9.77)
(10.95)
(11.71)
(10.81)
(11.23)
(11.28)
W.N (diag)
WN (full)
CV
d
d= 3
d=
d=2
d=
K=1 (LAX)
-7126.89
7133.78
71122.10
77124.63
7119.59
7111.55
(10.06)
(13.57)
(12.29)
(8.09)
(10.84)
(0.07)
K=2 (LO0)
-7121.13
-7125.60
-7121.80
-7114.99
-7130.40
-7116.14
(13.06)
(13.04)
(12.52)
(4.66)
(10.80)
(10.66)
K=2 (LO0+LAX)
-7120.90
-7120.96
49126.49
-7119.81
-7129.58
-7116.21
(10.10)
(13.08)
(12.01)
(9.78)
(10.24)
(10.75)
K=3 (L.0)
-7130.67
77124.54
71100.94
-7123.94
7112285
-7119.77
(10.67)
(12.20)
(13.11)
(13.07)
(11.96)
(1180)
K=3 (LO0+LAX)
71288.40
-7125.69
-7125.78
77155.19
-7120.66
-7116.09
(9.78)
(13.02)
(13.10)
(8.16)
(13.12)
(10.67)
We also summarize the frequency of bipartitions observed for tree topology samples in Fig. g The intermediate values between and in this frequency plot reveal topology diversities, where VBPI-GNN aligns more closely with MrBayes compared to GeoPhy. For DS1 and DS4, GeoPhy tends to take values near zero and one the slope region. For DS7, while GeoPhy traced the curve more accurately, fluctuations around this curve highlight potential room for improvement.
D.4 Runtime comparison
We measured CPU runtimes using a single thread 20-core Xeon processor. In Fig. T0 compared single-core runtime of MrBayes, VPBI-GNN, and GeoPhy (ours). As MrBayes is written in C and highly optimized its per-iteration computation time is much faster than current variational methods. In general. the per-step efficiency of GeoPhy is comparable to that of VBPI-GNN, Although the convergence of MILL estimates tends to be slow in terms of the number of consumed samples for larger K (Fig. Left), the total CPU-time slightly decreases as K increases (Fig. Left). We also included the estimated runtime across the eight datasets and the different model configurations, where the number of species ranges from 27 to 64 (Fig. [0 Right).
24

---

Normal Idiag d=2, K-1)
Normal (diag de3, K-1]
Normal (diag d-4, K-11
-7500
-7500
MLL extirete
-8000
-8000
-8500
-8500
-8500
-9000
LAX
-9000
-9000
LAX
0.0
0.2
0.4
0.6
0.0
1.0
(dlag
(olag dm3,
(dlag Km2)
-1500
-7500
-7500
MLL extrrvat
ectimoat
MLL esttrntte
-9000
-3000
-2000
NeCY LAX
NcCV LAX
NcCV
LAX
-8500
LGO
-8500
LOG
-8500
LOG
LOO
LO0
LOG-LAX
-9000
VINGO
-9000
VINCO
-0000
VICCO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
3.0
0.0
0.2
0.4
0.6
0.8
(dlag
(clag d=3.
Idlag
K=3)
-7500
-7500
-7500
MLL ectimate
MLL cotrette
8000
-3000
-3000
Rscy
NeEV LAX
NeCV
LAX
LAX
-8500
LOO
-8500
L00
-8500
LOC
1OO
LO0-LAX
LOC-LAX
TV
-9000
VINCO
-9000
VINCO
-9000
VIMCO
0.0
0.2
0.6
0.8
1.0
0.0
0.2
0.4
0.6
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Roomaal (full dm2, K-1)
Normal d-3, K-1)
Nomaa tfull dmt, K-1)
-7500
MLL eetenaaa
MLL catmaat
-8000
-8000
-8000
-8500
-8500
-8500
McCV
NoCV
RoCV
92000
-9000
LAX
-9000
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
2.0
0.0
0.2
0.4
0.6
0.8
1.0
8 J.ommal d=,, K-2)
Normal (full d=3, K-2)
Normal (full dmf, K-2)
MLL ext2vaio
MLL cstreete
-8000
-3000
-8000
LAX
-8500
LOO
-8500
-8500
100
-9000
VINCO
-9000
Virco
-9000
VIMCO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
3.0
0.0
0.2
0.4
0.6
0.8
1.0
Normal (fall d=2, K=3)
Normal (full d=3, K=31
Normal (full d=4, K=31
MLL catmatee
MLL esttmatte
MLL exttratte
-8000
-8000
-8000
-8500
88500
-8500
-9000
9000
-900
VIMCO
0.6
0.8
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.8
165
Iterations
1c6
Figure 5: Comparison of marginal log-likeliine likelioooo (MLL) estimates for DS1 dataset using different control variates and configurations of Qz)). We configured Q(2) to be an independent normal distribution of dimensions number of MC samples was set to K 1, 2, or 3. For a detailed explanation of the legend, refer to the caption of Fig.
2, 3, with either diagonal (diag) or full covariance matrix. The
25

---

Wropped Normal Idiag d-2, K-1)
Wrapped Normal (diag de3, K-1)
Wiapped Normal (ding d-4, K-3)
MLL extirete
-8000
-8000
-8500
-8500
-8500
-9000
LAX
-900
9000
LAX
0.0
0.2
0.4
0.6
0.0
1.0
Wrapped (diag K-2)
Wrapped Normal ldlag d-3,
Wrapped Idiag K=2)
-1500
-7500
-7500
MLL extrrvat
ectimott
MLL esttrntte
-9000
-2000
-2000
NeCY LAX
NcCV LAX
NcCV LAX
-8500
LGO
-8500
LOG
-8500
LOG
LOO
LO0
LOG-LAX
-9000
VINGO
-0000
VINCO
40000
VIMCO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
3.0
0.0
0.2
0.4
0.6
088
1.0
Wrapped (diag
Wrapped Idiag d=
Wrapped Idiag
K=3)
-7500
-7500
-7500
MLL ectimate
MLL cotrette
8000
-3000
-3000
Rscy
NeEV LAX
NetV
LAX
LAX
-8500
LOO
-8500
L0C
-8500
LOC
1OO
LGO 1LAX
LOC-LAX
TV
-9000
VINCO
-9000
VINCO
-9000
VIMCO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
04
0.6
1.0
0.0
0.2
0.4
0.6
08
1.0
Wrapped Normal (full 6-2, K-1)
Wrapped Normal (full d-3, K-1)
Wrapped Nomal (full dmt, K-11
-7500
MLL eetenaaa
MLL cotmant
-8000
-8000
-8000
-8500
-8500
-8500
McCV
NoCV LAX
RoCV
90000
-9000
-9000
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
2.0
0.0
0.2
0.4
0.6
0.8
1.0
Wrapped Normal (full d=2, K-2)
Wragped Normal (full d-3, K-2)
Wragped Normal (full dmt, K-21
MLL ext2vaio
MLL cxtnetee
-8000
-3000
-8000
LAX
-8500
LO0
-8500
-8500
100
LO0
LOO-LAX
99000
VINCO
-9000
Varco
-9000
VIMcO
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
066
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Wragped Normal (Fal d=2, K-31
Wrapped Normal (full d=3, K=31
Wrapped Normal Ifll d=4, K=3)
MLL cotmatte
MLL esttmatte
MLL estiratte
-8000
-8000
-8000
-8500
-8500-
-8500
-9000
-9000
-9000
VIMCO
0.6
0.8
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.8
165
Iterations
1c6
Comparison of marginal log- likeliholdoo (MLL) estimates for DS1 dataset using different control variates and configurations of (2)). We employed wrapped normal distributions for ((2) instead of normal distributions shown in Fig. [5] The dimensions 2, 3, or and the number of MC samples K 1, 2, or 3 are the same as in the referenced figure.
Figure 6:
26

---

Table 3: Extended results of Table comparing the marginal log-lik log-liielihood (MLL) estimates with different approaches in eight benchmark datasets. The MIL values for MrBayes SS and VBPI-GNN were obtained from [40]. while CSMC, VCSMC, and 6-CSMSC are referenced from [6]. We also included replicated results for MrBayes SS. The MLL values for our approach (GeoPhy) are presented for a variety of (2)) configurations, encompassing distribution types (normal .V or wrapped normal WW), embedding dimensions (2 or 4), and the covariance matrix (full or diagonal). Each result features various CVs: LAX with K = 1, LOO with K = 3 denoted as LOO(3), and combination of LOO and LAX, denoted as LO0(3)+. The figures highlighted in bold represent the highest values obtained with GeoPhy and the three CSMC-baset methods, all of which perform approximate Bayesian inference without the preselection of topologies. We have underlined MLL estimates where GeoPhy outperformed the other CSMC-based methods.
Dataset
DSI
DS2 6
DS3 [7
DS4
DSS 17
DS6
DS7
DS8
Reference
8
0
43
38
[00
#Tasa (N)
27
29
36
59
#Sites (1f)
1949
2520
1812
1137
378
1133
1824
1008
MrBayes SS
7108.42
26367.57
333335.4
13333006
-2244.5
6724.07
333332.7
8649.88
29 42
(0.18)
(0.48)
(0.5)
(0.54)
(0.28)
(0.86)
(2.42)
(1.75)
(Replication)
7107.81
26366.45
333332.7
13328.40
8209.17
6721.55
-37331.85
8646.18
(0.25)
(0.40)
(0.63)
(0.45)
(0.46)
(0.77)
(3.08)
(1.19)
VBPI
7108.41
26366.73 73
33733512
13329.94
-2144.64
6724.37
-37332.04
8650.65
(0.14)
(0.07)
(0.09)
(0.10)
(0.38)
(0.4)
(0.26)
(0.45)
CSMC
8306.76
27884.37
335981.01
15019.21
8940.62
8029.51
11013.57
84
(166.27)
(226.6)
-3721
(218.18)
(100,61)
(16.40)
(83,667)
(113.49)
VCSMC
-9180.34
288700.7
-17106.1
(362.74)
99499.65
92966666
24
(170.27)
(4892.67)
(397.97)
(2578.58)
(0046-7)
4-CSMC
7290.36
30066.49 19
13582.24
781383
92009. 18
16
(7.23)
(31.34)
(35.08)
(8.87)
(16.99)
(18.03)
(ddiag.2)
7126.79
26:10.52 5.
3381-198
13356.21
8251.99
6777.15
33756 41
8727
(10.06)
(28. 78)
(20.31)
(9.48)
(9.43)
(15.21)
(66 28)
(43.33)
LOO(3)
7130.60
26375.10 10
33737,71
13345,55
8236.90 99
6747 46
333375,93 93
8716
(10.66)
(11.75)
(2:32)
(3.26)
(6:29)
(6.90)
(28.86)
(26.32)
L00(3)+
7128.36
26369.93 93
337735.11
18346.08
8236.18
6761 .97
37430.82
(9.77)
(0.25)
(0. 13)
(4.54)
(6.71)
(11.24)
(67.19)
(10.70)
(ddiag.4
7123.95
26382.9 91
33762.46
13341 02
-211.07
6735,77 78
-33366 04
8679.4
(12.03)
(16.89)
(10.68)
(3.61)
(6.56)
(5.25)
(26.39)
(27.55)
LOO(3)
7120.88
26308.65 53
33736.01
13338 99
8238. 16
6735.59
33357.80 86
8665.85 54
(13.15)
(0.05)
(0.09)
(6.08)
(0.52)
(4.51)
(10.76)
(5.68)
LO0(3)+
7119.81
26368.49
33735.92
13339.79
8236.09
6736.74
373553.0
8065 .09
(11.71)
(0.10)
(5.15)
(4.55)
(4.70)
(3.55)
(16.97)
(7.53)
(fful),2)
=7 7129.70
26487.71
33807.06
13353.30
-8251.01
6750.00
-37487.40
8736.81
(6.14)
(64.79)
(22.97)
(5.92)
(10,34)
(11.91)
(50.43)
(52.38)
LOO(3)
77132.35
26391 00
33736.98 98
13347.17
8237.75
6752 16
337462.07
8681 38
(6.89)
(11.8)
(1,89)
(7.77)
(6.08)
(8.64)
(54.40)
(7.98)
L00()+
7122.76
26380.59 59
33736.93 93
13343.21
8239.96
6753.88
337119.02
8691 96
(10.81)
(14.39)
(2.10)
(2.14)
(L&1)
(14.30)
(35.94)
(13.51)
((full.4)
7120.03
26378.55
337533.0
13342.22 27
82377.3
66734.51
33333333
8662.53
(11.92
(1 (11.05)
(3.03)
(2.71)
(5.41)
(1.95)
(10.08)
(4.58)
LOO()
7121.62
26368.09 19
33736.03
13333.77
88
6673149
36317746
8606
(12.33)
(0.13)
(0. 16)
(1.71)
(6.11)
(3.14)
(11.93)
(7.86)
L00(3++
77123.37
26368.51 51
33733.99
13337.06
8211.25
67341.3
-37352 30
8606 39
(11.28)
(0.09)
(0.05)
(1.45)
(8.15)
(2.18)
(12.32)
(7.54)
WW(ding.2)
77126.89
26114844
33823.74 74
13358.16
8251.45
6745.60
33516.88 88
8719.44
(10.06)
(27.91)
(15.62)
(9.79)
(9.72)
(8.36)
(69.80
(60.544
L00(3)
7130.67
26380.41
337/97.75
13346.94
8239.36
67411.6
33338222
8690.41
(10.67)
(14.40)
(2.48)
(4.25)
(4.62)
(3.23)
(31.96)
(15.92)
L
7128 10
26375 28
33736.91 91
13347.32
-233.41
66422.10 40
337411.28
8683 22
(9.78)
(1 (11.78)
(1.91)
(4.42)
(6.70)
(1.94)
(56.74)
(13.13)
WW(ding.4)
77122.10
26381 81
337559.19
13342.81
-8213.92
6833.38
-37369.30 36
8666 85
(12.29)
(17.18)
(9.95)
(3.45)
(6.74)
(0.79)
(13.45)
(10.63)
LOD(3)
7120.94
26366.52 52
337735.8
13339.77
-8236.42
- 66735.1
333441.9
8673.15
(13.1T)
(0.03)
(0.08)
13.84)
(3.63)
(2.52)
(0.15)
(0.97)
L00()+
7125.78
26368.55 51
33736.00
13342.38
8235.03
673633322 20
33345,80
8606 68
(13.10)
(0.10)
(0.18)
(6.35)
(5.36)
(1.91)
(11.13)
(5.78)
W.((fll,2)
7121163
26:18 50
3330-163
13358.16 16
-251100
67488.777 78
-371-198 98
8717 2.27
(6.09)
(2244)
(20.76)
(1.20)
(7.46)
(7.38)
(34.91)
(28.49)
L0O(3)
77125.9 91
26391.0 02
33736.99 98
13341.13
8236.90
6753 86
337116.00
8681.0 1.90
(13.07)
(11.89)
(1.96)
(0.22)
(5.13)
(10.68)
(3.12)
(12.81)
100(3)+
71155.19
26385.2 21
39736.97
18343.95
-8299.55
66747.6
337331.76
8081.55
(8.16)
(13.77)
(1.93)
(1.76)
(4.72)
(6.87)
(43.65)
(3.57)
7111.56
26379. 18
33757.79
133342.7 71
-8210.87
66355.11 14
33377.88 86
8603
(0.07)
(11.60)
(8.07)
(1.61)
(9.80)
(2.64)
(29.
(6.85)
L0O(3)
-7119.77
26308-11
33736.01
13339.26
-8244.06
66733.91
333300777
8671.32
(11.0)
(0.13)
(0,03)
(3.19)
(7.53)
(0.57)
(11.74)
(5.99)
L00(3)+
7116.09
26368.54
33735.85
13337.42
8233.89
6735.00
373558 96
-8600.46
(10.67)
(0.12)
(0.12)
(1.32)
(6.63)
(1.13)
(13.06)
(0.78)
27

---

DS1
DS3
DS4
17.5
14
25
20
20
12.5
10
Cons RF
10.0
8
15
10
dm
7.5
6
10
5.0
2.5
2
0
-7110 -7110-7120- -7130
-26400-26450-26500 estimate
-33750 -33800
-13340 -13360
MLL estimate
MLL estimate
estimate
DS5
DS8
30
17.5
30
25
15.0
25
35
12.5
20
20
Cons RF
7.5
15
20
10
2.5
5
0.0
-8220 -8240 -8260
-6740
-676C
-37400
-37600
8700
8800
MLL estimate
MLL estimate
MLL estimate
MLL estimate
Figure 7: Marginal log-like log-likelihood (MLL) estimates and Cobinson-Foulle (RF) distance between the consensus trees obtained from topology distribution Q(T) and MrBayes (MCMC) for datasets DSI-DS8. Each dot denotes an individual experimental run. The term 'dim' refers to the dimension of wrapped normal distributions used for Q(2). To avoid dot overlap the plot, RF metrics are randomly jittered within =0.2. The dashed line illustrates the MLL value estimated by MrBayes.
DS1
DS7
1.0
VSPIC
VaPiOWK
0.8
0.6
0.6
0.6
Fre3
Frog
0.4
0.4
0.4
0.0
0.0
200
200
Bipartitions
Bipartitions
Bipartitions
Comparison of bipartition frequencies from the posterior tree topology distributions of MrBayes, VBPI-GNN, and GeoPhy (ours), for DS1, DS4, and DS7 datasets. The bipartitions were ordered by descending frequency seen in MrBayes.
Figure 8:
28

---

Table 4: Comparative results of tree topology diversity statistics for DS1, DS4, and DS7 datasets using MrBayes, VBPI-GNN, and GeoPhy. Statistics presented include (a) the Simpson's diversity index, (b) the frequency of the most frequent topology, and (c) the number of topologies accounting for the top 95% cumulative frequency.
Dataset
Statistics
MrBayes
VBPI-GNN
GeoPhy
0.874
0.891
0.362
DSI
(b) 0,268
0.313
0.793
(c) 42
11
(a) 0.895
0.895
0.678
DS4
(b) 0.277
0.285
0,554
(c) 208
90
58
(a) 0.988
0.991
0.996
DS7
(b) 0.022
0.028
0.022
(c) 753
315
553
Normal (full d-2)
Normal (ful d-2, K-3)
-7200
-7400
-7800
K62
K<8 (LOO]
K43 [LOO+LAX]
(LOO)
8200
K=100.001
IPGMA (1000
-8400
(LO0(LO
UPGMA [LDO +LAX]
0.00
0.25
0.50
0.75
1.00
0.00
0.25
0.50
0.75
1.00
Itarations
1e6
Iterations
le6
Figure Left: Comparison of convergence curves for different per-step MC samples (K 2,3,5, 5, 10) in dataset DS1. Right: Comparison of convergence in dataset DS1 for different link functions, NJ and UPGMA, that map coordinates into tree topologies r(2).
DS1
K
Mrtlayes 55
400
30
VBPL-GNM
Oun/Ke))
20
OursK
MI2.0001
WIGGLE WHIL-UDDLLL
25000 50000 75000
DSS
Runtime (sec)
Figure 10: Left: Single-core runtime comparison of MrBayes, VBPI-GNN, and GcoPhy (ours). MrBayes SS* represents one of ten runs used for MrBayes SS. VBPI-GNN trained for 400k iterations with K 10 MC samples. For VBPI-GNN* and GcoPhy (ours), we represent runtimes for training with IM MC samples. For GeoPhy, we employed a 2d-diagonal normal distribution for Q(z) and LOO+LAX for control variates (CVs). Right: CPU-time per 1k MC samples for two types of distribution (2)) across eight datasets. We compared LOO or LOO+LAX for CVs.
29


Reviews:
[{"review_id": "dHVPzrhjhw", "reviewer": "Reviewer_AJKq", "summary": "This paper proposed a family of implicit distribution over tree topologies which allows support free variational Bayesian phylogenetic inference. The distribution is constructed based on the neighbor joining algorithm which maps a distribution over the tip node vectors to the tree topology space. Both Euclidean and hyperbolic spaces are considered for the tip node distributions. As the resulting distribution over tree topologies is intractable, an auxiliary reverse distribution was introduced which leads to a looser lower bound for optimization. Various gradient estimators were investigated for the tree topology variational parameters optimization. Experiments on a benchmark of real data problems demonstrated the advantage over other variational approaches (mainly sequential monte carlo based) that does not require preselected tree topololgies.", "questions": "1. In equation 10 and equaiton 12, the expection of the last term is actually zero. Any idea why this term is kept there given that other control variates are considered?\n2. How many samples are used to compute the marginal likelihood estimates?\n3. The current construction of the reverse distribution $R(z|\\tau)$ seems not take the constraints $\\tau(z)=\\tau$ into consideration directly. Would this have a negative impact on performance?\n4. I note that only low dimension tip node vectors ($d\\leq 4$) are considered in this paper. Would higher dimension embedding be helpful?\n5. How about the tree topology approximation given by the implicit distribution? Although it may be hard to estimate the density, one can also report summary statistics from the sampled trees, for example.", "limitations": "yes, they do.", "rating": 5, "confidence": 5, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. The paper is written clearly and well organized. The problem of constructing flexible families of distributions over tree topologies is important for the phylogenetic communities and would to of interest to many practitioners.\n2. In addition to a careful derivation to the training objective, the authors also investigated different gradient estimators for tree topology variational parameters.", "weaknesses": "1. As admitted by the authors, the proposed variational approximation $Q(z)$ and the reverse distribution $R(z|\\tau)$ is simple, which may damage the overall approximation quality of the method. \n2. In the derivation of the lower bound, the reverse distribution $R(z|\\tau)$ needs to satisfy the constraint: the support of $R(z|\\tau)$ should be inside the region $\\{z: \\tau(z)=\\tau\\}$.\n2. In the experiments, it seems that the performance would be really sensitive to the choice of gradient estimators. I wonder if that is related to the small sample size $K$? The marginal likelihood estimates also tends to have large variance in most cases.\n3. It may be hard to assess the performance on the tree topology posterior estimation as the distribution is implicit and defined in a rather sutble way.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful review and feedback.\n\n> W1. As admitted by the authors, the proposed variational approximation $Q(z)$ and the reverse distribution $R(z | τ)$ is simple, which may damage the overall approximation quality of the method.\n\nWe believe that enhancing the expressiveness of $Q(z)$ and $R(z | \\tau)$ is a valuable and non-trivial direction for future work. However, despite the simplicity of our choice, we consider it a significant contribution that we can consistently outperform CSMC-based methods (Table 3).\n\n> W3. In the experiments, it seems that the performance would be really sensitive to the choice of gradient estimators. I wonder if that is related to the small sample size K? The marginal likelihood estimates also tends to have large variance in most cases.\n\nThank you for your detailed feedback.\nWhile the variance of our methods is still higher than that of goldstandard MCMC runs,\nwe have confirmed that we can achieve better optimization than CSMC-based methods in most cases \nby using either $K=1$ (LAX), $K=3$ (LOO), or $K=3$ (LOO+LAX) (Table 2). We are interested in whether increasing \n$K$ could help reduce variance. However, we observed that increasing $K$ requires more iterations to achieve the same accuracy (Fig. R3 Third).\n\n> Q1. In equation 10 and equaiton 12, the expection of the last term is actually zero. Any idea why this term is kept there given that other control variates are considered?\n\nThank you for your question.\nIt may be confusing because we refer to $\\nabla_\\theta Q_\\theta(z)$ as the score function.\nWhile the expectation of $\\nabla_z \\log Q_\\theta(z)$ is zero, the term with $\\nabla_theta$ remains. \nWe will include a note in the footnote for clarification.\n\n> Q2. How many samples are used to compute the marginal likelihood estimates?\n\nWe used 1000 MC samples. We will move the description from Table 1 caption to section 5.1 Experimental setup for more clarity.\n\n> W2. In the derivation of the lower bound, the reverse distribution $R(z | \\tau)$ needs to satisfy the constraint: the support of $R(z | \\tau)$ should be inside the region $z: \\tau(z)=\\tau$.\n> Q3. The current construction of the reverse distribution \n seems not take the constraints \n into consideration directly. Would this have a negative impact on performance?\n\nThank you for highlighting this important point.\nInitially, we believed that the expressiveness of $R$ is only beneficial when a more complex $Q$ is used.\nAt least given that $L[Q] \\geq L[Q, R]$,\nit is not expected that we overestimate the quality of $Q$.\nHowever, there might be cases like\n$L[Q] \\geq L[Q'] \\geq L[Q', R] \\geq L[Q, R]$,\nso even with a simple $Q$,\ndesigning $R$ to meet the conditions could prevent suboptimal trnasitions like $Q \\to Q'$.\nWe would like to include the discussion in our revised manuscript.\n\n> Q4. I note that only low dimension tip node vectors (\n) are considered in this paper. Would higher dimension embedding be helpful?\n\nThank you for your question. There is an advantage in using a low dimension as it reduces the number of parameters. To answer the question empirically, we have added experiments for $Q(z)$ with dimensions 5 and 6 (Fig. R2 Left to Right). While the difference between dims 2 and 4 is more pronounced, there appears to be a slight improvement in DS3 and DS7 when the dimension is increased to 5 or 6.\n\n> W4. It may be hard to assess the performance on the tree topology posterior estimation as the distribution is implicit and defined in a rather sutble way.\n> Q5. How about the tree topology approximation given by the implicit distribution? Although it may be hard to estimate the density, one can also report summary statistics from the sampled trees, for example.\n\nAs we can sample easily from $Q(z)$, sampling \nfrom $Q(\\tau)$ is straigtfoward.\nTherefore, evaluating the empirical distribution is not an issue.\nEach of the consensus tree used in Fig. R1 and R2 was obtained from 1000 samples obtained from $Q(\\tau)$"}, {"review_id": "vH4rFLEZmo", "reviewer": "Reviewer_MJWZ", "summary": "Authors proposed GeoPhy as a fully differentiable approach for phylogenetic inference, addressing a fundamental problem in phylogenetic inference. In experiments with real benchmark datasets, GeoPhy demonstrated its superior performance compared to other methods when considering all topological candidates. This approach is of interest to the general ML community.", "questions": "Major questions:\n\n1. In Section 3.1, the neighbor-joining algorithm is defined as the projection from the distance matrix to the tree topology. However, other standard algorithms, such as UPGMA, could be alternative approaches. Have you investigated if the choice of projection leads to different performance in the benchmark tests?\n\n2. Appendix C.5 provides a rough estimation of the running time, but it is unclear how the other baseline models included in the benchmark tests perform in terms of time. Furthermore, how does the algorithm scale as the number of tips and sites increase? Have you observed any trends in the time curve with respect to the number of taxa and sites?\n\n3. In Table 2, GeoPhy achieves superior performance in the benchmark tests based on likelihoods. Did the optimized topologies look significantly different across the different methods?\n\n4. How were the per-step Monte Carlo samples (K) determined in the benchmark tests? There is no clear winner based on Table 1.\n\n5. Could you elaborate on the approaches used to calculate the marginal log-likelihood (MLL) for all the listed methods?", "limitations": "Authors briefly mentioned two limitations of the proposed model: 1) Q(z) can be improved to enhance its expressive power, and 2) the computation cost per update step can be optimized to speed up the overall inference process. These practical concerns can be further addressed in future studies.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "Overall, the manuscript is well-written.", "weaknesses": "It is highly recommended to include an algorithm block summarizing the GeoPhy method from input to output for clarity. Additionally, there are some details in the experimental section that need clarification.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> W1. It is highly recommended to include an algorithm block summarizing the GeoPhy method from input to output for clarity. Additionally, there are some details in the experimental section that need clarification.\n\nThank you for your suggestion. We will include the algorithm block that summarizes the method workflow in our revised manuscript.\n\n> Q1. In Section 3.1, the neighbor-joining algorithm is defined as the projection from the distance matrix to the tree topology. However, other standard algorithms, such as UPGMA, could be alternative approaches. Have you investigated if the choice of projection leads to different performance in the benchmark tests?\n\nThank you for your suggestion. Our framework defined in equation (4) certainly allows an alternative link function to obtain a tree topology.\nWe have included an experiment with UPGMA with a promising optimization trajectory (Fig. R3 Fourth).\n\n> Q2. Appendix C.5 provides a rough estimation of the running time, but it is unclear how the other baseline models included in the benchmark tests perform in terms of time. Furthermore, how does the algorithm scale as the number of tips and sites increases? Have you observed any trends in the time curve with respect to the number of taxa and sites?\n\nThank you for your suggestion.\nwe have analyzed the runtime of MrBayes for MCMC and Stepping-Stone (SS) method, VPBI-GNN, and our methods with their standard use we have employed in this work for dataset DS1 (Fig. R3 First).\nAs MrBayes is written in C and highly optimized its per-iteration computation time is much faster than current variational methods.\nRuntimes of our methods per iteration (number of likelihood evaluations) were comparable to VBPI-GNN.\n\n> Q3. In Table 2, GeoPhy achieves superior performance in the benchmark tests based on likelihoods. Did the optimized topologies look significantly different across the different methods?\n\nIn our additional experiments, we have confirmed that the discrepancy measure (Robinson-Foulds distance) of the majority-rule consensus trees derived from the topological distribution $Q(\\tau)$ and MCMC (MrBayes) are highly correlated with the marginal log-likelihood (MLL) estimates (Fig. R2 First to Third).\nWe have not evaluated topology distributions obtained with methods whose MLL estimates were far deviated from the reference values.\nAlthough when the consensus tree of our inference runs exactly matches that of MCMC, it is still difficult to express the distribution of tree topologies faithfully due to the limitation of expressivity in $Q(z)$. \nWe exemplify the case with DS1 dataset in Fig. R2 Fourth, \nwhere VBPI-GNN shows more aligned bipartition frequencies to MrBayes than ours.\n\n> Q4. How were the per-step Monte Carlo samples (K) determined in the benchmark tests? There is no clear winner based on Table 1.\n\nThank you for your question.\nWe have observed that the smaller $K$ shows relatively faster convergence in terms of the number of iterations.\nWe included the corresponding experiment in Fig. R3 Third.\n\n> Q5. Could you elaborate on the approaches used to calculate the marginal log-likelihood (MLL) for all the listed methods?\n\nFor VBPI-GNN, the importance-weighted ELBO with 1000 iterations was used for the MLL estimates. \nFor CSMC-based methods, we will include the details of computations in our revised manuscript."}, {"review_id": "wQ5YXVgUsw", "reviewer": "Reviewer_Bru6", "summary": "The authors present a novel variational distribution for tree topologies, $Q(\\tau)$, in the context of variational inference (VI) in phylogenetics. They construct their $Q(\\tau)$ by introducing a continuous distribution $Q(z)$ in hyperbolic space and define $Q(\\tau)$ by an expectation over the support of $Q(z)$ of an indicator function acting on a mapping from the hyperbolic space to the topology space; this allows sampling from $Q(z)$ and reconstructing \\tau through the neighbour joining algorithm and, by the differentiable $Q(z)$ distribution, enables Monte-Carlo gradient estimation to maximize the Evidence Lower Bound (ELBO). \n\nAn ablation study is conducted to evaluate different $Q(z)$ distributions (Normal and Wrapped-Normal) and control variates. The model is then compared to other VI methods and MrBayes across 8 datasets.", "questions": "The paper is a significant, novel and promising contribution to the field and is in its current format acceptable for publication.\n\nHowever, based on my concerns in the weakness section, addressing the following points in the rebuttal could make me increase my score:\n1. A discussion regarding the low variance in Table 2 and limited knowledge of the representative power $Q(\\tau)$.\n2. Clearer reporting of the runtime of the algorithm in the main paper.\n\nAddressing the following points in the rebuttal could increase my score even further:\n1. A experiment demonstrating the representative power/limitations of $Q(\\tau)$\n\n\nSome misprints: \"Monte-Calro\" figure text of Figure 2 line 4.", "limitations": "Limitations has been somewhat discussed (see weaknesses for omitted discussion of possible limitation).", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "The problem of phylogenetic tree reconstruction is paramount to understanding evolution, e.g., the evolution of species and cancer. The grand size of the tree topology space induced by a set of taxa constitutes a major obstacle for inference; any method successful in efficient exploration of this space is a significant contribution to the field. Here, the authors construct a differentiable $Q(\\tau)$, enabling the use of gradients to indirectly maneuver the tree topology space. Furthermore, the proposed framework does not require restricting the tree topology space based on pre-processing steps as in VBPI. \n\nThe novel, differentiable approach to $Q(\\tau)$, the avoidance of restrictions of the tree support and reporting strong results on datasets DS1-DS8 in terms of the ELBO(Q,R), makes the paper a significant and original contribution to the field of Bayesian phylogenetics.", "weaknesses": "The standard deviations of the proposed method reported in Table 2 are low for multiple datasets; this may be due to the consistency of the optimization across seeds, but may also indicate that the support of $Q(z)$ collapses to regions containing one/few unique $\\tau(z)$. Since there is no analysis/discussion regarding this, the paper presents limited information regarding the representative power of $Q(\\tau)$, which is a major set back of the paper, especially as the MrBayes \"golden-run\" posteriors over the tree topologies for DS1-DS8 support multiple topologies (Whidden and Matsen 2015, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4395846/ ).\n\nThe use of $Q(z)$ makes the method extra susceptible to limited topology support by the inherent mode-seeking behaviour in VI; while a mode in the topology space could represent multiple tree topologies, a mode in Z-space might only support reconstruction of, at worst, one tree topology. \n\nThe paper would benefit from reporting runtime compared to MrBayes and other VI methods, especially as one of the strengths in VI is speed. The current report on runtime is limited and is relegated to Appendix C.\n\nFurthermore, an obvious weakness is the need to introduce a lower bound of the ELBO.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thoughtful suggestion and feedback.\n\n> W1. The standard deviations of the proposed method reported in Table 2 are low for multiple datasets; this may be due to the consistency of the optimization across seeds, but may also indicate that the support of $Q(z)$ collapses to regions containing one/few unique $\\tau(z)$.\n\n> Q1-1. A discussion regarding the low variance in Table 2 and limited knowledge of the representative power $Q(\\tau)$.\n\nThank you for your question. We expect that the main source of relatively low variance in Table 2 of our methods comes from that the marginal log-likelihoods (MLLs) are sufficiently close to the reference values.\nIndeed, although VBPI-GNN needs the reasonable preselection of candidate topologies, it exhibits the lowest variance among the methods (Table 2) and still shows aligned bipartition frequencies to those of MCMC (Fig. R2 Fourth).\nThe future work should address the expressivity of $Q(z)$ to represent a more diverse topological distribution $Q(\\tau)$.\n\n> The paper would benefit from reporting runtime compared to MrBayes and other VI methods, especially as one of the strengths in VI is speed. The current report on runtime is limited and is relegated to Appendix C.\n\n> Q1-2. Clearer reporting of the runtime of the algorithm in the main paper.\n\nThank you for your suggestion. We have compiled the runtime of our methods in comparison with MrBayes and VBPI-GNN (Fig. R3 First).\nWhile we expect that more efficient per-step computation will be an important future work, current performance for the fixed number of iterations is comparable to that of VBPI-GNN.\nWe also included the estimated runtime across the eight datasets \nand the different $Q(z)$ model configurations,\nwhere the number of species ranges from $N=27$ to $64$ (Fig. R3 Second).\n\n> Q2-1. A experiment demonstrating the representative power/limitations of Q(τ)\n\n> The use of Q(z) makes the method extra susceptible to limited topology support by the inherent mode-seeking behaviour in VI; while a mode in the topology space could represent multiple tree topologies, a mode in Z-space might only support reconstruction of, at worst, one tree topology.\n\nThank you for your constructive suggestion.\nWe confirm that the mode of the topological distribution $Q(\\tau)$ is well matched to that of MCMC when the MLL values are close to goldstandard, through the visualization of consensus trees obtained from $Q(z)$ and MCMC shown in Fig. R1 and\nthe evident correlation between Robinson-Foulds (RF) distance and MLL estimates for multiple datasets in Fig. R2 First to Thrid.\nFor $Q(\\tau)$ in DS1 dataset, we showcase bipartition frequencies observed in the posterior tree distribution of MCMC, VBPI-GNN, and our methods (Fig. R2 Fourth).\nThe divergence index of our $Q(\\tau)$ in Fig. R2 Fourth is 0.36 > 0.\nHowever, the divergence is considerably lower than those of MCMC (0.87) and VBPI-GNN (0.86),\nand it is still difficult for our method to represent the tree distributions around the mode (the consensus tree).\n\n> Some misprints: \"Monte-Calro\" figure text of Figure 2 line 4.\n\nThank you for your suggestion. We will correct the spelling."}, {"review_id": "tGtk4skfjA", "reviewer": "Reviewer_xhDC", "summary": "The authors propose a method for learning phylogenetic trees from sequence data. \n\nThe key idea is borrowed from [16], which is to represent the tree topology in terms of an embedding of the leafs of the trees in a continuous space, $z \\in \\mathcal{Z}$ from which the topology is extracted via a mapping $\\tau \\;:\\; \\mathcal{Z} \\to \\mathcal{T}$, where $\\mathcal{T}$ is the space of binary trees, using a distance-based technique (in particular the popular Neighbor-Joining (NJ) method). The fact that the representation in terms of $z$ is continuous has benefits in a variational representation of the posterior distribution over trees.\n\nThe theoretical contribution of the present paper is resolving some issues with \"an unevaluated Jacobian determinant between $B_\\tau$ and $z$\" (see Sec. 4 \"Related work\", p. 6).\n\nThe resulting method is demonstrated through simulations, showing that the method yields high-scoring models (high log-likelihoods), but not as high as existing variational and MCMC techniques (Table 2).", "questions": "My main question is what are the benefits compared to, e.g., MrBayes, VBPI-GNN, and the method proposed in [33], the first two of which are referred to as gold standard. I'm guessing the proposed method is computationally more efficient and scalable, but this should be clearly demonstrated (or am I missing something?).\n\nI should point out that I'm no expert of variational Bayes, so I didn't check the derivations from that point of view. I'm counting on the other reviewers to comment on that side.\n\ndetailed comments:\np. 7: \"marginal log-likelihood (MLL) estimates\": Does this mean MLL values? I mean, \"maximum (log-) likelihood estimate\" usually refers to an estimate of a parameter that is obtained by likelihood maximization, so it'd be good to be unambiguous.\np. 7: typo \"Monte-Calro\"\n ", "limitations": "I think so (but see my question about the benefits wrt. existing methods above).", "rating": 3, "confidence": 2, "soundness": 4, "presentation": 3, "contribution": 2, "strengths": "- I like the idea (which really borrowed from [16], so not really a new idea) of the embedding in terms of continuous coordinates\n- solid theoretical derivation of the variational estimators\n- promising results\n- quite well written (although a bit dense)", "weaknesses": "- not clear what the benefits are compared to existing variational and MCMC methods (which are presented as the gold-standard)", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your constructive review and comments.\n\n> My main question is what are the benefits compared to, e.g., MrBayes, VBPI-GNN, and the method proposed in [33], the first two of which are referred to as gold standard. I'm guessing the proposed method is computationally more efficient and scalable, but this should be clearly demonstrated (or am I missing something?).\n\nThank you for your question. Our main contribution is the development of variational Bayesian inference methods for phylogenetics without the preselection of tree topologies, which VBPI-GNN requires.\nFor the evaluation of the computational efficiency, we have included the performance comparison in Fig. R3 First, where the runtimes of our methods are comparable to VBPI-GNN for the same number of likelihood evaluations.\nAlthough MrBayes is a highly optimized and practical implementation for MCMC-based Bayesian phylogenetic inference,\nthe extensibility, variational Bayesian methods is a promising alternative that facilitates a gradient-based optimization, evaluation of model evidence, such as $\\ln P(Y)$, based on a tractable approximate posterior distribution, and the extensibility of the algorithms to more complex and large models.\n\n> detailed comments: p. 7: \"marginal log-likelihood (MLL) estimates\": Does this mean MLL values? I mean, \"maximum (log-) likelihood estimate\" usually refers to an estimate of a parameter that is obtained by likelihood maximization, so it'd be good to be unambiguous. \n\nThank you for your question. In the context of the variational inference, the marginal log-likelihood (MLL) estimate is used as a quality metric of the posterior approximation $Q(\\tau, B_\\tau)$.\nThis is because the expectation value of MLL estimates $L^{(K)}[Q(\\tau, B_\\tau]$ is the lower bound of the true MLL values, which coincides when the $Q(\\tau, B_\\tau) = P(\\tau, B_\\tau | Y)$.\n\n> p. 7: typo \"Monte-Calro\"\n\nThank you for your suggestions. We will correct the spelling of our revised manuscript. \n\n"}, {"review_id": "snFwLQZfyK", "reviewer": "Reviewer_Wg7V", "summary": "This work presents a new robust and scalable method for inferring phylogenetic trees based on (variational) Bayesian inference.", "questions": "Good performance compared to alternatives is shown but is it possible to demonstrate more on practical relevance for this?", "limitations": "No source code.\n\nNo assessment of societal implications.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "Originality\n- the originality of this work is in providing a robust and rigorous solution to an important application problem where the application of Bayesian methodology has so far been relatively limited\n\nQuality\n- the work is well motivated and appropriately implemented including the relevant derivations and experimental benchmarking \n\nClarity\n- language is fluent, and intuitive illustrations support reading; the technical details are described in sufficient detail and complemented by intuitive explanation. Conclucions and limitations are clearly stated. \n\nSignificance\n- significant, general and timely application problem addressed\n- improved robustness, while maintaining scalability that allows practical application\n\n\n", "weaknesses": "Treatment of the benchmark data sets is limited compared to the potential of the method for real applications. This could be expanded to highlight the relevance of the work.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "> 1. Treatment of the benchmark data sets is limited compared to the potential of the method for real applications. This could be expanded to highlight the relevance of the work.\n\nThank you for your constructive comments.\nRegarding the improved treatment of the current experiments, we introduced an intuitive understanding and comparison of topology spread and consensus trees through visualization analysis, as shown in Fig. R1. Additionally, to gain a deeper understanding of the estimated topology distribution Q(τ), we analyzed the relationship between the marginal log-likelihood (MLL) estimates and the fidelity of the majority consensus tree obtained from $Q(\\tau)$ in Fig. R2.\n\n> 2. Good performance compared to alternatives is shown but is it possible to demonstrate more on practical relevance for this?\n\nThank you for your question. In Fig. R2, we confirm that the performance of the marginal log-likelihood (MLL) estimates aligns well with the quality of the consensus tree obtained from $Q(\\tau).\nWe believe that the estimation of evolutionary parameters and demographic history, which are more practical issues, will become possible in the future by expanding the simple phylogenetic tree model and deriving variational algorithms for it.\n\n> 3. No assessment of societal implications.\n\nThank you for your insightful comments.\nRegarding the societal implication of our work,\nwe recognize the importance of addressing this aspect.\nScalable phylogenetic inference methods, like the ones we propose, have the potential to greatly enhance our understanding of the evolution, origin, and spread mechanisms of viruses and bacteria. Such insights could have profound societal impacts, especially in the areas of public health and disease control.\nWe will include a discussion on the potential applications and the future perspectives related to these methods in the 'Limitations and Future work' section.\n\n> 4. No source code.\n\nWe have included our source code in the supplementary materials.\nWe will ensure that our updated manuscript provides clear references to it for more clarity."}]

Author rebuttal:
We thank all the reviewers for taking the time to provide thorough and insightful feedback on our manuscript.
Your constructive comments have greatly enhanced the quality of our work.
In response to the points raised, we have engaged in a detailed discussion on the expressivity, performance, and limitations of our study with Figures R1 to R3 in the attached PDF file.

Meta-review: Thank you for engaging with the reviewers in the discussion. The reviewers and I are mostly in agreement that this paper should be accepted, given the the interesting methodology and clear exposition. However, the reviewers expressed some concerns, and after reading and considering the discussion I echo a number of these:

- Reviewer AJKq noticed that the derivation of the lower bound is only correct when $R(z|\tau)$ generates draws that respect the constraints of the tree construction. More technically, we require that $R(z|\tau)$ *dominates* $Q(z|\tau)$ for $KL(Q||R)$ to be well defined. The problem is that $R(z|\tau)$ is a distribution with full support on $z$, while $Q(z|\tau)$ only has support on a degenerate space where $\tau = \tau(z)$. Based on the authors' responses, I was not sure whether this issue was fully appreciated. However, at the end of the math display after line 141, the final expression just involves the ratio of $Q(z)$ and $R(z | \tau)$, which should be fine (both likely have full support, and are continuous densities without atoms). The authors should ideally find a proof of proposition 1 that doesn't involve a dubious step with the (probably) ill-defined KL. In the worst-case, the authors should clearly point out this limitation.
- The paper should be clear about its connection to [16] and standard variational lower bounds / stochastic gradient optimization.
- The effect of the control variates in Table 1 and 2 is somewhat inconclusive. The variation of results with respect to the various control variates is small compared to its performance with respect to CSMC, VCSMC, and $\phi$-CSMC. It would probably clarify the paper to just pick one any of these and defer this comparison to the appendix.
- The experiments section should provide an unbiased and accurate assessment of the results. In particular, the comparison to VBPI (and MrBayes) is underwhelming. Both VBPI and MrBayes generally outperform the method in terms of marginal log-likelihood bound. The authors argue that VBPI requires tree pre-selection, which is true, but the argument is not very strong because VBPI uses automated methods for the pre-selection step.
- The experiments in the submitted draft only show MLL comparisons, but not anything more particular (recovery of MCMC trees, for example). The authors did address this issue in the rebuttal with RF comparisons, but note that these results seem to show that VBPI/MrBayes are again competitive. The camera ready paper should include an unbiased assessment of these results.

Please make sure to explicitly address these issues in the camera ready.
Decision: Accept (poster)