PhyloGFN: Phylogenetic inference with generative flow networks
Authors: Ming Yang Zhou
Zichao Yan
Elliot Layne
Nikolay Malkin
Dinghuai Zhang
Moksh Jain
Mathieu Blanchette
Yoshua Bengio
Venue: ICLR 2024 poster (2024)
Keywords: Phylogenetic Inference
GFlowNets
Bayesian Inference
Deep Generative Modeling
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We use generative flow networks as amortized samplers of phylogenetic trees, achieving strong results in both Bayesian and parsimony-based phylogenetic inference.

Abstract: Phylogenetics is a branch of computational biology that studies the evolutionary relationships among biological entities. Its long history and numerous applications notwithstanding, inference of phylogenetic trees from sequence data remains challenging: the high complexity of tree space poses a significant obstacle for the current combinatorial and probabilistic techniques. In this paper, we adopt the framework of generative flow networks (GFlowNets) to tackle two core problems in phylogenetics: parsimony-based and Bayesian phylogenetic inference. Because GFlowNets are well-suited for sampling complex combinatorial structures, they are a natural choice for exploring and sampling from the multimodal posterior distribution over tree topologies and evolutionary distances. We demonstrate that our amortized posterior sampler, PhyloGFN, produces diverse and high-quality evolutionary hypotheses on real benchmark datasets. PhyloGFN is competitive with prior works in marginal likelihood estimation and achieves a closer fit to the target distribution than state-of-the-art variational inference methods.

Published as a conference paper at ICLR 2024
PHYLOGFN: PHYLOGENETIC INFERENCE WITH GENERATIVE FLOW NETWORKS
Mingyang Zhou' McGill University
Zichao Yan
Elliot Layne
McGill University Université de Montréal, Mila
McGill University, Mila
Nikolay Malkin Université de Montréal, Mila
Dinghuai Zhang Université de Montréal, Mila
Moksh Jain
Université de Montréal, Mila
Mathieu Blanchette McGill University, Mila
Yoshua Bengio Université de Montréal, Mila
[ming. zhcu, zichao. yan, ellict.......... 1ayneiem.ayeeemmill..... wwl.lmmmmmll... blanchemêcs .magiiise.mmmmimmmmoo. 11 ca (nikolay coolay.malkee malkin, dinghuai. zhang, moksh. jain, yosiaa.eenooaeeeee.......
ben.iio quebbe
ABSTRACT
Phylogenetics is a branch of computational biology that studies the evolutionary relationships among biological entities. Its long history and numerous applica- tions notwithstanding, inference of phylogenetic trees from sequence data remains challenging: the extremely large space poses significant obstacle for the cur- rent combinatorial and probabilistic techniques. In this paper, we adopt the frame- work of generative flow networks (GFlowNets) tackle two problems in phylogenetics: parsimony-based and Bayesian phylogenetic inference. Because GelowNets are well-suited for sampling complex combinatorial structures, they are a natural choice for exploring and sampling from the multimodal posterior distribution tree topologies and evolutionary distances. We demonstrate that our amortized posterior sampler, PhyloGFN, produces diverse and high-quality evolutionary hypotheses real benchmark datasets. PhyloGFN is competitive with prior works in marginal likelihood estimation and achieves a closer fit to the target distribution than state-of-the-ar variational inference methods.
INTRODUCTION
Phylogenetic inference has long been central problem in the field of computational biology. Ac- curate phylogenetis inference is critical for a number of important biological analyses, such as un- derstanding the development antibiotic resistance (Aminov & Mackic 2007 Ranjbar al. 2020
Layne et al. 2020). assessing the risk of invasive species (Hamelin et al. 2022, Dort et al. 2023
and many others. Accurate phylogenetic trees can also be used improve downstream computa- tional analyses, such as multiple genome alignment (Blanchette et al. 2004), ancestral sequence reconstruction (Ma et al., 2006), protein structure and function annotation (Celniker al. 2013).
Despite its strong medical relevance and wide applications in life science, phylogenetic inference has remained a standing challenge, in part due to the high complexity of tree space - for species, (2n - 5)!! unique unrooted bifurcating branches of phylogenetic inference; both maximum-likelitunn
topologies exist. This poses common obstacle to all
and maximum-parsimony tree recon- struction NP-hard problems (Day) 1987 Chor & Tuller 2005). Under the Bayesian formulation of phylogenetics, the inference problem further compounded by the inclusion of continuous vari- ables that capture the level of sequence divergence along each branch of the
One line of prior work considers Markov chain Monte Carlo (MCMC)-based approaches, such as MrBayes (Ronquist et al.] 2012). These approaches have been successfully applied to Bayesian phylogenetic inference. However, a known limitation of MCMC is scalability to digh-dimeensionl
Equul contribution. *CIFAR Senior Fellow.

---

Published as a conference paper at ICLR 2024
distributions with multiple separated modes (Tjelmeland & Hegstad 2001), which arise in larger phylogenetic datasets. Recently, variational inference (Vi)-based approaches have emerged. Among these methods, some model only limited portion of the space of tree topologies, while others are weaker in marginal likelihood estimation due to simplifying assumptions. In parsimony analysis, dfateoof-the-ar methods such as PAUP* (Swofford (1988 have extensively relied on heuristic search algorithms that are efficient but lack theoretical foundations and a guarantees.
Coming from the intersection of variational inference and reinforcement learning is the class of mod- els known generative flow networks (GFlowNets: Bengio et al. (021). The flexibility afforded by GFlowNets to learn sequential samplers for distributions over compositional objects makes them a promising candidate for performing inference over the posterior space of phylogenetic tree topolo- gies and evolutionary distances.
In this work, we propose PhyloGFN, the first adaptation of GFlowNets to the task of Bayesian and parsimony-based phylogenetic inference. Our contributions are as follows:
(I) We design an acyclic Markov decision process (MDP) with fully customizable reward functions, by which our PhyloGFN can be trained to construct phylogenetic trees bottom-up fashion. (2) PhyloGFN leverages novel tree representation inspired by Fitch and Felsenstein's 's algorithms to represent rooted trees without introducing additional learnable parameters to the model. Phy- loGFN is also coupled with simple yet effective training techniques such using mixed on- policy and fithered-policy rollouts, replay buffers and cascading cmmpeaature-anneaiing
(3) PhyloGFN has the capacity to explore and sample from the entire phylogenetic tree space, achieving a balance between exploration in this vast space and high-fidelity modeling of the modes. While PhyloGFN performs on par with the state-of-the-art MCMC- and VI-based meth- ods in the summary metric of marginal log-likelihood, it substantially outperforms these ap- proaches in terms of its ability to estimate the posterior probability of suboptimal trees.
2 RELATED WORK
Markov chain Monte Carlo (MCMC)-based algorithms are commonly employed for Bayesian phy- logeneties, with notable examples including MrBayes and RevBayes (Ronquist et al. 2012 Höhna
ct al. 2016). which are considered state-of-the-art in the field. Amortized variational inference (VI) alternative approach that parametrically estimates the posterior distribution. VBPI-
GNN Zhhnng 2023) employs subsplit Bayesian networks (SBN) (Zhang & Matsen IV (0118)) to
model tree topology distributions and uses graph neural networks to learn tree topological embed- dings (Zhang (023). While VBPI-GNN has obtained marginal log likelihood competitive with MrBayes in real datasets, it requires pre-generated set of high-quality tree topologies constrain its action space for tree construction, which ultimately limits its ability to model the entire tree space.
There exist other VI approaches that do not limit the space of trees. VaiPhy (Koptagel et al. 2022) approximates the posterior distribution in the augmented space of tree topologies, edge lengths, and ancestral sequences. Combined with combinatorial sequential Monte Carlo (CSMC; Moretti et al., 2021). the proposed method enables faster estimation of marginal likelihood. GeoPhy (Mimori & Hamada (023) models the tree topology distribution in continuous space by mapping continuous- valued coordinates to tree topologies, using the same technique as VBPI-GNN to model tree topo- logical embeddings. While both methods model the entire tree topology space, their performance on marginal likelihood estimation underperforms the state of the art.
For the optimization problem underpinning maximum parsimony inference, PAUP* is one of the most commonly used programs (Swofford. 1998); it features several fast, greedy, and heuristic al- gorithms based local branch-swapping operations such as tree bisection and reconnection.
GFlowNets a family of methods for sampling discrete objects from multimodal distributions, such as molecules (Bengio et al. 2021) and biological sequences (Jain et al.) (022), and are used to solve discrete optimization tasks (Zhang al. 2023a b). With their theoretical foundations laid out
in Bengio et al. (2023); Lahlou et al. (2023) and connections to variational inference established in Malkin et al. (2023). GFlowNets have been successfully applied to tackle complex Bayesian inference problems, such as inferring latent causal structures in gene regulatory networks (Deleu
al. 2022 2023), and parse trees in hierarchical grammars (Hu al. 2023).
2

---

Published as a conference paper at ICLR 2024
3 BACKGROUND
3.1 PHYLOGENETICC INFERENCE
Here we introduce the problems of Bayesian and parsimony-based phylogenetic inference. A weighted phylogenetic tree is denoted by (2, where represents the tree topology with its leaves labeled by observed sequences, and represents the branch lengths. The tree topology can be either rooted binary tree bifurcating unrooted tree. For tree topology 2, let L(z) denote the labeled sequence set and (zz) the set of edges. For an edge E E(z), let b(e) denote its length. Let Y 22..)....Hu E Exxm be set of n observed sequences, each having m characters from alphabet Σ, e.g. GA, for DNA sequences. We denote the ش site of all sequences by y' = [i]. ... V,[]]. In this work, make two assumptions that are common in the phy- logenetic inference literature: (i) sites evolve independently; (ii) evolution follows time-reversible substitution model. The latter implies that unrooted tree has the same parsimony score or like- lihood as its rooted versions, and thus the algorithms we introduce below (Fitch and Felsenstein) apply unrooted trees well.
3.1.1 BAYESIAN INFERENCE
In Bayesian phylogenetic inference, we are interested sampling from the posterior distribution over weighted phylogenetic trees (z, b), formulated as:
P(z. Y) P(z. P(Y z)
P(Y)
where P(Y b) is the likelihood, P(Y) is the intractable marginal, and P(z, b) is the prior density over tree topology and branch lengths. Under the site independence assumption, the likelihood can be factorized as: P(Y) z, b) П., P(Y' b), and each factor is obtained by marginalizing over all internal nodes a and their possible character assignment:
2n-2
n
Σ S-
P(azx-)
I jun+1
Pldylaaun b(e,)) I
Plralilldaay be(e)) km1
where a' "xt'".- al 2n- represent the internal node characters assigned to site and a(i) represent the parent of node P(a2--1) is distribution at the root node, which is usually assumed to be uniform over the vocabulary, while the conditional probability P(atiaa' agy b(e,)) is defined by the substitution model (where e is the edge linking a to a(a))).
Felsenstein's 's algorithm The likelihood of a given weighted phylogenetic tree can be calculated efficiently using Felsenstein's pruning algorithm (Felsenstein) 1973) in bottom-up fashion through dynamic programming. Defining 1,, as the leaf sequence characters at site i below the internal node u, and given its two child nodes and w, the conditional probability P(LÍ, laf,) can be obtained from P(L(la() and P(L ):
P(L d',))
Σ ab.cker
P(a I b(e.))) d.)((a, a b(e.)) d.,). (1)
The dynamic programming. or recursion, essentially a post-order traversal of the tree and P(L,, a')) is calculated at every internal node u, and we use one-hot encoding of the sequences to represent the conditional probabilities at the leaves. Finally, the conditional probability for each node at site i is stored in a data structure fú E [0. || fúlc] = P(L,la" = c). and we call it the Felsenstein feature for node u. Note that the conditional probability at the root P(Y (ain-1) P(Y
is used
to calculate the likelihood of the tree: P(Y' b) Edi,,,e ((
3.1.2 PARSIMONY ANALYSIS
The problem of finding the optimal tree topology under the maximum parsimony principle is com- monly referred as the Large Parsimony problem, which harhard. For a given tree topology z, the parsimony score is the minimum number of character changes between sequences across branches obtained by optimally assigning sequences
internal nodes. Let M (z(Y) be the parsimony score of tree topology 2 with leaf labels Y. Due to site independence, M(z = Σ, (z(Y'). The trees with the lowest parsimony score, or most parsimonious trees, are solutions the Large Parsimony problem. Note that the Large Parsimony problem is limiting case of the maximum likelihood tree inference problem, where branch lengths are constrained be equal and infinitesimally short.
3

---

Published as a conference paper at ICLR 2024
Fitch algorithm Given rooted tree topology the Fitch algorithm assigns optimal sequences to internal nodes and computes the parsimony score in linear time. At each node u, the algorithm tracks the set of possible characters labeling for node that can yield most parsimonious solution for the subtree rooted at u. This character set can be represented by a binary vector fü E {0. 1) | for site We label this vector the Fitch feature. As in Felsenstein's 's algorithm, this vector a one-hot encoding of the sequences at the leaves and is computed recursively for non-leaves. Specifically, given a rooted tree with root and two child trees with roots V and w, fú is calculated as:
S A fi S V f otherwise
if ft f #
fu =
where A and V are element-wise conjunctions and disjunctions. The algorithm traverses the tree two times, first in post-order (bottom-up) to calculate the character set each node, then in pre- order (top-down) to assign optimal sequences. The total number of character changes between these optimal sequences along the tree's edges counted as the parsimony score.
3.2 GELOWNETS
Generative flow networks (GFlow/eets) are algorithms for learning generative models of complex distributions given by unnormalized density functions over structured spaces. Here, we give a con- cise summary of the the GFlowNet framework.
Setting A GelowNet treats generation of objects x lying in a sample space X as a sequential makini problem an acyclic deterministic MDP with set of states S X and set of actions R C S x S. The MDP has a designated initial state so. which has no incoming actions, and set of terminal states (those with no outgoing actions) that coincides with X. Any reached from by sequence of actions so - - ... - = (with each (s,, 52+1) E A). Such sequences are called complete trajectories, and the of all complete trajectories is denoted T.
X can be
A (forward) policy PF is a collection of distributions Pp(s') s) over the children of each nontermi- nal state s E S \ X. A policy induces a distribution over T:
PF(T (sD
s,))
[ PF(S st).
1=0
A policy gives a way to sample objects in X, by sampling a complete trajectory t ~ PF and returning its final state, inducing a marginal distribution P over X: PJ is the sum of PP(T) over all complete trajectories r that end in x (a possibly intractable sum).
The goal of GelowNet training is to estimate a parametric policy PF(- | :; such that the induced P proportional a given reward function R X Eso, i.e.
(x R(x) V
where Z = Exex R(x) is the unknown normalization constant (partition function).
Trajectory balance objective The direct optimization of PF's parameters 0 is impossible since it involves an intractable sum over all complete trajectorics. Instead, we leverage the trajectory balance (TB) training objective (Malkin et al. (002),, which introduces two auxiliary objects: an estimate Ze of the partition function and a backward policy. In our experiments, we fix the backward policy to uniform, which results in simplified objective:
LTB(T) = log Zo 11/50 Pr(S 1;;0))
[x) =
(3)
R(x)Po(T (x)
PP((s))
where Pa(s) denotes the set of parents of s. By the results of Malkin et al.; (2022). there exists a unique policy PF and scalar Zo that simultaneously make LTB(T) = 0 for all r E T, and at this optimum, the policy Pr satisfies (2) and Zo equals the true partition function Z. In practice, the policy Pr(- s; 0) is parameterized as a neural network that outputs the logits of actions - s' given a representation of the state y as input, while Zo is parameterized in the log domain.
(TB(T) is minimized by gradient descent on trajectories 7 chosen from a behaviour policy that can encourage exploration to accelerate mode discovery (see our behaviour policy choices in $4.3)
4

---

Published as a conference paper at ICLR 2024
Crestee
Drexding
Partria
Med
Traelig
simplry
Figure I: Left: PhyloGFN's state space on a four-sequence dataset. Initial state so comprises leaf nodes. Successive steps merge pairs of trees until single unrooted tree remains. Right: Policy model for PhyloGFN- Bayesian. Transformer encoder processes tree-level features xi ((2),b)) b1),bb))... ... (z). by)). Pairwise features eij are derived and used by MLPs seleet tree pairs for merging and sample branch lengths.
4 PHYLOGENETIC INFERENCE WITH GFTSOWNETS
4.1 GFLOWNETS FOR BAYESIAN PHYLOGENETIC INFERENCE
This section introduces PhyloGFN-Bayesann netic inference. Given a set of observed sequences Y, PhyloGFN-Bayeste
GFlowNet-based method for Bayesian phyloge- learns a sampler over the joint space of tree topologies and edge lengths X = {(2,( (z) = Y) such that the sampling
probability of (z, b) E X approximates its posterior (z, b) P(z. blY).
We follow the same setup as Koptagel et al. (2022); Zhang 2023): : Mimori & Hamada (202)): (i)
uniform prior over tree topologies; (ii) decomposed proor P(z. b) P(2) (ii) exponential
(t= 10) prior over branch lengths; (iv) Jukes-Cantor substitution model.
GFlowNet state and action space The sequential procedure of constructing phylogenetic trees is illustrated in Fig. [] The initial state is set of rooted trees, each containing single leaf node labeled with an observed sequence. Each action chooses a pair of trees and joins their roots by a common parent node, thus creating new tree. The number of rooted trees in the set is reduced by at every step, so after n - steps a single rooted tree with n leaves is left. To obtain an unrooted tree, we simply remove the root node.
Thus, state consists of set of disjoint rooted trees U, L(z,) = Y. Given a nonterminal state with > I trees, a transition action consists of two steps: (i) choosing a pair of trees to join out of the (5) possible pairs; and (ii) generating branch lengths for the two introduced edges between the new root and its two children. The distribution over the pair of branch lengths is modeled jointly as a discrete distribution with fixed bin size. Following the initial submission of the paper, we conducted additional experiments by employing continuous distribution to model branch lengths. Further details can be found in of the appendix.
(1z,,b))...... bi)), I S and
Reward function We define the reward function as the product of the likelihood and the edge length prior: R(z, b) = P(Y|z, b) P(b), implicitly imposing a uniform prior over tree topologies. By training with this reward, PhyloGFN learns to approximate the posterior, since P(z, b|Y) R(z. b) P(z) and P(z), P(Y) are both constant.
P(Y
It is worth emphasizing that in our bottom-up construction of trees, the set of possible actions the steps that select two trees to join by new common root is never larger than n² even though the size of the space of all tree topologies- all of which can be reached by our sampler - superexponential in n. This stands in contrast to the modeling choices of VBPI-GNN (Zhang 2023), which constructs trees in top-down fashion and limits the action space using a pre-generated set of trees, therefore also limiting the set of trees it can sample.
State representation To represent a rooted tree in a non-terminal state, we compute features for each site independently by taking advantage of the Felsenstein features (33.1.)) Let (z, b) be weighted tree with root which has two children and w. Let f f(. ft, E [0. [E] be the Felsen-
5

---

Published as a conference paper at ICLR 2024
stein feature nodes u, v, at site i. The representation p for site i is computed as following:
Pu fu
11 ((b=(e))
(4)
cEE(Z)
where (e) llebb P(b(e)) is the edge length prior. The tree-level feature is the concatena- tion of site-level features p= [p''...p'' p... A state s (z)... which a collection of rooted trees, is represented by the set pr).
Representation power Although the proposed feature representation does not capture all the information of tree structure and leaf sequences, we show that p indeed contains sufficient informa- tion to express the optimal policy. The proposition below shows that given optimal GelowNet with uniform PB. two states with identical features set share the same transition probabilities.
Proposition 1. Let s1 {(11,))) (22, (o))...EE1 bI)) and $2 {Z2, b1), (2/2, b2) ... ... (2G.b)) b)) be two non-terminal states such that $ # 82 but sharing the same features pi = P. Let a be any sequence of actions, which applied 81 and 82, respectively, results in full weighted frees x (z, b2,x),,' = (z',b') with two partial trajectories T == (S1 PF is the policy of optimal GFlowNet with uniform PB, then P Pr( PF(T').
_ xx. = (82 - ... - x'). If
All proofs are Appendix A The proposition shows that our proposed features have sufficient representation power for the edge length priors are used calculating reward by Felsenstein's 's algorithm. Therefore, computing these features does not introduce any additional variables, and computation overhead is minimized.
policy. Furthermore, Felsenstein features and
4.2 GELSWNETS FOR PARSIMONY ANALYSIS
This section introduces PhyloGFN PhyloGFN-Parsimont Parsimony, our GFlowNet- method for parsimony analysis. We treat large parsimony analysis. a discrete optimization problem, as a sampling problem from the
(2 defined over tree topologies. Here, M (z(Y) is the parsimony energy distribution exp T
score of and T is pre-defined temperature term control the smoothness of distribution. With sufficiently small T, the most parsimonious trees dominate the energy distribution. To state our goals learns a sampling policy Pr over the Migh0
formally, given observed sequences PhyloGFN-Pariimon
space of tree topologies (z)((z) = Y} such that (z(2) & e 7 As - O, this target distribution approaches uniform distribution over the set of tree topologies with minimum parsimony scores.
PhyloGFN-Parsing
can be seen as reduced version of SyyloGFN-Bayesiinn generation procedure is the same as before, but we no longer generate branch lengths. The reward defined as R(2) = exp C-M(zz)) where C is an extra hyperparameter introduced for stability to T
The tree shape
offsct the typically large (=|Y) values. Note that C can be absorbed into the partition function and has no influence on the reward distribution.
Similar to MaoooGGN--Byyeiiinn a state (collection of rooted trees) is represented by the set of tree features, with each tree represented by concatenating its site-level features. With z the rooted tree topology with root u, we represent the tree at site i by its root level Fitch feature fú defined in $3.1.2 The proposition below, analogous to Proposition [ shows the representation power of the proposed feature.
Proposition Let 81 (21, 22........ 2) and 82 (2 2.....z/ be two non-terminal states such that # 82 but sharing the same Fitch features fz. = F2 Vi. Let be any sequence of actions, which, applied to and respectively, results in tree topologies x,x' E Z. with two partial trajectories (s)
x' (S2 3 - x'). If Pr the policy of an optimal GFlowNet with uniform Pn. then PT( Pr(r')
This shows that the Fitch features contain sufficient information for DoyyloGFN-Parsimonn Further- more, the Fitch features are used the computation of the reward by Fitch's 's algorithm, so their use in the policy model does not introduce additional variables or extra computation.
PhyloGEN The temperature T controls the trade-off between sample diversity and parsimony scores. Following Zhang et al. (2023a), we extend PhyloGFN-Parsimont
cemperature-condiiionn
by conditioning the policy on T, with reward R(z:T) exp
C (-M(c))) and we learn sampler
T
such that PT (T;;)) R)z;T). See Appendix for more details.

---

Published as a conference paper at ICLR 2024
Table 1: Marginal log-lik likelihhodd estimation with different methods on real datasets DSI-DS8. PhyloGFN
outperforms d-CSMC across all datasets and GeoPhy most. *VBPI-GNN uses predefined tree topologies in training and is directly comparable.
MCMC
ML-based / amortized, full tree space
Dataset
MrBayes SS
VBPI-GNN
<-CSMC
GeoPhy
PhyloGFN
DSI
-7108.42
=0.18
-7108.41
+0.14
-7290.36
+7.23
-7111.55 +0.07
-7108.95 +0.06
DS2
-26367.57
10.48
-26367.73
10.07
-30568.49
131.34
-26368 44 10.13
-26368 90 10.28
DS3
-33735.44
20.5
-33735.12
10.09
-33798.00
16.62
-33735.88
10.12
-33735. 6 10.35
DS4
-13330,06
20.54
-13329,94
10.19
-13582.24
135.08
-13337.42
11:32
-13331.3 .83 10.19
DS5
-8214.51
=0.28
-8214.64
+0.38
-8367.51
+8.87
-8233.89
+6.63
-8215. +0.2
DS6
-6724.07
=0.86
-6724.37
=0.4
-7013.83
+16.99
-6733.91
+0.57
-6700.6 +0.54
DS7
-37332.76
+2.42
37332.04
+0.12
-37350.77
+11.74
-37359.99
+1.14
DS8
-8649.88 11.75
-8650.65
20.45
-9209.18
118.03
-8660.48
10.78
-8654.76 10.19
4.3 MODEL ARCHITECTURE AND TRAINING
Parameterization of forward transitions We parameterize the forward transitions of tree topol- ogy construction using Transformer-based neural network, whose architecture is shown in Fig. T
We select Transformer because the input is a set and the model needs to be order-equivariant For state consisting of trees, after n embeddings are generated from the Transformer encoder, (5) pairwise features are created for all possible pairs of trees, and a common MI.P generates probability logits for joining every tree pair. ByyloGFN-Bayesinnn additionally requires generating edge lengths.
Once the pair of trees to join is selected, another MI.P is applied to the corresponding pair feature to generate probability logits for sampling the edge lengths. See more details in D
Off-policy training The action model P s; 0) is trained with the trajectory balance objective. Training data are generated from two sources: (i) A set of trajectories constructed from the currently trained GFlowNet, with actions sampled from the policy with probability 1 and uniformly at random with probability e. The e rate drops from a pre-defined Extars to near 0 during the course of training. (ii) Trajectories corresponding to the best trees are sampled backward from these high-reward trees with uniform backward policy.
date (replay buffer). Trajectories
Temperature annealing For PhyloGFN-Parssee Parsimony, it is crucial to choose the appropriate tem- perature T. Large 7 defines flat target distribution, while small 7 makes the reward landscape less smooth and leads to training difficulties. We cannot predetermine the ideal choice of 7 before training. do not know priori the parsimony score for the dataset. Therefore, initialize the training with large T and reduce periodically during training. See Appendix for details.
4.4 MARGINAL LOG-LIKELIHOOM ESTIMATION
To assess how well the GFlowNet sampler approximates the true posterior distribution, the following importance-weightee variational lower bound on the marginal log-likelihood (MLL):
log P(Y) 2 E TPP log P(z)
Σ
(5)
PF(T))
T2150 (.50)
Our estimator is computed by sampling batch of K trajectories and averaging P(a))P(rlE. PP(T) over the batch. This expectation of this estimate is guaranteed to be a lower bound on log P(Y) and its bias decreases as K grows (Burda et al. (016).
models branch lengths with discrete multinomial distributions, while in reality these are continuous variables. To properly estimate the MLL and compare with other methods de- fined over continuous space, we augment our model to a continuous sampler by performing random perturbations over edges of trees sampled from PhyloGFN-Bayesian. The perturbation follows the uniform distribution U|50.5..0..5.) where w is the fixed bin size for edge modeling in PhyloGFN- Bayesian. The resulting model over branch lengths then piecewise-constant continuous distri- bution. We discuss the computation details as well as the derivation of (5) in Appendix
PhyloGFN-
5 EXPERIMENTS
We evaluate PhyloGFN on a suite of real-world benchmark datasets (Table ST in Appendix C that standard in the literature. These datasets feature pre-aligned sequences and vary in difficulty (27

---

Published as a conference paper at ICLR 2024
to 64 sequences; 378 to 2520 sites). In the following sections we present results and analysis on Bayesian and parsimony- based phylogenetic inference.
5.1 BAYESIAN PHYLOGENETIC INFERENCE
CPN aaddm brees
VB7% 221 0% modom steee
PhyloGFN is compared with variety of baselines in
terms of sampling-bass estimates of marginal log- likelihood (MLL; see details in $4.4). The baselines we compare to are MrBayes SS, Stepping-Stone sampling algorithm implemented in MrBayes (Ron- guist et al. 2012). and three variational inference methods: VBPI-GNN (Zhang) 2023). 6-CSMC pro-
-4370
4475
-500
4335
210
20
301 nencee
24h
201
VSP 10%
218
20
291
-0000
-4300-
posed in VaiPhy (Koptagel et al.) (022). and Geo- Phy (Mimori& & Hamada 2023). The sampling setup for MrBayes follows Zhang & Matsen IV (2018b)
-000
8150-
3300
and otherwise show the highest MLL reported in the
-3000
respective papers. PhyloGFN MLL is estimated fol-
1i4 H5
212
265
353
lowing the formulation in $4.4. with mean and stan-
6900
dard deviation obtained from 10 repetitions, each us-
ing 1000 samples. See Sections B and G for addi- tional training details, hardware specifications, and running time comparison. Additional results from two repeated experiments are in Table S5
2000
3330
-200-
Sandingg Daiicly
Sanuiing Bliidyy
The results are summarized in Table S4 Our Phy- loGEN is markedly superior to Φ-CS 6-CSMC across all datasets and outperforms GeoPhy on most, with the exception of DS2 and DS3 where the two perform similarly, and DS7 where GeoPhy obtains a better result. VBPI-GNN, is the only machine learning- based method that performs on par against MrBayes, the current gold standard in Bayesian phylogenetic inference. However, it should be emphasized that VBPI-GNN requires a set of pre-defined tree topologies that are likely to achieve high likelihood, and as a consequence, its training and inference are both constrained to a small space tree topologies. On the other hand, PhyloGFN operates on the full space of tree topologies and, in fact, achieves a closer fit to the true posterior distribution. To show this, for each dataset, we created three sets of phylogenetic trees with high/medium//oww posterior probabilities and obtained the corresponding sampling probabilities from PhyloGFN and VBPI-GNN. The three classes of trees are generated from VBPI-GNN by randomly inserting uni- formly sampled actions into its sequential tree topology construction process with 0%c 30%, 50% probability, respectively. which cireumvents VBPI-GNN's limitation of being confined to a small tree topology space. Table 2 and Fig. show that PhyloGFN achieves higher Pearson correlation between the sampling log- log-probability and the unnormalized ground truth posterior log-density for the majority of datasets and classes of trees. In particular, while VBPI-GNN performs better on high-probability trees, its correlation drops significantly on lower-probability trees. On the other hand, PhyloGFN maintains a high correlation for all three classes of trees across all datasets, the only exception being the high-probabiliiyy trees in DS7. See Appendix K for details and extended results and Appendix for short explanation of the significance of modeling suboptimal trees.
Figure 2: Model sampling log log-density vs. unnormalized posterior log-density for high/meddium//owprrbbbbiil We highlight that PhyloGPN-Bayesian per- forms significantly better on medium- and low-probability trees, highlighting its superiority in modeling the entire data space.
trees on DS1.
Continuous branch length modeling Following the initial submission, additional experiments involving continuous branch length modeling demonstrate PhyloGEN's ability to achieve state-of- the-art Bayesian inference performance. For more details, please see Appendix F
5.2 PARSIMONY-BASEE PHYLOGENETIC INFERENCE
As a special case of Bayesian phylogenetic inference, the parsimony problem is concerned with finding the post-paassmooious trees task which is also amenable to PhyloGFN. Here, we compare to the state-of of-te-of-the-art parsimony analysis software PAUP* (Swofford), (998). For all datasets, our PhyloGFN and PAUP* are able to identify the same set of optimal solutions to the Large Parsimony problem, ranging from a single optimal tree for DSI to 21 optimal trees for DS8.
Although the results are similar between PhyloGFN and PAUP*, once again we emphasize that PhyloGFN is based on rigorous mathematical framework of fitting and sampling from well-defined
8

---

Published as conference paper at ICLR 2024
Table 2: Pearson correlation of model sampling log-density and ground truth unnormalized posterior log- density for each dataset high/mediumm/oww posterior density trees generated by VBPI-GNN. PhyloGFN achieves good fit on both high and low posterior density regions.
No random
30% random
50% random
Dataset
PhyloGFN
VBPI-GNN
PhyloGFN
VBPI-GNN
PhyloGFN
VBPI-GNN
DS1
0.994
0.955
0.961
0.589
0.955
0,512
DS2
0.930
0.952
0.948
0.580
0.919
0,538
DS3
0.917
0.968
0.963
0.543
0.950
0.499
DS4
0.942
0.960
0.945
0.770
0.966
0.76
DS5
0.969
0,965
0.937
0.786
0.939
0.773
DS6
0.993
0.887
0.973
0,816
0.934
0,702
DS7
0.624
0.955
0.787
0,682
0.764
0,678
DS8
0.978
0,955
0.913
0.604
0,901
0.463
posterior distributions over tree topologies, whereas PAUP relies on heuristic algorithms. To put it more concretely, we show in Fig. S2 that PhyloGFN is able to (i) learn a smooth echelon of sampling probabilities that distinguish the optimal trees from suboptimal ones; (ii) learn similar sampling probabilities for trees within each group of parsimonious trees; and (iii) fit all 3 rooted that belong to the same unrooted tree equally well.
Finally, Fig. 3 shows that single temperature- conditioned PhyloGFN can sample phylogenetic trees of varying ranges parsimony scores by pro- viding suitable temperatures 7 input to the model. Also, PhyloGFN is able to sample proportionately from the Boltzmann distribution defined at differ- temperatures and achieves high correlation be- tween sampling log- log-probability and log-reward. Al- though the temperature-cone innperaaure-coonditioeedd PhyloGFN has only been trained on small range of temperatures between 4.0 and 1.0, Fig. Fig.33ssoo 3/sows it can also approx- imate the distribution defined at temperature 8.0. Further results are presented in Appendix L
1200
owket T=B.C
1000-
samp GF owhet samp es T=10
800-
COUPS2
4350
parsimony Score
No02 Mogg0ti BUT mont scure
B
-4050
41.00
4150
DISCUSSION AND FUTURE WORK
4200
4200
4150
4250
In this paper, propose PhyloGFN, GFlowNet- based generative modeling algorithm, to solve parsimony-based and Bayesian phylogenetic in- ference. We design an intuitive yet effective tree construction procedure to efficiently model the entire tree topology space. We propose a novel tree representation based Fitch's and Felsen- stein's 's algorithms to represent rooted trees without introducing additional learnable parameters, and we show the sufficiency of our features for the purpose of phylogenetic tree generation. We apply our algorithm on eight real datasets, demonstrating that PhyloGFN is competitive with superior to prior works in terms of marginal likelihood estimation, while achieving closer fit to the target distribution compared to state-of-thh--ar variational inference methods. While initial experiments with continuous branch length sampling have demonstrated notable performance enhancements, there remains a need for future research to address training efficiency. In addition, we plan explore the of conditional GFlowNets amortize the dependence on the sequence dataset itself. This would allow a single trained GFlowNet to sample phylogenetic trees for sets of sequences that were not during training.
Pearsonk 0.84
Pearson 0.47
8 x
4350
Ghowwket probabiity Te8,0
Ghiowitee probabiity Tu4.0
100000 SA Aapeeeo scurs
-4500
4000
4100
4100
4150
4220
4200
4250
-3000
Pearsant: 0.555
4300
Pearsont:
-4350
probability
T=2.0
Te2.0
Figure 3: A temperature-conndiionee PhyloGFN is trained on DS1 using temperatures sampled be- tween 4.0 and 1.0. (A) Parsimony score distribu- tion with varying statistical temperature input (8.0 to 1.0) to PhyloGFN policy (10,000 trees sampled per temperature). (B) PhyloGFN achieves high Pearson correlation at each temperature.

---

Published as a conference paper at ICLR 2024
REPRODUCIBILITY
Code and data are available for detailed description.
See Appendix g
ACKNOWLEDGMENTS
We thank Oskar Kviman and Jens Lagergren for the insightful discussions on GFlowNet and phy- logenetic inference. We thank Tianyu Xie and Cheng Zhang for discussions on continuous branch lengths modeling. We acknowledge the support of Mingyang Zhou by the NSERC Discovery grant awarded Mathieu Blanchette, and the funding from CIFAR, NSERC, IBM, Intel, Genentech and Samsung to Yoshua Bengio. Finally, we thank the Digital Alliance of Canada for providing the computational resources.
REFERENCES
Rustam I Aminov and Roderick Mackie. Evolution and ecology of antibiotic resistance genes. FEMS microbiology letters, 271(2)147-161,, 2007.
Emmanuel Bengio, Moksh Jain, Maksym Korablyov, Doina Precup, and Yoshua Bengio. Flow net- work based generative models for non-iterative diverse candidate generation. Neural Information Processing Systems (NeurIPS), 2021.
Yoshua Bengio, Salem Lahlou, Tristan Deleu, Edward J Hu, Mo Tiwari, and Emmanuel Bengio. GFlowNet foundations. Journal of Machine Learning Research, (24):1-76, 2023.
Mathieu Blanchette, W James Kent, Cathy Riemer, Laura Elnitski, Arian FA Smit, Krishna M Roskin, Robert Baertsch, Kate Rosenbloom, Hiram Clawson, Eric D Green, et al. Aligning mul- tiple genomic sequences with the threaded blockset aligner. Genome research, 14(4):708-715, 2004.
Patricia Buendia, Brice Cadwallader, and Victor DeGruttola. phylogenetic and Markov model approach for the reconstruction of mutational pathways of drug resistance. Bioinformatics, 25 (19):2522-2529, 08 2009.
Yuri Burda, Roger Baker Grosse, and Ruslan Salakhutdinov. Importance weighted autoencoders. International Conference on Learning Representations (ICLR), 2016.
Gershon Celniker, Guy Nimrod, Haim Ashkenazy, Fabian Glaser, Eric Martz, Itay Mayrose, Tal Pupko, and Nir Ben-Tal. Consurf: using evolutionary data to raise testable hypotheses about protein function. Israel Journal of Chemistry, 53(3-4):199-26,, 2013.
Benny Chor and Tamir Tuller. Maximum likelihood of evolutionary trees is hard. In Annual Inter- national Conference on Research Computational Molecular Biology, pp. 296-310. Springer, 2005.
William HE Day. Computational complexity of inferring phylogenies from dissimilarity matrices. Bulletin of mathematical biology, 49(4):461-467, 1987.
Tristan Deleu, António Góis, Chris Emezue, Mansi Rankawat, Simon Lacoste-Julien, Stefan Bauer, and Yoshua Bengio. Bayesian structure learning with generative flow networks. Uncertainty in Artificial Intelligence (UAI), 2022.
Tristan Delcu, Mizu Nishikawa- Toomey, Jithendaraa Subramanian, Nikolay Malkin, Laurent Char- lin, and Yoshua Bengio. Joint Bayesian inference of graphical structure and parameters with single generative flow network. arXiv preprint arXiv:2305. 2023.
Erika Dort, Elliot Layne, Nicolas Feau, Alexander Butyaev, Bernard Henrissat, Francis Martin, Saject Haridas, Asaf Salamov, Igor Grigoriev, Mathieu Blanchette, and Hamelin Richard. Large- scale genomic analyses with machine learning uncover predictive patterns associated with fungal phytopathogenic lifestyles and traits. Scientific Reports, 13(1), 2023.
Alexei Drummond, Simon Y. W Ho, Matthew J Phillips, and Andrew Rambaut, Relaxed phyloge- netics and dating with confidence. PLOS Biology, 4(5):null, 03 2006. doi: 10. 1/.13n/ 0040088.
10

---

Published as a conference paper at ICLR 2024
Joseph Felsenstein. Maximum likelihood and minimum-steps methods for estimating evolutionary trees from data on discrete characters. Systematic Biology, 22(3):240-245 1973.
R. Fisher, L. Pusztai, and C. Swanton. Cancer heterogeneity: implications for targeted therapeutics. British Journal of Cancer, 108(3):479 485, Feb 2013.
J R Garey, J Near, M R Nonnemacher, and $ A Nadler. Molecular evidence for acanthocephala as subtaxon of rotifera. J Mol Evol, 43(3):287-299 September 1996.
Richard Hamelin, Guillaume Bilodeau, Renate Heinzelmann, Kelly Hrywkiw, Arnaud Capron, Erika Dort, Angela L. Dale, Emilie Giroux, Stacey Kus, Nick C Carleson, et al. Genomic bio- surveillance detects sexual hybrid in the sudden oak death pathogen. Communications Biology, 5(1):477, 2022.
S B Hedges, K D Moberg, and L Maxson. Tetrapod phylogeny inferred from 18S and 28S ribosomal RNA sequences and a review of the evidence for amniote relationships. Molecular Biology and Evolution, 7(6):607-633, 11 1990.
Daniel A Henk, Alex Weir, and Meredith Blackwell. Laboulbeniopsis termitarius, an ectoparasite of termites newly recognized as a member of the laboulbeniomycetes. Mycologia, 95(4):561-564, July 2003.
Diep Thi Hoang, Olga Chernomor, Arndt Von Haeseler, Bui Quang Minh, and Le Sy Vinh. Ufboot2: improving the ultrafast bootstrap approximation. Molecular biology and evolution, 35(2):518- 522, 2018.
Edward J Hu, Nikolay Malkin, Moksh Jain, Katie Everett, Alexandros Graikos, and Yoshua Bengio. GFlowNet-EM for learning compositional latent variable models. International Conference on Machine Learming (ICML), 2023.
Sebastian Höhna, Michael Landis, Tracy A. Heath, Bastien Boussau, Nicolas Lartillot, Brian R. Moore, John P. Huelsenbeck, and Fredrik Ronquist. RevBayes: Bayesian Phylogenetic Inference Using Graphical Models and Interactive Seddi-pSpcciiicttionn Language. Systematic Biology, 6564)7726-736, 05 2016. ISSN 1063-5157. doi: 10. 10.1093//yssboosyy002. URLhttps:
org arg/10.100 1101110011993 sysbio sywbss/ysooooo111111 syw021
Moksh Jain, Emmanuel Bengio, Alex Hernandez-Gareia. Jarrid Rector-Brooks, Bonaventure F.P. Dossou, Chanakya Ekbote, Jie Fu, Tianyu Zhang, Micheal Kilgour, Dinghuai Zhang, Lena Simine, Payel Das, and Yoshua Bengio. Biological sequence design with GFlowNets. Inter- national Conference Machine Learning (ICML), 2022.
Hazal Koptagel, Oskar Kviman, Harald Melin, Negar Satinianaini, and Jens Lagergren. VaiPhy: variational inference based algorithm for phylogeny. Neural Information Processing Systems (NeurIPS), 2022.
Salem Lahlou, Tristan Deleu, Pablo Lemos, Dinghuai Zhang, Alexandra Volokhova, Alex Hernández-Garefa, Léna Néhale Ezzine, Yoshua Bengio, and Nikolay Malkin. A theory of con- tinuous generative flow networks. International Conference on Machine Learning (ICML), 2023.
Clemens Lakner. Paul Van Der Mark, John P Huelsenbeck, Bret Larget, and Fredrik Ronquist. Efficiency of Markov chain Monte Carlo tree proposals in Bayesian phylogenetics. Systematic biology, 57(1))86-100 2008.
Elliot Layne, Erika N Dort, Richard Hamelin, Yue Li, and Mathieu Blanchette. Supervised learning on phylogenetically distributed data. Bioinformatics, 36(Supplement 2):1895-i902, 2020.
Jian Ma, Louxin Zhang, Bernard B Suh, Brian J Raney, Richard C Burhans, W James Kent, Mathieu Blanchette, David Haussler, and Webb Miller. Reconstructing contiguous regions of ancestral genome. Genome research, 16(12):1557-1566) 2006.
Nikolay Malkin, Moksh Jain, Emmanuel Bengio, Chen Sun, and Yoshua Bengio. Trajectory balance: Improved credit assignment in GelowNets. Neural Information Processing Systems (NeurIPS), 2022.
11

---

Published as a conference paper at ICLR 2024
Nikolay Malkin, Salem Lahlou, Tristan Deleu, Xu Ji, Edward Hu, Katie Everett, Dinghuai Zhang, and Yoshua Bengio. GFlowNets and variational inference. International Conference on Learning Representations (ICLR), 2023.
Nicholas McGranahan, Francesco Favero, Elza C. de Bruin, Nicolai Juul Birkbak, Zoltan Szallasi, and Charles Swanton. Clonal status of actionable driver events and the timing of mutational processes in cancer evolution. Science Translational Medicine, (8223):223r555-228rra5 2015.
Takahiro Mimori and Michiaki Hamada. Geophy: Differentiable phylogenetic inference via geo-
metric gradients of tree topologies. arXiv preprint arXiv:2307.03673 2023.
Bui Quang Minh, Minh Anh Thi Nguyen, and Arndt Von Haeseler. Ultrafast approximation for phylogenetic bootstrap. Molecular biology and evolution. 3065):1188-1155, 2013.
Bui Quang Minh, Heiko A Schmidt, Olga Chernomor, Dominik Schrempf, Michael D Woodhams, Arndt Von Haeseler, and Robert Lanfear. Iq-tree 2: new models and efficient methods for phylo- genetic inference in the genomic era. Molecular biology and evolution, 3735):150- 2020.
Antonio Khalil Moretti, Liyi Zhang, Christian A. Naesseth, Hadiah Venner, David Blei, and It- sik Pe'er. er. Variational combinatorial sequential Monte Carlo methods for Bayesian phylogenetic inference. Uncertainty in Artificial Intelligence (UAI), 2021.
Reza Ranjbar, Sedigheh Nazari, and Omid Farahani. Phylogenetic analysis and antimicrobial resis- tance profiles of escherichia coli strains isolated from uti-suspected patients. Iranian Journal of Public Health, 49(9):1743, 2020.
Fredrik Ronquist, Maxim Teslenko, Paul Van Der Mark, Daniel L Ayres, Aaron Darling, Sebastian Höhma, Bret Larget, Liang Liu, Marc A Suchard, and John P Huelsenbeck. MrBayes 3.2: effi- cient Bayesian phylogenetic inference and model choice across large model space. Systematic biology, 61(3):339-542, 2012.
Amy Y. Rossman, John M. McKemy, Rebecca A. sado--Shhuttheis,, and Hans-Josef Schroers. Molecular studies of the bionectriaceae using large subunit rdna sequences. Mycologia, 93(1): 100-110, 2001.
David L Swofford. Phylogenetic analysis using parsimony. 1998.
Hakon Tjelmeland and Bjorn Kare Hegstad. Mode jumping proposals in meme. Scandinaviar journal of statistics, 28( 28(1):00--223, 2001.
Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N. Gomez, Lukasz Kaiser, and Illia Polosukhin. Attention is all you need. Neural Information Processing Systems (NIPS), 2017.
Ziheng Yang and Anne D. Yoder. Comparison of Likelihood and Bayesian Methods for Estimat- ing Divergence Times Using Multiple Gene Loci and Calibration Points, with Application Radiation of Cute-Looking Mouse Lemur Species. Systematic Biology, 52(5):705-716, 10 2003.
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.
Cheng Zhang. Learnable topological features for phylogenetic inference via graph neural networks. arXiv preprint arXiv:2302.008840. 2023.
Cheng Zhang and Frederick A Matsen IV. Generalizing tree probability estimation via bayesian networks. Neural Information Processing Systems (NIPS), 2018a.
Cheng Zhang and Frederick Matsen IV. Variational bayesian phylogenetic inference. Interna- tional Conference on Learning Representations (ICLR), 2018b.
David Zhang, Corrado Rainone, Markus Peschl, and Roberto Bondesan. Robust scheduling with GFlowNets. International Conference Learning Representations (ICLR), 2023a.
12

---

Published as a conference paper at ICLR 2024
Dinghuai Zhang, Nikolay Malkin, Zhen Liu, Alexandra Volokhova, Aaron Courville, and Yoshua Bengio. Generative flow networks for discrete probabilistie modeling. International Conference on Machine Learning (ICML), 2022.
Dinghuai Zhang, Hanjun Dai, Nikolay Malkin, Aaron Courville, Yoshua Bengio, and Ling Pan. Let the flows tell: Solving graph combinatorial problems with GFlowNets. arXiv preprint 11X10 2023b.
Ning Zhang and Meredith Blackwell. Molecular phylogeny of dogwood anthraenose fungus (discula destructiva) and the diaporthales. Mycologia, 93(2):355-365, 2001.
13

---

Published as a conference paper at ICLR 2024
A PROOFS OF PROPOSITIONS AND 2
Before proving proposition [ we first prove the following two lemmas. First, we show that for two states sharing the same tree features, applying the same action to the states results in two new states still sharing the same features.
Lemma 1. Let $ = b1),,b),, (z2, b2) (2,, and s2 (11211,1) (22. ... (2gp b' be two con-terminating states sharing the same features Pi = Let a be the action that joins the trees with indices (v, w) to form new tree indexed a with edge lengths (b(ewr), b(eww)). By applying 52, we join (zy. bb) and (zinr byw) to form new tree (Zy. bqg). Then the new trees' features are equal: Pu = Pu.
on St, we join (Zv, by) and (Zw. bw) to form new tree (Zu. bu). By applying
Proof. We show that Pu can be calculated from Pv and pw:
P ((b(em))) P(b(euw))
Σ
P(a = jla, k. (p
KeE
x > jlab k, ()
k
since Pv = Ply, Pw Pw Pu Pu
p and (b(env), b(euw)) are new branch lengths for both two trees. Therefore
Next, we show that for two states sharing the same tree features, applying the same action sequences results in two phylogenetic trees with the same reward.
{(21. b1). (22,b2)... (21, bi)} and non-terminating states sharing the same features Pt p Let a be any sequence of actions to apply on $ and U to form full trees = (x,b,),x' = b(), R(z, b) (()',b').
Lemma 2. Let
{(z], (2), ... (2 by} be two
Proof. Let Pu denote the tree feature for (z, b), Plu fú II P(b(e)). We first show that the reward can be directly calculated from the root feature Pu:
[ P(azn-1) Pu
|| P(ble)
[]
P(azn-1) fú
-(b))(()lz,b)
= b(,,b),
where P(azn-1) is the constant root character assignment probability. As applied to si and s2 in sequential manner, at every step we obtain two state swith the same tree features (by Lemma until, at the final state, Pu = Pu- It follows that R(z. b) R(z', b').
We are now ready to prove the propositions.
Proposition Let 8i {(21, b1), (zz, b2) (2,, bi)) and 82 = ((zg,b)) bi)bb()(gggg (Z) b2).bb............. ... (E) b)} be two non-terminal states such that si # 52 but sharing the same features pi = p/p Let be any sequence of actions, which applied to $1 and 82, respectively, results in full weighted frees x (z, b,) = ('),b'), with two partial trajectories T (81 Pe the policy of an optimal GFlowNet with uniform PB. then PF(T) PF(T').
x). (s2
x') If
Proof of Proposition 7 Let Gay Goy be sub-graphs of the GFlow Net state graph G = (S defined by reachable states from sı and s2 in G. Since si and s2 have the same number of trees, Gs, and G. have the graph structure. Let X1, S X be the terminal states reachable from s and s2. There is thus bijective correspondence between X1 and X2: for every action set applying on ls to obtain x E X1, we obtain x' E X2 by applying on $2. Let and r' be the partial trajectories created by applying on 5i and 52, ( (T''')'') Moreover, R(x) = R(x') since 5j and 52 share the same set of features. We have the following expressions for PF(T) and PF(T'):
((x)))(T))) PF(T)= Σnj,x, R(x,) P((,)))(
R(x')n(F'(lx) PF(T') = Σxj,, R(x))P(())x))
Hence Pr(r) Pr(r')
14

---

Published as a conference paper at ICLR 2024
Proposition 2. Let $1 (21, 22.... 2)) and 52 (2 22..... ...)....) be two non-terminal states such that $1 # $2 but sharing the same Fitch features fa = f2; Vi. Let a be any sequence of actions, which, applied to $1 and 52, respectively, results in tree topologies x, E Z. with two partial trajectories (51
X' = (52 -x x'). If PF is the policy of an optimal GFlowNet with uniform Pg. then (T( Pp(T))
We use the same notation as in the proof of of Proposition [ Since and 52 have the same number of trees, Gs and Gs- have the same graph structure. Let X1.X2 be the terminal states reachable from 5j and $2. There is bijective correspondence between X1 and X2: for every action and r' be the partial trajectories created by applying a
Proof of Proposition 2
applying on 5₁ obtain x X1, obtain x' X2 by applying on $2. Let and $2. Ps(T(x) = Ps(t'lx').
For simplicity, we denote (x) (x)) It is worth to note that for two tree topologies sharing the same Fitch feature, their parsimony scores do not necessarily equal. However, when applying on two states 8 and sharing the Fitch feature, the additional parsimony scores introduced are the same:
((x) Σ M(2)= MMx) Σ M(Z))
We have the following expressions for (F(T) and Pr(T'):
-Mix
Pb(T(X)
PF(T)) e PB(T(x))
-Mix'
(F(T)= = Erjay e Pb(T,lxj)
-Mixg)
-Mixh
Σx,o,e PR(T)]x')
Eparis))
Mi Mief
We multiply
by (T( and
by PF(T') and obtain:
Mics
MEN
Mi PB(T(x)
Ei Mi S1 Mix'
PB(T(x')
PF(T)= M Mix)
PF(T')
Pb(T)lxj)
Erj,y,
PB(T)}x)
Mic
PB(TT)x)) 2g Mia
Since e
Pb(Tjlx,) for all PF(()) =PF(T).
B MARGINAL LOG-LIKELIHOOD ESTIMATION
Estimating the sampling likelihood of terminal state In the PhyloGFN state space (in both the Bayesian and parsimony-butt based settings), there exist multiple trajectories leading to the same termi- nal state x, hence the sampling probability of is calculated as: PT(x) = E - (T))). This sum is intractable for large-scale problems. However, sampling (Zhang et al.) 2022):
estimate the sum using importance
P
Σ
PF(T) PB(T)xx)
(S1)
where the trajectories t; sampled from Pn(T; x). The logarithm of the right side of (11) is. in expectation, a lower bound on the logarithm of the true sampling likelihood on the left side of (ST).
Estimating the MIL The lower bound on MLL can be evaluated using the importance-weightee bound log P(Y) 2 Ex....x-Pp log R Σ P(x,))
(x) (Burda et al. 2016). However, we cannot use it for
PhyloGFN since we cannot compute the exact pf (x), only get a lower bound on it using S1
Therefore, propose the following variational lower bound:
log P(Y) 2 ET Tu- PF log P(z) K
Σ
PB(Tlen R() bi) PF(T))
Ty/00
(z)
15

---

Published as a conference paper at ICLR 2024
We show the derivation of the estimator and thus prove its consistency:
P(Y P(Y|2,b)P(blaP((
r.b
R(z,b)P(2)
2.b
R(z.b)P(zz)
Σ PB(T|2.b)
J..b
[ R(2,b)P(2)
PF(T)
Σ
PB(T|2,b) b)
J.b
PF(T)
T: (2.b)
P PB b)
Σ
PF(T)
Jz.b (.....
PF(T)
PB(TIZr, R(r.)((..
PF(T)
PF(T)
PB(T|2r,,, by R
P(2)Br-Pr
PT(T)
Pn(TiZzi R(),( bi)
P(z)
(T(T))
T(..0...11..b)))
One can show, in manner identical to the standard importance-weighted bound, that this estimate a lower bound on log P(Y).
models edge lengths using discrete distributions. To estimate our algo- rithm's 's MLL, we augment the sampler to continuous sampler by modeling branch lengths with a piecewise-constant continuous distribution based on the fixed-bin multinomial distribution of Phy- loGFN. We can still use the above formulation to estimate the lower bound. However, each trajectory now has one extra step: 80 - domly perturbing each branch length by adding noise from U[-0.5w. 0.5w], where w is the bin size used for PhyloGFN. Let = so _ can compute Pp(T'),P(r)) from PF(T),
PhyloGFN-Bayese
(z, b) (z, by) where (z, b;) obtained from z, by ran-
(z, b) be the original trajectory in the discrete PhyloGFN,
PF(T') = (T) [FFE =a(T))=PBC)
The term i Eिe is introduced in Pr because we additionally sample over a uniform range o for all [E()] edges. The backward probability PB stays unchanged because given (z, b) generated from the discrete GFN, for any (z, b') resulting from perturbing edges, (z, b) is the the only possible ancestral state.
16

---

Published as a conference paper at ICLR 2024
C
DATASET INFORMATION
Table S1: Statistics of the benehmark datasets from DSI to DS8.
Dataset
Species
Sites
Reference
DS1
27
1949
Hedges et al.
1990)
DS2
29
2520
Garey al.
1996
DS3
36
1812
Yang & Yoder
(2003
DS4
41
1137
Henk ct al. (2003)
DS5
50
378
Lakner ct al (2008)
DS6
50
1133
Zhang & Blaakwwel]
2001
DS7
59
1824
Yoder & Yang
2004
DS8
64
1008
Rossman ct 2001
D MODELING
Given the character set Σ, we one-hot encoding represent each site in sequence. To deal with wild characters in the dataset, for parsimony analysis we consider them as one special character in Σ, while in Bayesian inference, we represent them by vector of 1l.
For both resentation feature is the set (pI...pr), where p a vector of length m| For example, for DNA sequences of 1000 sites, each p would have length 4000. Therefore, before passing these features to the Transformer encoder, we first use a linear transformation to obtain lower-dimensional embed- dings of the input features.
and PhyloGFN-Bayesian, given a state with / rooted trees, its rep-
We use the Transformer architecture /Vaswani et al.) 2017) to build the Transformer encoder. For a state with trees, the output is set 7 I features (ex, .1....,ei where e, denotes the summary feature (i.e., the [CLS] token of the Transformer encoder input).
To select pairs of trees to join, we evaluate tree-pair features for every pair of trees in the state and pass these tree-pair features as input to the tree MLP generate probability logits for all pairs of trees. The tree-pair feature for tree pair (i with representations éi Si is the concatenation of e,te; with the summary embedding of the state, i.e., the feature is [[g: (eteetee, where [+: | denotes
vector direct sum (concatenation). For state with trees, (6) I(1-1) such pairwise features are generated for all possible pairs.
To generate edge lengths for the joined tree pair (i. j), we pass [exie;tej] the concatenation of the summary feature with the tree-level features of trees and as input to the edge MLP. For unrooted tree topologies we need to distinguish two scenarios: (i) when only two rooted trees are left in the state (i.e., at the last step of PhyloGFN state transitions), we only need to generate single edge; and (ii) when there are more than two rooted trees in the state, a pair of edges is required. Therefore, two separate edge MLPs are employed, edge length is modeled by & discrete bins, the edge MLP used at the last step has an output size of & (to model a distribution over a single edge length) whereas the other edge MIP would have an output size of 2 (to model a joint distribution over two edge lengths).
For the temperature-conee conditiooee PhyloGFN, as then temperature T is passed to our PhyloGFN as an input, two major modifications are required: (i) the estimation of the partition Zo is now a function of Zo (7), which modeled by a MLP: (ii) the summary token to the Transformer encoder also captures the temperature information by replacing the usual [CLS] token with temp MLP that accepts 7' input.
E TRAINING DETAILS
Here, describe the training details for PhyloGFN. For PhyloGFN- Bayesian, our models trained with fixed 500 epochs. For ity of sampling the optimal trees, or the most parsimonious trees PhyloGFN has above 0.001. Each training epoch consists of 1000 gradient update steps using a batch of 64 trajec- tories. For e-greedy exploration, the value is linearly annealed from 0.5 to 0.001 during the first 240 epochs. All common hyperparameters for PhyloGFN-Bayesian and shown in Table S3
our models are trained until the probabil- so far, is
are
17

---

Published as a conference paper at ICLR 2024
Temperature annealing For PhyloGFN-Baaesiann Bayesian, the initial temperature is set to 16 for all exper- iments. For WavloGFN-Parsimmonn scheme, 7 is reduced by half per every 80 epochs of training. For ThyloGFN-Bayesinrt 7 is always reduced to and then fixed at 1, whereas for PhyloGFN-Parsimen Parsimony. T is only reduced when the most parsimonious trees seen by our PhyloGFN far are sampled with probability below 0.001.
7 is initialized at 4. Under the cascading temperature annealing
Hyperparameter selection For PhyloGEF-Parsimonny, Parsimony, the reward is defined R(x)
C-M(xiX) where C is hyperparameter introduced for training stability and it controls the
exp 7
magnitude of the partition function Σx R(x). Given that we cannot determine the precise value of C prior to training since we do know the value of Z as a priori, we use the following heuristic to choose C: 1000 random tree topologies are generated via stepwise-addition, and we set C to the 10% quantile of the lowest parsimony score.
Similarly, C' is employed for MhyloGFN-Bayesian under the same goal of stabilizing the training and reducing the magnitude of the partition function. Recall the reward function R(z, b) defined in
when T 1. Note that exp { can be absorbed into the partition function. For selecting the C. once again we randomly sample 1000 weighted phylogenetic trees via stepwise-addition and with random edge length, followed by setting C to the 10% quantile of the lowest - log |(Y P(b).
$4.1 it can be rewritten as R(z, b) exp C-
PhyloGFN is introduced that single trained PhyloGFN can be used sample from series of reward distribu- tions defined by different T. We modify the exxedttemperature
Pemperature-connitton
PhyloGFN-Parsimony
The temperature-condition
ThyloGFNNParsimoon by introducing
T as input in 3 places: (i) the reward function Rx;:)); (ii) the forward transition policy Py (x::; and (iii) the learned partition function estimate Ze To train the temperature--cnniiionn Phy- loGEN, the TB objective also needs be updated accordingly:
LTB(T; = log (T) I1 PF(5 (,;
PB(T = [ |Pa(s:))'
(n(T x)
Note that PB is unaffected.
During training, values for 7' are randomly selected from the range [Tmin- T. When training with state from the replay buffer, the temperature used for training is resampled and may be different than the one originally used when the state was added to the buffer.
We also employ scheme of gradually reducing the average of sampled values through the course of training: 7's are sampled from truncated normal distribution with a fixed pre-defined variance and a moving mean T that linearly reduces from Tmax to min-
Modeling branch lengths with discrete multinomial distribution When a pair of trees are joined at the root, the branch lengths of the two newly formed edges are modeled jointly by discrete multinomial distribution. The reason for using joint distribution is because under the Jukes-Cantor evolution model, the likelihood at the root depends on the sum of the two branch lengths.
We use different set of maximum edge length, bin number and bin size depending each dataset, and by testing various combinations we have selected the set with optimal performance. The maximum branch length is chosen among {0.05, 0. 0.2), and bin size w is chosen among (0.001, 0.002, 0.004). Table S2 shows the our final selected bin size and bin number for each dataset.
Table $2: Bin sizes and bin numbers used model branch lengths for DS1 to DS8.
Dataset
Bin Size
# Bins
DS1
0.001
50
DS2
0.004
50
DS3
0.004
50
DS4
0.002
100
DS5
0.002
100
DS6
0.001
100
DS7
0.001
200
DS8
0.001
100
18

---

Published as a conference paper at ICLR 2024
Table Common hyperparameters PhyloGFN- Bayesian and MayyooFN--Pasimman
Transformer encoder
hidden size number of layers
128
6
number of heads
4
learning rate (model)
5e-5
learning rate (Z)
5e-3
tree MLP
hidden size number of layers
256
3
edge MLP
hidden size number of layers
256
3
MLP (in temperaturre-ooddiionned PhyloGFN)
hidden size
128
number of layers
1
temp MLP (in lemperauure-condiioneed PhyloGFN)
hidden size
256
number of layers
3
19

---

Published as a conference paper at ICLR 2024
F CONTINUOUS BRANCH LENGTH MODELING WITH PHYLOGFN
While the original PhyloGFN-Bayesian only samples discrete branch lengths, we have experimented with new version of PhyloGFN that uses mixture of Gaussian to model and sample branch lengths effectively treating them as continuous variables. We refer to this new version as PhyloGFN- Continuous and we point out that the edge modeling is the only implementational detail that differs from the previous models.
Specifically, the edge MLP of PhyloGFN-Continn Continuuus now outputs the parameters of the Gaussian mixture which models the logarithm of branch length. These include (1) logits for selecting the components of the Gaussian mixture, (2) mean of the log branch length in each mixture, and (3) log variance of the log branch length in each mixture.
Although PhyloGFN-Coninn Continuuus continues to employ two separate edge MLP, one for the last step of state transition where only single edge needs be sampled and the other for sampling pair of edges at the intermediate step, BeoloGFN-Continoou now samples each edge independently as we have found this benefit training stability.
also worth pointing out that we no longer perform e-greedy random exploration on branch lengths, again for the sake of training stability, and the mean of the log branch length is initialized -4.0 instead of 0.0, which is more reasonable value for log branch length. There are five components in each Gaussian mixture.
Results are shown below. We are delighted to report that ChyloGFN-Continee has eliminated the quantization eπror introduced in the earlier discrete PhyloGFN Bayesian and therefore, our new model is effectively performing on par with the state of the art MrBayes and VBPI-GNN models. In particular, PhyloGFN-Contimun all datasets, with the only exception being DS7 where BoyyloFN--Continnn closely matches the performance of VBPI-GNN.
has the smallest standard deviation for the MLL estimation across
Note that ConttinnN-Cotssuuus has undergone the same training routine that is described in as is PhyloGFN-Bayesian.
Table S4: Marginal log-likelihood estimation with different methods on real datasets DSI-DS8. PhyloGFN- B(ayesian) C(ontinuous) now outperforms & A-CSMC, GeoPhy and PhyloGFN- PhyloGFN-Blayesaan) across all datascts and it is effec- tively performing par with the state of the arts MrBayes and VBPI-GNN.
MCMC
ML-based amortized, full tree space
Dataset
MrBayes
VBPI-GNN*
d-CSMC
GcoPhy
PhyloGFN-B
PhyloGFN-C
DSI
-7108.
-7108.41
-7290.36
-7111.55
-7108.95
-7108.44
DS2
-26367.57
:048
-26367.73
60.07
-30568.49
#31.34
-26368.4
10.13
-26368.90 10.38
-26367.70
+0.04
DS3
-33735.44
-0.5
-2
-0.09
-33798.06
+6.62
-33735.85
10.12
-33735.6 10.35
-33735. 11 +0.00
DS4
-13330.06
20.54
13329.94=0.08
-13582.24
=35.08
-13337.42
=1.32
-13331.83
+0.19
-13329.91
+0.00
-8214.51
=0.25
-8214.64
8367.51
=8.87
-8233.89
26.63
-8215.15 #0.2
-8214.38 #0.16
DS6
-6724.07
-6724.37
-7013.83
#16.99
-6733.91
s0.57
-6730.68
-6724.17
DS7
-37332.76
-37332.04
-37350.77
33355.99
-37331.89
+0.14
DS8
8649.88 =1.75
8660.65 =0.45
-9209.18 =18.03
-8660.48
=0.78
-665..76
20.19
-8650.46 =0.05
20

---

Published as a conference paper at ICLR 2024
G RUNNING TIME AND HARDWARE REQUIREMENTS
PhyloGFN is trained on virtual machines equipped with CPU cores and 10GB RAM for all datasets. We use one V100 GPU for datasets DS1-DS6 and one A A100 GPU for DS7-DS8, although the choice of hardware not essential for running our training algorithms.
For Bayesian inference, the models used for the MIL estimation in Table S4 are trained on a to- tal of 32 million examples, with training wall time ranging from 3 to 8 days across the eight datasets. However, our algorithm demonstrates the capacity to achieve similar performance levels with significantly reduced training data. The (able(S5] below presents the performance of PhyloGFN with 32 million training examples (PhyloGFN Full) and with only 40% of the training trajectories (PhyloGFN Short). Each type of experiment is repeated 3 times. Table Saalee/Soooce compares running time of the full experiment with the shorter experiment. The tables show that the shorter runs exhibit comparable performance our full run experiments, and all conclude within 3 days.
We compare the running time of PhyloGFN with VI baselines (VBPI-GNN, Vaiphy, and GeoPhy) using the DS1 dataset. VBPI-GNN and GeoPhy are trained using the same virtual machine config- uration as PhyloGFN (10 cores, 10GB ram, IxV100 GPU). The training setup for both algorithms mirrors the one that yielded the best performance as documented in their respective papers. As for VaiPhy, we employed the recorded running time from the paper on machine with 16 CPU cores and 96GB RAM. For PhyloGFN, we calculate the running time of the full training process (PhyloGFN- Full) and four shorter experiments with 40%. 24%, 16% and 8% training examples. The table SS7 documents both the running time and MIL estimation for each experiment. While our most com- prehensive experiment, PhyloGFN Full, takes the longest time to train, our shorter runs all of which conclude training within a day show only a marginal degradation in performance. Remark- ably, even our shortest run, PhyloGFN 8%, outperforms both GeoPhy and d d-CSMC, achieving this superior performance with approximately half the training time of GeoPhy.
Table S5: PhyloGFN Bayesian MLL estimation on 8 datasets. We repeat both full experiment and short exper- iment (with 40% training examples) 3 times
Experiment
PhyloGFN Full
PhyloGFN Short (40% Training data)
Repeat
2
2
DSI
-7108.95 20.36
-7108.97
10.05
-7108.94
20.05
-7108.97
a0.14
-7108.94
.0.22
-7109.04
10.08
DS2
-26368.9
40.28
-26368.77
-0.43
-26368.89
+0.29
-26368.9
+0.39
-26369.03
10.31
-26368.88
60.32.
DS3
-33735.6 +0.35
-33735.60
-0.40
-33735.68
+0.64
-33735.9
+0.91
-33735.83
+0.62
-33735.76
=0.75
-13331.83
+0.19
-13331.80
20.31
-13331.94
20.42
-13332.04
=0.57
-13331.87
#0.31
-13331.78
=0.37
DSS
-8215.15
-8214.92
20.27
8214.85
20.28
-8215.38
20.27
-8215.37
20.26
-8215.38
DS6
-6730.68
=0.54
-6730.72
+0.26
-6730.89
=0.22
-6731.35
+0.31
-6731.2
+0.4
-6731.1
=0.38
DS7
-37359.96
+1.14
-37360.59
=1.62
-37361.51
+2.89
-37362.03
=5.2
-37363.43
#2.2
-37362.37
=2.65
-8654.76
20.19
-8654.67
20.39
-8654.86
20.15
-865588 20.95
-8655.65
=0.37
-8654.96
20.46
Table S6: PhyloGFN-Bayesian training time
Dataset
PhyloGFN Full
PhyloGFN Short
DS1
62h40min
20h40min
DS2
69h16min
28h
DS3
80h20min
35h40min
DS4
103h54min
44h30min
DS5
127h50min
51h40min
DS6
135h10min
53h10min
DS7
174h3min
60h20min
DS8
190h25min
61400in
H ABLATION STUDY
H.1 BRANCH LENGTHS MODEL HYPERPARAMETERS
PhyloGFN models branch lengths using discrete multinomial distributions. When estimating MLL, the branch length model is transformed into piecewise-constane continuous form, introducing small quantization error. Two hyperparameters define the multinomial distribution: bin size and bin number. This analysis investigates the impact of bin size and bin number MIL estimation.
21

---

Published as a conference paper at ICLR 2024
Table S7: Running time of WawylooFN--Byyesiin and VI baseline methods DS1
Running Time
MLL
VBPI-GNN
16h10min
-7108.41
#0.14
GeoPhy
12h50min
-7111.55
+0.07
*-CSMC**
-2h
-7290.36
+7.23
PhyloGFN Full
62h40min
-7108.97
#0.05
PhyloGFN 40%
20h40min
-7108.97
#0.14
PhyloGFN 24%
15h40min
-7109.01
+0.15
PhyloGFN 16%
10h50min
-7109.15 +0,23
PhyloGFN 8% 5h10min
-7110.65
+0.39
For the fixed bin size of 0.001, assess four sets of bin numbers: 50, 100, 150, and 200, and for the fixed bin number 50, we evaluate three sets bin sizes: 0.001. 0.002, and 0.004. For each setup, train PhyloGFF-Bayesian model on the dataset DSI using 12.8 millions training examples.
Table S8 displays the MLL obtained in each setup. A noticeable decline in MIL. estimation occurs as the bin size increases, which is expected due to the increased quantization error. However, MLL estimation also significantly declines the number of bins increases over 100 under the fixed bin size of 0.001. We conjecture this is because the size of the action pool for sampling the pair of branch lengths increases quadratically by the number of bins. For example at 200 bins, the MLP head that generates branch lengths has 40,000 possible outcomes, leading to increased optimization challenges.
H.2 EXPLORATION POLICIES
PhyloGFN employs the following techniques to encourage exploration during training:
e-greedy annealing: set of trajectories are generated from the GFlowNet that is being trained, with actions sampled from the GFlowNet policy with probability 1 - and uni- formly at random with probability The e rate drops from pre-defined Extarr to near through the course of training.
Replay Buffer: a replay buffer used store the best trees seen date, and use them for training, random trajectories are constructed from these high-reward trees using the backward probability Pn-
the training of GFlowNet begins large temperature T and is divided half periodically during training.
Temperature annealing:
To assess the effectiveness of various exploration methods, we train PhyloGFN-Bayesiee on DS1 with the following setups:
1. All trajectories are generated strictly on-policy training without any exploration methods (On-policy).
2. All trajectories are generated with epsilon-greedy annealing (E).
3. Half of trajectories are generated with e-greedy annealing, and the other half are con- structed from the replay buffer (e RP).
4. The same setup as 3, with the addition of temperature annealing. T initialized at 16 and reduced by half per every 1.5 million training examples until reaching I (E + RP T Cascading).
5. The same setup 4, except the temperature drops linearly from 16 (E + RP+ T Linear).
For each setup, we train model using 12.8 millions training examples. Table $9 displays the MLL estimation for each setup. On-policy training without any additional exploration strategies results in the poorest model performance. Epsilon- -reedggreyy annealing significantly enhances performance, and combining all three strategies yields the optimal results. There no significant difference between cascading and linear temperature drops.
I SIGNIFICANCE OF MODELING SUBOPTIMAL TREES WELL
Several downstream tasks that rely on phylogenetic inference require not only the identification of optimal (maximum likelihood most parsimonious) trees, but also the proper sampling of suitably
22

---

Published as a conference paper at ICLR 2024
Table S8: MLL estimation of PhyloGFN-Bayeste DSI with different combinations of bin size and bin number model edge lengths.
Bin Size
Bin Number
MLL
0.001
50
-7108.98 #0.14
0.001
100
-7108.96 +0.18
0.001
150
-7109.17 +1.14
0.001
200
-7109.3 +7.23
0.002
50
-7109.95 +0.29
0.004
50
-7114.48
#0.52
Table 89: MLJ. estimation of PhyloGFN-Bayesinn on DS1 with different exploration policies
Training Policies
MLL
On policy
-7307.99 =0.26
-7109.92 a0.18
é + RP
-7109.95 20.22
RP + T Cascading
7108.98 20.14
E + RP + T Linear
7108.96 20.15
weighted suboptimal trees. In evolutionary studies, this includes the computation of branch length support (i.c., the probability that given subtree present in the true tree), as well as the estima- tion of confidence intervals for the timing of specific evolutionary events (Drummond et al. p006 These are particularly important the very common situations where the data only weakly define the posterior tree topologies (e.g. small amount of data per species, high degree of divergence between sequences). In such cases, because suboptimal trees vastly outnumber optimal trees, they contribute to non-negligible extent the probability mass of specifie marginal probabilities. Full modeling of posterior distributions on trees is also critical in studying tumor or virus evolution within patient (McGranahan al. 20155 Buendia et al.) 2209), e.g., to properly assess the prob- ability of the ordered sequence of driver or passenger mutations that may have led to metastasis or drug resistance (Fisher et al. 2013))
23

---

Published as a conference paper at ICLR 2024
J REPRODUCIBILITY
In this section, we describe the procedure to reproduce our results. Our most recent implementation involving continuous branch length modeling can be accessed through the following link: https: //githab. com/ zmy1 6/phy 1 ogfn/ For parsimony analysis and bayesian inference with discrete modeling, refer to code and data provided in the supplementary material.
J.1 MRBAYES
To reproduce the exact MLL results from (Zhang Matsen IV 2018a), MrBayes needs be in- stalled with Anaconda.
conda install bioconda: :mrbayes
With MrBayes built from source, the resulted M.L.score differs slightly, as pointed out (Mimori Hamada 2023). We use the following script compute MLL for DS1-DS8.
execute INPUT_FOLDER/ds1. nex 1set nst=1 prset statefrregpr=fixeeeeeue prset bp
ngen=10000000 nruns=10 nchains=4 printfreq=1000 \ samplefreq=100 savebrlens-yer fiLenam=QOTPOT_FPDDDE
J.2 VBPI-GNN
We use the script provided in the paper's GitHub page to generate models for DSI-DS8:
python main.py py -ddataset DS1 --brlen_mddel gmn --gnn_type edge
\
-- 2 --hdim 100 --maxTter 400000 --empF req --psp
VBPI-GNN requires mre-genraation of trees. In the original work, the trees topologies are generated from the algorithm UFBoot (Minh al. [0013] Hoang et al.) 2018). We use the program iqtree (Minh et al. 2020) to run UFBoot on all datasets with following command:
iqtree -s DS1 -b 10000 wbb
JC69
-redo
To generate randomly perturbed trees to compare log sampling density versus log posterior density. Run the following notebook in our supplementary data folder:
ipynb
J.3 GEOPHY
We the following script to GeoPhy data set DS1. The script is created based on the best configuration recorded in Mimori & Hamada (2023)
-ip ds ds-data/ds1/D$1.o. -op resultss/dsl/dd__ppeee -c sconfig/default.yyan -s -ua
python //scripts/zun_p...
-es lorentz -edt full -ed -eqs le-1 -ast 100_1000 -mes
-ci 1000 -ms 1000_000 -ul
J.4 PHYLLOGEN-BAYESIAN
Training configuration files for DS1-DS8 in the following folder in the supplementary data:
bayyssian
Run the following script to train the PhyloGFN-Bayesian model:
cd phylo_gfn
python run ru catfig///ennnnrrr....ateee/r/rrs yyaml
$output_path --nb_device 2 \
-racc --acc_gen_sampline ing - 1000 et
24

---

Published as a conference paper at ICLR 2024
To compute log MLL, plot sampled log density versus log posterior density and compute the Pearson correlation, run the following notebook in supplementary data:
phylo_gfn/phyt http://yoogotoaayiaiiaaaeeeeeeee
ipynb
Models for repeated experiments of full training and short training are trained with memory effi- cient PhloGFN implementation in folder phylo gfn bayesian memoryefficient
The full training scripts are folder
The short training scripts are folder
Training scripts for ablation study are in folder
J.5
PIYLLOGEN-PARSIMOYY
Training configuration files for DS1-DS8 are the following folder the supplementary data:
Run the following script to train the PhyloGHN-Parsimony model:
ed phylo_gfn
python run, -y http:///wg///aarcccc...a.aatttt/ cklkk Soutput_path --nb_devis 2
yaml \
--acc_gee__same ing --random_sanpe ing_batch_siize 1000 --quiet.
generate the plots Fig. S2
run the following script:
ed phylo_gfn python parsimony_eval.pe
@oodel_foldez
csequences_pachh
J.6 TEMPERATURE-CNNDITONEDD PHYLOGFN
To reproduce the result in Fig. 3 first run the following script train the model:
ed phylo_gfn
python run.py py dataset/benchmark, httpee//eeeeaaaaaccaatttttteee/e waw.ddaattatsssssce pickle $output_pan path --nb device 2 \
http::///onnnntt....aaaattttrrrs cotp://on//ddsccoooiiiiolaooaal. yaml
--nb_device
--ace_gen sesee_samlinng -- sampoon_ammpling bize 1000 --quiet
To generate the plots, run the following notebook in the supplementary data:
dsl httph//ioiccccniiionnnn.....ll.. ipynb
25

---

Published as a conference paper at ICLR 2024
K ASSESSING SAMPLING DENSITY AGAINST UNNORMALIZED POSTERIOR
DENSITY
Here, we assess PhyloGFN-Bayssinn 's ability to model the entire phylogenetic tree space by com- paring sampling density against unnormalized posterior density over gegh/meddim/loww-pobaaee regions of the tree space. To generate trees from medium and low probability regions, we use trained VBPI-GNN model, the state-of-the-art variational inference algorithm.
We first provide brief introduction to how VBPI-GNN generates a phylogenetic tree in two phases: (i) it first uses SBN to sample tree topology: (ii) followed by GNN-based model to sample edge lengths over the tree topology. SBN constructs tree topologies in sequential manner. Starting from the root node and set of sequences Y to be labeled at the leaves, the algorithm iteratively generates the child nodes by splitting and distributing the sequence set among the child nodes. The action at each step is thus to choose how to split a set of sequences into two subsets.
To introduce random perturbations during tree construction, at every step, with probability e, the action is uniformly sampled from the support choices. Given a well-trained VBPI-GNN model, we sample from high/mediuum/low www.dddiu//ww-pbbbbblllyyyyy regions with 0%, 30% and 50% probability. We apply PhyloGFN-Bayen
and VBPI-GNN on these sets to compute sampling density and compare with the ground truth unnormalized posterior density computed as P(Y zz, Pz, b). Note that VBPI- GNN models the log-edge length instead of the edge length. Hence we perform change of variables when computing both the sampling probability and the unnormalized prior.
Fig. shows scatter plots of sampling density versus ground-truth unnormalized posterior density of datasets DS2 to DS8, complementing the result on DSI shown in Fig. [ We can see that while VBPI-GNN performs well on high-probale probbabbbbyyy regions, our method is better at modeling the target distribution overall. Fig. shows that our model performs poorly at modeling the high-probability region for DS7. The result aligns with our MLL estimation, shown in Table S4 Further work is required to investigate the cause of PhyloGFN's s poor modeling on DS7.
trees
50% trees
Unverno bost Log Pooter Denvety
-26155
20150
26700
26259
-26255
25150
-26170
-26175
28453
-21550
-2400
-24400
-26199
-26152
PasistrRud 948
Pearonsttdd 914
120
210
-102
50
100
150
200
-100
-30
100
150
200
250
VBP GNN on 076 random trees
VIIR-GNN on 30% random trees
VIPI-GNN on 30% random trees
wwworrrreeddoo Log Rootstor Diee2
-25153
28203
20253
-25383
25250
-26175
26193
24550
26380
22400
-2440
25193
25450
190
200
250
db 160 sio
200
zio 160
350
50 103 150 záa zia 160
350
Sampling Donity
Samsling Danily
Samplin Log Duralty
(a) DS2
Figure SI: Sampling log density vs. ground truth unnormalized posterior log density for DS2-DS8
26

---

Published as a conference paper at ICLR 2024
31500
33150
39415
31250
-33160
33000
-33000
13650
-33160
-33700
23730
-33305
-31000
VSHI-CNN on 0% random trees
VBII-GNN on 50% random trees
VBPI-ENI on 50% random trees
-55410
13500
-31350
-30175
33100
-3300
-33485
-31150
-33150
31650
33100
-32700
-31
-3350
-33205
Pearsee..n
Powmonkonn Sc3
11800
230
263
25c
260
20
50
160
150
200
X6
30d
360
tin
50
100
340
200
350
350
Samplinn Devdit
Samsling log Dondty
Sampl no Seniity
(b) DS3
12000
Unnomeeteed Rastene Centry
15230
11000
11550
-11090
17000
-15100
-13330
-13120
-13159
-13240
15340
AvaroeR-n 602
Poymonki 996
-13200
Anarswn
240
aio
180
136
225
23c
sie
VBPI-GNN un 0% random trees
VBIN-GNN en 30% randim Cruss
VSPI-GNII on 50% ran Dres
-10000
-15230
-3377
Uivemaa bac Puttets
13000
113220
-13000
-13125
13240
11340
-11175
13250
Pearown 666
-11180
Premtonk 6.77)
283
200
300
110
30n
353
330
225
750 275 300 305 350
375
200
750
100
305
Samoliing Dendity
Samsing log Dandty
Sampl Sanctyy
(c) DS4
GPN
GFN 50%
-7900
3600
-79.0
5000
-720
-7300
-3550
-4100
252
203
152
202
VBPI GNN 0% random treas
VBPL-NNN 30% random treas
VBPI-GNN on 50% random trees
2425
-7550
-7500
7715
-7920
-000
7930
-0225
-200
-8230
-$300
-0275
Pransrfi-6.155
-8100
-225
Perracnt-0.007
$150
Pr 3.72
320
Samplng
Sumplnng Dersity
Semp Log Daraty
(d) DS5
Figure S1: (cont.)
27

---

Published as a conference paper at ICLR 2024
50%
-0770
4100
0350
-6375
-1400
-4100
-140
-5150
5450
-6320
6310
-0300
-8340
330
420
202
252
420
200
250
VBPI-GNN 8% random trees
VRP&-NN 30% random trees
VBPI-GNN
50% random trees
-4350
-010
-4400
-0000
-010
-630
-5200
-1500
6030
Prensrfi-..157
te 2.73
250
460
sio
Sempling Dora
Sumplnng Deraiy:
Semp Log Duntty
(e) DS6
GFN en 098 dom troes
GFN on 30% random trees
GFN on 50% random trees
59100
36900
30220
17020
29100
37100
37150
-37200
337300
faamanke2.787
-17420
300
325
250
375
405
425
450
300
200
303
400
-120
0
100
200
200
400
VSPI GNN random
30%
GNN
50%
35330
30050
-35340
-34550
-27130
-27100
-37200
25370
-37550
-37320
-35000
420 Sump in bereay
Comp in Beeeity
Sampll Lse Deniity
(1) DS7
on 0% random trees
random Diees
GFN on 50% random
-1200
5150
8120
-0200
-3250
-1300
-8130
8170
PramorR: 6.578
Peerennk 6413
Premrnnk 6000
483
500
510
560
212
303
350
400
<to
Sto
200
253
300
350
VBR--NNN
VBPI-GNN
-0120
8150
-1200
1220
-3200
-0120
01:0
-1250
-5400
-2550
$450
zla
560
sio
400
gg) DS8
28

---

Published as a conference paper at ICLR 2024
L PARSIMONY ANALYSIS RESULTS
-6.0
GFlowNet log probability
6.5
-7.0
-7.5
-8.0
-8.5
4026
4027
4028
4028
4029
4029
4029
4029
4030
4030
4030
4030
4030
parsimony
single most-parsimonious tree with score of 4026 been identified for DS1.
6.5
-7.0
GF owNet log probabi ity
-8.0
-8.5
-9.0
-9.5
6223
6224
6225
6225
6225
6226
6226
6226
6226 6227
6227
6227
6227
parsimony scores
(b) A single most-parsimonious tree with score of 6223 has been identified for DS2.
6.5
-7.0
GFlowNet log probability
-7.5
-9.0
-9.5
-10.0
6659
6659
6661
6661
6662 6662 6662
6663
6663
6663
6663
6663
parsimony scores
(c) Two most-parsimonious trees with score of 6659 have been identified for DS3.
Figure S2: For each dataset, the sampling probabilities for the most parsimonious well several suboptimal estimated from learned PhyloGFNs. Each hoxplot represents single unique unrooted tree and shows a distribution log-probabilite over its 3 rooted versions where n denotes the number species. The dashed bounding box groups the set of purlll--prsimooiouss trees.
29

---

Published as a conference paper at ICLR 2024
GFloweet log probablity
-10
-12
14
2424
2424
2424
2424
2425
2425
2425
2425
2425
2426
2426
2426
2426
2426
2427
2427
2427
2427
2427
2428
2428
2428
2428
2428
parsimony scores
(d) Four most-parsimonions trees with score of 2424 have been identified for DS4,
GFlowNet log probabiltty
-10
-12
14-
-18
1491
1491
1492
1492
1492
1492
1492
1493
1493
1493
1493
1493
1494
1494
1494
1494
1494
parsimony scores
(c)
Two most-parsimonious trees with a score 1491 have been identified for DS5.
GFlowNet log probabiltty
-10
-12
-14
16:
-18
879
879
879
879
879
880
880
880
880
880
881
881
881
881
B81
882
882
882
862
882
parsimony scores
(0 12 most-parsimonious trees with score of 879 have been identified for DS6.
Figure S2: (cont.)
30

---

Published as a conference paper at ICLR 2024
11
11
m
GF owNet log probability
12
14
-16
18
-2
-22
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7150
7151
7151
7151
7151
7152
152
152
152
7152
7153
7153
7153
7153
7153
parsimony scores
g
most-parsimonious trees with score 7150 have been identified DS7.
GFlowNet log probabllty
-10
-12
-14
-6
-88
1461
1461
1461
1461
1461
1461
1461
1467
1461
1467
1467
467
467
1461
1463
1463
1463
1463
1463
1464
1464
1464
1464
1464
parsimony scores
(h) 21 most- ponstmaaniooious frees with score of 1461 have been identified for DS8,
Figure $2: (cont.)
31

---

Published as a conference paper at ICLR 2024
M TRAINING CURVES
The plot in Figure $3 illustrates the training curves for PhyloGFN-Bayesian across all eight datasets. In each dataset plot, the left side displays the training loss for every batch of 64 examples, while the right side shows the ML.. computed for every 1.28 million training examples.
Loss
MLL
101
-7150
-2000
10a
-7250
10
-1300
-7350
107
-7400
10a
-7450
2.0
3.0
0.5
3.0
Training examples
Training examples
(a) DS1
Loss
MLL
30
-26400
Loos
26600
-26700
-26800
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Training examples
1e7
examples
1e7
(b) DS2
Loss
MLL
10a
-33800
34200
Lo0
34400
10a
34600
10a
10a
-35000
Training examples
1e7
Training examples
lel
(c) DS3
Figure S3: Training curves for DS1-DS8
32

---

Published as a conference paper at ICLR 2024
-13400
13800
107
-14000
107
14200
14400
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
255
3.0
examples
le7
examples
le7
DS4
Loss
MLL
-8200
107
-8400
LOG3
10a
-8800
10a
9000
-9200
2.0
2.5
3.0
0.0
0.5
1.0
2.5
Training examples
1e7
Training examples
1e7
(c) DS5
Loss
MLL
-6800
-7200
1033
30a
-7600
10a
-7800
10a
1.0
1.5
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
3.0
Training examples
Training examples
DS6
33

---

Published as a conference paper at ICLR 2024
Loss
MLL
10a
39000
-40000
101
Lo05
-41000
10a
-43000
-44000
107
0.5
Training examples
1e7
Training examples
1e7
(g)
Loss
107
10c
8800
-9000
Loo5
304
10a
-9400
-9600
101
10F
0.0
0.5
1.0
1.5
2.0
2.5
3.0
0.0
0.5
1.0
1.5
2.0
2.5
3.0
Training examples
1e7
Training examples
1e7
(h) DS8
34


Reviews:
[{"review_id": "t3wjERFqei", "reviewer": "Reviewer_tiKV", "summary": "The paper presents an application of GFlowNets for phylogenetic inference. The approach encompasses both Bayesian posterior inference and parsimony-based inference (MLE with the assumption of site independence). To address the continuous nature of branch lengths, they are transformed into discrete bins, making the phylogenetic inference problem discrete. GFlowNets are parameterized by neural networks, teaching them a policy to generate phylogenetic trees by selecting and connecting pairs of subtrees. With its probabilistic policy, the system can produce a distribution of phylogenetic trees suitable for posterior inference. In the MLE scenario, the model employs a temperature parameter that, when annealed, ensures the GFlowNet samples align with the MLE samples. The training of GFlowNets utilizes standard trajectory balance. The results suggest that the approach matches the performance of variational inference-based methods for Bayesian posterior inference, and in the MLE case, it is comparable to traditional greedy heuristic-based search algorithms.", "questions": "- *Time Efficiency*: Given the comparable performance of GFlowNets and PAUP*, can the authors provide information on the computational efficiency (wall-clock time) of GFlowNets to discern its advantages or disadvantages?\n- *Utility of Suboptimal Structures*: Can the authors elucidate the practical implications of having estimates for suboptimal structures in the context of phylogenetics? Is there a tangible benefit in real-world applications?", "limitations": "", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- *Innovative Application*: Utilizing GFlowNets for phylogenetic inference is a novel idea, showcasing the versatility of GFlowNets in unique problem settings.\n- *Competitive Performance*: The approach not only matches up to VBPI-GNN-based methods in Bayesian posterior inference but offers capabilities beyond them, like generating from arbitrarily filled-in subtrees which preceding methods could not tackle.\n- *Estimating probability of suboptimal structures*: The ability to outperform VBPI methods in estimating probabilities of suboptimal structures is a notable accomplishment.", "weaknesses": "- *Performance vs. Efficiency*: While GFlowNets might perform comparably to PAUP* in the parsimony-based inference setting, the real differentiator would be computational efficiency on a new inference task. Unfortunately, no wall-clock time data is provided, making it challenging to discern any advantages of GFlowNets in this scenario.\n- *Methodological Novelty*: The paper does not seem to bring forth significant machine learning methodological advancements, with much of the methodology being straightforward applications without any notable novel ML innovations for the particular problem at hand. The biggest methodological contribution is the use of discretized bins for branch lengths, which is a compromised solution.\n- *Utility of Results*: While the GFlowNets surpass VBPI methods in estimating probabilities for suboptimal structures, the practical significance of this in the context of phylogenetic applications remains unexplained. The experiments focus solely on the accuracy of probability estimates and fail to provide insights into their relevance for the broader application.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "8BzFnC0cwW", "reviewer": "Reviewer_iXKc", "summary": "This paper tackles the related problems of Bayesian phylogenetic inference and maximum parsimony phylogenetic inference.  These problems have recently seen a lot of development in the machine learning space with a flurry of variational approaches.  The novelty in this paper is framing the problem as a Markov Decision Problem, where one needs to sequentially build a tree from the bottom up by joining the roots of rooted trees (starting with each leaf being its own rooted tree with one node).  This paper tries to learn a good policy for the MDP (i.e., a generative model for building trees) using the recently developed GFlowNets.  The paper puts some effort into finding provably good features to use when learning the policy, and showing that the optimal policy would, in fact, result in a distribution over trees equivalent to the posterior.  The authors apply their approach to a number of standard benchmarking datasets and find reasonable performance, particularly for low probability tree topologies.", "questions": "* I am being a bit of a devil's advocate here, and it is a minor point, but, who cares about extremely unlikely trees?  Difference in log-likelihood on Figure 2 suggests that the trees with lowest posterior support are about $10^{-66}$ less likely than the trees with highest posterior support.  Does saying that those trees have zero probability really matter? Is there a real world use-case where knowing the posterior probability of those trees more accurately would be useful?\n* Would it be possible to use the machinery presented in this paper just for topologies in the Bayesian setting by somehow marginalizing out branch lengths at each step?  This is obviously not necessary for the present manuscript, I am just curious.\n* A preprint that was posted shortly after the ICLR deadline is highly relevant: https://arxiv.org/abs/2310.09553.  Since that was posted after the ICLR deadline, it did not influence my review of this paper.  The method in that preprint, ARTree, is based on reinforcement learning and also seeks to learn a posterior via framing a tree as the outcome of sequential tree building.  Yet, the training objectives and MDP formulation are substantially different, and so I do not see ARTree as reducing in any way the novelty of the present work.  I just bring this up in the hopes that it is interesting/helpful to the authors.", "limitations": "", "rating": 8, "confidence": 4, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "* The paper is clear and well-written.\n* The approach is interesting, conceptually simple, and provides a nice, efficient way to generate distributions over trees that support all of tree space (as opposed to relying on a pre-defined subset of tree space like VBPI).\n* The theoretical results and connections to Felsenstein's and Fitch's algorithms are really nice.\n* The performance across the full space of trees is promising.", "weaknesses": "* I know that it is common in the Bayesian Phylogenetics field, but I am uncomfortable with using the Marginal log-likelihood (MLL) as a measure of the accuracy of the posterior.  As the authors note, taking $K \\to \\infty$ in equation (5) or the log of (S1)results in the true MLL regardless of the distribution used.  For finite $K$, both the bias and variance of the estimated MLL will depend on the learned posterior and it is incredibly difficult for me to compare methods.  For example in Table 1, on DS7, GeoPhy is bolded for having the smallest mean MLL, but its standard error is huge, which suggests to me that by some measures it may not have learned a very good posterior.  Is there a different task that could get at the quality of the posterior in a more interpretable way?\n* I also find it somewhat surprising how much VBPI-GNN outperforms PhyloGFN on the MLL estimation task.\n* Saying that VBPI-GNN is \"severely\" limited in its applicability to postulating alternative phylogenetic theories feels like far too strong of a statement, especially in light of Table 2 -- for trees with non-negligible posterior probability VBPI-GNN is quite accurate.  I agree that VBPI being unable to put mass on all of tree space is conceptually displeasing, but I don't think the presented evidence supports the claim that VBPI would result in any real-world failures of inference.\n* Are the axes flipped on Figure 2? If I understand correctly the x-axis should be the unnormalized posterior under the long run of MrBayes, and the y-axis is the unnormalized posterior for either GFN or VBPI.  Why do the points fall on different x-axis ranges for the two columns but very similar y-axis ranges?  (Similarly for Figure S1)\n* It would be nice to include some information about runtime and to include some training curves.  Are these models difficult to train?  How sensitive is training to the choice of distribution over $\\tau$?  Etc...\n\nTypos:\n* \"Given a set of observed sequence\" --> \"Given a set of observed sequences\"\n* \"Each action chooses a pair of trees and join them\" --> \"Each action chooses a pair of trees and joins them\"\n* Is it important or a typo in Proposition 1, Lemma 1, and Lemma 2 that only the branch lengths of the first tree are different?  I.e., it is $(z_1, b_1)$ and $(z_1', b_1')$ but then $(z_2, b_2)$ and $(z_2', b_2)$ and so on.  I don't see that used anywhere.\n* In the statement of Lemma 1 I believe it should be $(b(e_{uv}), b(e_{uw}))$ not $(b(e_{uv}), (b(e_{uv}))$.\n* This sentence needs substantial rewording: \"Note that $R(x) \\ne R(x') even $s_1$ and $s_2$ share the same Fitch feature because two trees can have different parsimony score when their root level Fitch feature equals.\"\n* In the proof of proposition 2, I believe that one must multiply by $\\frac{\\exp \\frac{\\sum_i M(z_i)}{T}}{\\exp \\frac{\\sum_i M(z_i)}{T}}$, not $\\frac{\\frac{\\sum_i M(z_i)}{T}}{\\frac{\\sum_i M(z_i)}{T}}$.\n* At the bottom of the first paragraph on p. 17 should it be $-\\log P(\\mathbf{Y} | z,b) P(b)$ not $-\\log P(\\mathbf{Y} | z,b) P(z)$ in order to match what is in the reward function?\n* \"ground-trueth\" --> \"ground-truth\"", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "VoVjvYS6NR", "reviewer": "Reviewer_SMRb", "summary": "The paper develops a method for inferring phylogenetic trees, i.e. graphical representations of the evolutionary relationships of species, using Generative Flow Networks (GFlowNets). The GFlowNet treats the tree building processes as a reinforcement learning problem, where the action set corresponds to joining the roots of the subtrees existing in the pool. The paper uses a transformer architecture to encode the input states that correspond to a set of features extracted from subtree structures.", "questions": "Are the \\pm values given in Table 1 standard deviation or standard error? What are these standard deviations/error over? Experiment repetitions? If standard deviation across repetitions, the results may be good enough. They may be alarming otherwise.\n\nOverall this is an interesting piece of work with decent potential for impact. I also give sincere value to the effort the authors make to use advanced machine learning methods for such hard applications of natural sciences. However, the paper requires some work to improve, especially in terms of presentation and clarification before being ready for publication. The case looks to me like borderline at present but has potential to improve towards an accept after the rebuttal.\n\n---\n\nThe rebuttal has addressed my concerns, so I raise my score to an accept.", "limitations": "", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* Both the proposed model architecture and the use of GFlowNets for the chosen application are novel and interesting.\n\n * The paper conducts a large body of well-planned experiments, compares against a properly chosen set of baselines and reports results favorable to the proposed method. \n\n * The fact that the comparison is not made against many alternative methods is understandable, as there probably are not many modern machine learning methods addressing the same problem.\n\n* The design choices used in the model architecture are well justified, for instance the one given at the end of Page 6 for the transformer makes perfect sense.", "weaknesses": "* The biggest weakness looks to me like the results in Table 1. The log-likelihood scores look very similar to each other. For instance -7108.95 for PhyloGFN and 7108.95 VBPI-GNN. Similar for other data sets.\n\n * It is great that the paper makes lots of effort to ensure the reproducibility of the results. However, It looks to me like the main paper lacks some essential details about the experiments, making it a bit hard for the reader to evaluate the results. See my question below. \n\n * Likewise, the paper would be more readable by the machine learning community if the used tree-level features are explained. As far as I was able to detect, the main paper mentions only that they are Fitch and Felsenstein features, which may be obvious to an evolutionary biologist but they do not tell anything to me as a machine learning researcher. Now that this is not a biology journal but a machine learning research venue, an introduction to such basic concepts somewhere in the paper could be beneficial.\n\n * The paper reports log-likelihood results, which measures model fit. It measures parsimony, which appears to be about the computational aspect. The log densities reported in Figure 2 is an indicator of diversity. It may be better to have a more direct and interpretable score of the discovery performance, for instance prediction accuracy of evolution tree links or ancestral relationship detection between pairs of species in cases where there is agreement on the ground truth. The current result landscape looks a bit too exploratory. While the results make intuitive sense, they still leave many gray areas in their detailed interpretation.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "FCi2nnFmvy", "reviewer": "Reviewer_7zzf", "summary": "This is a nice application of GFlowNets to the phylogeny problem in computational biology. There is a main issue that needs to be resolved, namely the performance. First the marginal likelihood comparison with other methods, are the PhyloGFN better in terms of lower bound and what is the running time required to get these results. Second, regarding the better estimates of posterior probabilities for low posterior  phylogenetic tree (which may be the main selling point for this paper), what is the gold standard and how does your running time compare with that of he gold standard and those of the other methods? Running time doesn't have to be measured as wall-clock or using a formal analysis, but could be based on the number of so-called Felsenstein peeling operations (or whatever that make it likely that this approach may with sufficient effort yield a efficient method in the future, if that isn't already the case). The proper resolution of these issues should make this become a good contribution.", "questions": "high complexity of tree space\n\n■ Page 1\nIt is large! complex is less clear.  \n \n\n \n\nPhyloGFN is competitive with prior works in marginal likelihood estimation and achieves a closer fit to the target distribution than state-of-the-art variational inference methods\n\n■ Page 1\n\n \n\nI don’t know, but the better fit for lower part of the posterior is worth mentioning. It may be key.\n\n \n\n \n\ncontinuous vari- ables that capture the level of sequence divergence along each branch of the tree.\n\n■ Page 1\n\n \n\nAlso the ml version has those.\n\n \n\n \n\nComing from the intersection of variational inference and reinforcement learning is\n\n■ Page 1\n\n \n\nReformulate.\n\n \n\n \n\nPhyloGFN leverages a novel tree representation inspired by Fitch and Felsenstein’s algorithms to represent rooted trees without introducing additional learnable parameters\n\n■ Page 2\n\n \n\nMake this clearer. In particular, in addition to what?\n\n \n\n \n\nexplores\n\n■ Page 2\n\n \n\nNo it has capacity or potential to do this.\n\n \n\n \n\nstate-of-the-art MCMC\n\n■ Page 2\n\n \n\nThe comparison must be made relative to resources since you use mcmc as gold standard\n\n \n\n \n\nWith their theoretical foundations laid out in Bengio et al. (2023);\n\n■ Page 2\n\n \n\nIs this the gfn sota? Then point out this fact.\n\n \n\n \n\nThe tree topology can be either a rooted binary tree or a bifurcating unrooted tree.\n\n■ Page 2\n\n \n\nThis is a potential weakness. Can you restrict you method to binary trees? Does this mean that the posterior support of a subsplit consist of both binary and multifurcating trees? On an earlier reading i got the impression that you considered binary trees. Please clarify.\n\n \n\n \n\nThis equation can be used to recursively compute 𝑃(𝐿𝑖 𝑢|𝑎𝑖 𝑢) at all nodes of the tree and sites 𝑖. The algorithm performs a post-order traversal of the tree,\n\n■ Page 3\n\n \n\nThis is called dynamic programming\n\n \n\n \n\nrecursively computed by (1)\n\n■ Page 3\n\n \n\nAgain DP\n\n \n\n \n\nThe algorithm traverses the tree two times, first in post-order (bottom-up) to calculate the character set at each node, then in pre- order (top-down) to assign optimal sequences.\n\n■ Page 4\n\n \n\nAlso this is DP\n\n \n\n \n\nGenerative flow networks (GFlowNets) are algorithms for learning generative models of complex distributions given by unnormalized density functions over structured spaces. Here, we give a con- cise summary of the the GFlowNet framework.\n\n■ Page 4\n\n \n\nThis section could provide more insight\n\n \n\n \n\nThe direct optimization of 𝑃𝐹’s parameters 𝜃 is impossible since it involves an intractable sum over all complete trajectories.\n\n■ Page 4\n\n \n\nThis argument is incorrect. With polynomial length trajectories, which you need, exponentially many trajectories implies exponentially many states. The ladder also yields your optimization infeasible, in worst case.\n\n \n\n \n\nEach action chooses a pair of trees and join them at the root, thus creating a new tree\n\n■ Page 5\n\n \n\nThis suggests binary trees. “At the root” is a poor formulation.\n\n \n\n \n\nhus, a state 𝑠 consists of a set of rooted trees\n\n■ Page 5\n\n \n\nDisjoint rooted …\n\n \n\n \n\nits two children\n\n■ Page 5\n\n \n\nThis should close the case!\n\n \n\n \n\nall of which can be reached by our sampler\n\n■ Page 5\n\n \n\nAre all probabilities always non zero ?\n\n \n\n \n\nState representation To represent a rooted tree in a non-terminal state, we compute features for each site independently by taking advantage of the Felsenstein features (§3.1.1).\n\n■ Page 5\n\n \n\nCan you motivate this choice?\n\n \n\n \n\nAlthough the proposed feature representation 𝜌 does not capture all the information of tree structure and leaf sequences, we show that 𝜌 indeed contains sufficient informa- tion to express the optimal policy\n\n■ Page 6\n\n \n\nWhat is the optimal policy? How do you show this?\n\n \n\n \n\nProposition 1. Let 𝑠1 = {(𝑧1, 𝑏1), (𝑧2, 𝑏2) . . . (𝑧𝑙, 𝑏𝑙)} and 𝑠2 = {(𝑧′ 1, 𝑏′ 1), (𝑧′ 2, 𝑏2) . . . (𝑧′ 𝑙, 𝑏𝑙)} be two non-terminal states such that 𝑠1 ≠ 𝑠2 but sharing the same features 𝜌𝑖 = 𝜌′ 𝑖. Let 𝒂 be any sequence of actions, which applied to 𝑠1 and 𝑠2, respectively, results in full weighted trees 𝑥 = (𝑧, 𝑏𝑧 ), 𝑥′ = (𝑧′, 𝑏′), with two partial trajectories 𝜏 = (𝑠1 → · · · → 𝑥), 𝜏′ = (𝑠2 → · · · → 𝑥′). If 𝑃𝐹 is the policy of an optimal GFlowNet with uniform 𝑃𝐵, then 𝑃𝐹 (𝜏) = 𝑃𝐹 (𝜏′\n\n■ Page 6\n\n \n\nPlease explain the importance.\n\n \n\n \n\n−7108.42 ±0.18 −7108.41 ±0.14 −7290.36 ±7.23 −7111.55 ±0.07 −7108.95 ±0.06\n\n■ Page 7\n\n \n\nAt least the second method is a lower bound and it entire interval os above yours. Isn’t that better?\n\n \n\n \n\nPearson correlation of model sampling log-density and ground truth unnormalized posterior log-density for each dataset o\n\n■ Page 8\n\n \n\nHow is the ground truth obtained?\n\n \n\n \n\n589\n\n■ Page 8\n\n \n\n \n\n512\n\n■ Page 8\n\n \n\n \n\n0.624\n\n■ Page 8\n\n \n\nFails miserably\n\n \n\n \n\nBAYESIAN PHYLOGENETIC INFERENCE\n\n■ Page 8\n\n \n\nWhat is the computational resources required for these methods?\n\n \n\n \n\nThe base- lines we compare to are the MCMC-based MrBayes combined with the stepping-stone sampling tech- nique (Ronquist et al., 2012),\n\n■ Page 8\n\n \n\nHow do you obtain that and are you sure of how it is computed?\n\n \n\n \n\nTo select pairs of trees to join, we evaluate tree-pair features for every pair of trees in the state and pass these tree-pair features as input to the tree MLP to generate probability logits for all pairs of trees. The tree-pair feature for a tree pair (𝑖, 𝑗) with representations 𝑒𝑖, 𝑠𝑗 is the concatenation of 𝑒𝑖 + 𝑒 𝑗 with the summary embedding of the state, i.e., the feature is [𝑒𝑠; 𝑒𝑖 + 𝑒 𝑗], where [·; ·] denotes vector direct sum (concatenation). For a state with 𝑙 trees, 𝑙 2 = 𝑙(𝑙−1) 2 such pairwise features are generated for all possible pairs\n\n■ Page 16", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The strength is really the novelty and that this method both obtains good performance and is not limited to a predefined set of trees.", "weaknesses": "It is not really trying to explain the GFN approach to the reader, but rather enumerating required items. The running time and performance in general require further explanation. The parsimony case is clearly less interesting and could partially be move to the appendix.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (poster)