Probabilistic Path Hamiltonian Monte Carlo
Vu Dinh Arman Bilge Cheng Zhang Frederick A. Matsen IV


[
Probabilistic Path Hamiltonian Monte Carlo















equal*


Vu Dinhequal,fredhutch
Arman Bilgeequal,fredhutch,UW
Cheng Zhangequal,fredhutch
Frederick A. Matsen IVfredhutch


fredhutchProgram in Computational Biology, Fred Hutchison Cancer Research Center, Seattle, WA, USA
UWDepartment of Statistics, University of Washington, Seattle, WA, USA

Frederick, A. Matsen IVmatsen@fredhutch.org




Hamiltonian Monte Carlo, leapfrog algorithm, phylogenetics

0.3in
]










 









Hamiltonian Monte Carlo (HMC) is an efficient and effective means of sampling posterior distributions on Euclidean space, which has been extended to manifolds with boundary.
However, some applications require an extension to more general spaces.
For example, phylogenetic (evolutionary) trees are defined in terms of both a discrete graph and associated continuous parameters; although one can represent these aspects using a single connected space, this rather complex space is not suitable for existing HMC algorithms.
In this paper, we develop Probabilistic Path HMC (PPHMC) as a first step to sampling distributions on spaces with intricate combinatorial structure.
We define PPHMC on orthant complexes, show that the resulting Markov chain is ergodic, and provide a promising implementation for the case of phylogenetic trees in open-source software.
We also show that a surrogate function to ease the transition across a boundary on which the log-posterior has discontinuous derivatives can greatly improve efficiency.







§ INTRODUCTION


Hamiltonian Monte Carlo is a powerful sampling algorithm which has been shown to outperform many existing MCMC algorithms, especially in problems with high-dimensional and correlated distributions <cit.>.
The algorithm mimics the movement of a body balancing potential and kinetic energy by extending the state space to include auxiliary momentum variables and using Hamiltonian dynamics.
By traversing long iso-probability contours in this extended state space, HMC is able to move long distances in state space in a single update step, and thus has proved to be more effective than standard MCMC methods in a variety of applications.
The method has gained a lot of interest from the scientific community and since then has been extended to tackle the problem of sampling on various geometric structures such as constrained spaces <cit.>, on general Hilbert space <cit.> and on Riemannian manifolds <cit.>.

However, these extensions are not yet sufficient to apply to all sampling problems, such as in phylogenetics, the inference of evolutionary trees.
Phylogenetics is the study of the evolutionary history and relationships among individuals or groups of organisms.
In its statistical formulation it is an inference problem on hypotheses of shared history based on observed heritable traits under a model of evolution.
Phylogenetics is an essential tool for understanding biological systems and is an important discipline of computational biology.
The Bayesian paradigm is now commonly used to assess support for inferred tree structures or to test hypotheses that can be expressed in phylogenetic terms <cit.>.

Although the last several decades have seen an explosion of advanced methods for sampling from Bayesian posterior distributions, including HMC, phylogenetics still uses relatively classical Markov chain Monte Carlo (MCMC) based methods.
This is in part because the number of possible tree topologies (the labeled graphs describing the branching structure of the evolutionary history) explodes combinatorially as the number of species increases.
Also, to represent the phylogenetic relation among a fixed number of species, one needs to specify both the tree topology (a discrete object) and the branch lengths (continuous distances).
This composite structure has thus far limited sampling methods to relatively classical Markov chain Monte Carlo (MCMC) based methods.
One path forward is to use a construction of the set of phylogenetic trees as a single connected space composed of Euclidean spaces glued together in a combinatorial fashion <cit.> and try to define an HMC-like algorithm thereupon.



Experts in HMC are acutely aware of the need to extend HMC to such spaces with intricate combinatorial structure: <cit.> describes the extension of HMC to discrete and tree spaces as a major outstanding challenge for the area.
However, there are several challenges to defining a continuous dynamics-based sampling methods on such spaces.
These tree spaces are composed of Euclidean components, one for each discrete tree topology, which are glued together in a way that respects natural similarities between trees.
These similarities dictate that more than two such Euclidean spaces should get glued together along a common lower-dimensional boundary.
The resulting lack of manifold structure poses a problem to the construction of an HMC sampling method on tree space, since up to now, HMC has just been defined on spaces with differential geometry.
Similarly, while the posterior function is smooth within each topology, the function's behavior may be very different between topologies.
In fact, there is no general notion of differentiability of the posterior function on the whole tree space and any scheme to approximate Hamiltonian dynamics needs to take this issue into account.

In this paper, we develop Probabilistic Path Hamiltonian Monte Carlo (PPHMC) as a first step to sampling distributions on spaces with intricate combinatorial structure (Figure <ref>).
After reviewing how the ensemble of phylogenetic trees is naturally described as a geometric structure we identify as an orthant complex <cit.>, we define PPHMC for sampling posteriors on orthant complexes along with a probabilistic version of the leapfrog algorithm.
This algorithm generalizes previous HMC algorithms by doing classical HMC on the Euclidean components of the orthant complex, but making random choices between alternative paths available at a boundary.
We establish that the integrator retains the good theoretical properties of Hamiltonian dynamics in classical settings, namely probabilistic equivalents of time-reversibility, volume preservation, and accessibility, which combined together result in a proof of ergodicity for the resulting Markov chain.
Although a direct application of the integrator to the phylogenetic posterior does work, we obtain significantly better performance by using a surrogate function near the boundary between topologies to control approximation error.
This approach also addresses a general problem in using Reflective HMC <cit.> for energy functions with discontinuous derivatives (for which the accuracy of RHMC is of order 𝒪(ϵ), instead of the standard local error 𝒪(ϵ^3) of HMC on ℝ^n).
We provide, validate, and benchmark two independent implementations in open-source software.




§ MATHEMATICAL FRAMEWORK




 §.§ Bayesian learning on phylogenetic tree space


A phylogenetic tree (τ, q) is a tree graph τ with N leaves, each of which has a label, and such that each edge e is associated with a non-negative number q_e.
Trees will be assumed to be bifurcating (internal nodes of degree 3) unless otherwise specified.
We denote the number of edges of such a tree by n=2N-3.
Any edge incident to a leaf is called a pendant edge, and any other edge is called an internal edge.
Let 𝒯_N be the set of all N-leaf phylogenetic trees for which the lengths of pendant edges are bounded from below by some constant e_0 > 0.
(This lower bound on branch lengths is a technical condition for theoretical development and can be relaxed in practice.)

We will use nearest neighbor interchange (NNI) moves <cit.> to formalize what tree topologies that are “near” to each other.
An NNI move is a transformation that collapses an internal edge to zero and then expands the resulting degree 4 vertex into an edge and two degree 3 vertices in a new way (Figure <ref>a).
Two tree topologies τ_1 and τ_2 are called adjacent topologies if there exists a single NNI move that transforms τ_1 into τ_2.

We will parameterize 𝒯_N as Billera-Holmes-Vogtmann (BHV) space <cit.>, which we describe as follows.
An orthant of dimension n is simply ℝ_≥ 0^n; each n-dimensional orthant is bounded by a collection of lower dimensional orthant faces.
An orthant complex is a geometric object 𝒳 obtained by gluing various orthants of the same dimension n, indexed by a countable set Γ, such that: (i) the intersection of any two orthants is a face of both orthants, and (ii) each x ∈𝒳 belongs to a finite number of orthants.
Each state of 𝒳 is thus represented by a pair (τ, q), where τ∈Γ and q ∈ℝ_≥ 0^n.
Generalizing the definitions from phylogenetics, we refer to τ as its topology and to q as the vector of attributes.
The topology of a point in an orthant complex indexes discrete structure, while the attributes formalize the continuous components of the space.

For phylogenetics, the complex is constructed by taking one n-dimensional orthant for each of the (2n-3)!! possible tree topologies, and gluing them together along their common faces.
The geometry can also be summarized as follows.
In BHV space, each of these orthants parameterizes the set of branch lengths for a single topology (as a technical point, because we are bounding pendant branch lengths below by e_0, we can take the corresponding entries in the orthant to parameterize the amount of branch length above e_0).
Top-dimensional orthants of the complex sharing a facet, i.e. a codimension 1 face, correspond to (NNI) adjacent topologies.

For a fixed phylogenetic tree (τ,q), the phylogenetic likelihood is defined as follows and will be denoted by L(τ, p) <cit.>.
Let ψ = (ψ_1, ψ_2,...,ψ_S) ∈Ω^N × S be the observed sequences (with characters in Ω) of length S over N leaves.
The likelihood of observing ψ given the tree has the form



    L(τ, q) = ∏_s=1^S∑_a^sη(a_ρ^s)∏_(u,v)∈ E(τ, q)P^uv_a^s_ua^s_v( q_uv)

where ρ is any internal node of the tree, a^s ranges over all extensions of ψ^s to the internal nodes of the tree, a^s_u denotes the assigned state of node u by a^s, E(τ, q) denotes the set of tree edges, P_ij(t) denotes the transition probability from character i to character j across an edge of length t defined by a given evolutionary model and η is the stationary distribution of this evolutionary model.
For this paper we will assume the simplest <cit.> model of a homogeneous continuous-time Markov chain on Ω with equal transition rates, noting that inferring parameters of complex substitution models is a vibrant yet separate subject of research <cit.>.

Given a proper prior distribution with density π_0 imposed on the branch lengths and on tree topologies, the posterior distribution can be computed as 𝒫(τ, q) ∝ L(τ, q) π_0(τ, q).





 §.§ Bayesian learning on orthant complexes


With the motivation of phylogenetics in mind, we now describe how the phylogenetic problem sits as a specific case of a more general problem of Bayesian learning on orthant complexes, and distill the criteria needed to enable PPHMC on these spaces.
This generality will also enable applications of Bayesian learning on similar spaces in other settings.
For example in robotics, the state complex can be described by a cubical complexes whose vertices are the states, whose edges correspond to allowable moves, and whose cubes correspond to collections of moves which can be performed simultaneously <cit.>.
Similarly, in learning on spaces of treelike shapes, the attributes are open curves translated to start at the origin, described by a fixed number of landmark points <cit.>.

An orthant complex, being a union of Euclidean orthants, naturally inherits the Lebesgue measure which we will denote hereafter by μ.

Orthant complexes are typically not manifolds, thus to ensure consistency in movements across orthants, we assume that local coordinates of the orthants are defined in such a way that there is a natural one-to-one correspondence between the sets of attributes of any two orthants sharing a common face.

[Consistency of local coordinates]

Given two topologies τ, τ' ∈Γ and state x = (τ, q_τ) = (τ', q_τ') on the boundary of the orthants for τ and τ', then q_τ= q_τ'.


We show that BHV tree space can be given such coordinates in the Appendix.
For the rest of the paper, we define for each state (τ, q) ∈𝒳 the set 𝒩(τ, q) of all neighboring topologies τ' such that τ' orthant contains (τ, q).
Note that 𝒩(τ, q) always includes τ, and if all coordinates of q are positive, 𝒩(τ, q) is exactly {τ}.
Moreover, if τ' ∈𝒩(τ, q) and τ' τ, we say that τ and τ' are joined by (τ,q).
If the intersection of orthants for two topologies is a facet of each, we say that the two topologies are adjacent.

Finally, let 𝒢 be the adjacency graph of the orthant complex 𝒳, that is, the graph with vertices representing the topologies and edges connecting adjacent topologies.
Recalling that the diameter of a graph is the maximum value of the graph distance between any two vertices, we assume that


The adjacency graph 𝒢 of 𝒳 has finite diameter, hereafter denoted by k.


For phylogenetics, k is of order 𝒪(N log N) <cit.>.

We seek to sample from a posterior distribution 𝒫(τ, q) on 𝒳.
Assume that the negative logarithm of the posterior distribution U(τ, q) :=-log P(τ, q) satisfies:



U(τ, q) is a continuous function on 𝒳, and is smooth up to the boundary of each orthant τ∈Γ.


In the Appendix, we prove that if the logarithm of the phylogenetic prior distribution π_0(τ, q) satisfies Assumption <ref>, then so does the phylogenetic posterior distribution.
It is also worth noting that while U(τ, q) is smooth within each orthant, the function's behavior may be very different between orthants and we do not assume any notion of differentiability of the posterior function on the whole space.



 §.§ Hamiltonian dynamics on orthant complexes

The HMC state space includes auxiliary momentum variables in addition to the usual state variables.
In our framework, the augmented state of this system is represented by the position (τ, q) and the momentum p, an n-dimensional vector.
We will denote the set of all possible augmented state (τ, q, p) of the system by 𝕋.

The Hamiltonian is defined as in the traditional HMC setting: H(τ, q, p) = U(τ, q) + K(p),
where K(p) = 1/2p^2.
We will refer to U(τ, q) and K(p) as the potential energy function and the kinetic energy function of the system at the state (τ,q, p), respectively.

Our stochastic Hamiltonian-type system of equations is:

    dp_i/dt =    - ∂ U/∂ q_i(τ, q)         if   q_i > 0 
    
    p_i  - p_i;      τ∼  Z(𝒩(τ, q))         if   q_i = 0  
    dq_i/dt =     p_i

where Z(A) denotes the uniform distribution on the set A.

If all coordinates of q are positive, the system behaves as in the traditional Hamiltonian setting on ℝ^n.
When some attributes hit zero, however, the new orthant is picked randomly from the orthants of neighboring topologies (including the current one), and the momenta corresponding to non-positive attributes are negated.

Assumption <ref> implies that despite the non-differentiable changes in the governing equation across orthants, the Hamiltonian of the system along any path is constant:

 H is conserved along any system dynamics.
 





 §.§ A probabilistic “leap-prog” algorithm

In practice, we approximate Hamiltonian dynamics by the following integrator with step size ϵ, which we call “leap-prog” as a probabilistic analog of the usual leapfrog algorithm.
This extends previous work of <cit.> on RHMC where particles can reflect against planar surfaces of ℝ^n_≥ 0.

In the RHMC formulation, one breaks the step size ϵ into smaller sub-steps, each of which correspond to an event when some of the coordinates cross zero.
We adapt this idea to HMC on orthant complexes as follows.
Every time such an event happens, we reevaluate the values of the position and the momentum vectors, update the topology (uniformly at random from the set of neighboring topologies), reverse the momentum of crossing coordinates and continue the process until a total step size ϵ is achieved (Algorithm <ref>).
We note that several topologies might be visited in one leap-prog step.

If there are no topological changes in the trajectory to time ϵ, this procedure is equivalent to classical HMC.
Moreover, since the algorithm only re-evaluates the gradient of the energy function at the end of the step when the final position has been fixed, changes in topology on the path have no effect on the changes of position and momentum.
Thus, the projection of the particles (in a single leap-prog step) to the (q,p) space is identical to a leapfrog step of RHMC on ℝ^n_≥ 0.



Here FirstUpdateEvent(τ, q,p, t) returns x, the position of the first event for which the line segment [q, q+ tp] crosses zero;  e, the time when this event happens; and I, the indices of the coordinates crossing zero during this event.
If q_i and p_i are both zero before FirstUpdateEvent is called, i is not considered as a crossing coordinate.
If no such an event exists, ∅ is returned.




§ HAMILTONIAN MONTE CARLO ON ORTHANT COMPLEXES


Probabilistic Path Hamiltonian Monte Carlo (PPHMC) with leap-prog dynamics iterates three steps, similar to that of classical HMC.
First, new values for the momentum variables are randomly drawn from their Gaussian distribution, independently of the current values of the position variables.
Second, starting with the current augmented state, s=(τ, q, p), the Hamiltonian dynamics is run for a fixed number of steps T using the leap-prog algorithm with fixed step size ϵ.
The momentum at the end of this trajectory is then negated, giving a proposed augmented state s^* = (τ^*, q^*, p^*).
Finally, this proposed augmented state is accepted as the next augmented state of the Markov chain with probability r(s,s^*)= min(1, exp(H(s) - H(s^*))).

PPHMC has two natural advantages over MCMC methods for phylogenetic inference: jumps between topologies are guided by the potential surface, and many jumps can be combined into a single proposal with high acceptance probability.
Indeed, rather than completely random jumps as used for MCMC, the topological modifications of HMC are guided by the gradient of the potential.
This is important because there are an enormous number of phylogenetic trees, namely (2n-3)!! trees with n leaves.
Secondly, HMC can combine a great number of tree modifications into a single step, allowing for large jumps in tree space with high acceptance probability.
These two characteristics are analogs of why HMC is superior to conventional MCMC in continuous settings and what we aimed to extend to our problem.



 §.§ Theoretical properties of the leap-prog integrator

To support the use of this leap-prog integrator for MCMC sampling, we establish that integrator retains analogs of the good theoretical properties of Hamiltonian dynamics in classical settings, namely, time-reversibility, volume preservation and accessibility (proofs in Appendix).

We formulate probabilistic reversibility as:
[Reversibility]
For a fixed finite time horizon T, we denote by P(s, s') the probability that the integrator moves s to s' in a single update step.
We have

    P((τ, q, p), (τ', q', p') ) =P((τ', q', -p'), (τ, q, -p) ).

for any augmented states (τ', q',p') and (τ, q, p) ∈𝕋.



The central part of proving the detailed balance condition of PPHMC is to verify that Hamiltonian dynamics preserves volume.
Unlike the traditional HMC setting where the proposal distribution is a single point mass, in our probabilistic setting, if we start at one augmented state s, we may end up at countably many end points due to stochastic HMC integration.
The equation describing volume preservation in this case needs to be generalized to the form of Equation (<ref>), where the summations account for the discreteness of the proposal distribution.
[Volume preservation]
For every pair of measurable sets A,B ⊂𝕋 and elements s,s' ∈𝕋, we denote by P(s,s') the probability that the integrator moves s to s' in a single update step and define

    B(s) = {  s' ∈ B: P(s, s') >0 }

and

    A(s') = {  s ∈ A: P(s, s') >0 }.

Then

    ∫_A∑_s' ∈ B(s)P(s,s') ds = ∫_B∑_s ∈ A(s')P(s',s) ds'.



If we restrict to the case of trajectories staying in a single topology, A(s) and B(s') are singletons and we get back the traditional equation of volume preservation.
We also note that the measure ds in Equation (<ref>) is the Lebesgue measure: when there is no randomness in the Hamiltonian paths, (3.1) becomes the standard volume preservation condition, where volumes are expressed by the Lebesgue measure.

Typically, accessibility poses no major problem in various settings of HMC since it is usually clear that one can go between any two positions in a single HMC step.
In the case of PPHMC, however, the composition of discrete and continuous structure, along with the possible non-differentiability of the potential energy across orthants, make it challenging to verify this condition.
Here we show instead that the PPHMC algorithm can go between any two states with k steps, where k denotes the diameter of adjacency graph 𝒢 the space 𝒳 and each PPHMC step consists of T leap-prog steps of size ϵ.

[k-accessibility]
For a fixed starting state (τ^(0), q^(0)), any state (τ', q') ∈𝒳 can be reached from (τ^(0), q^(0)) by running k steps of PPHMC.



The proof of this Lemma is based on Assumption <ref>, which asserts that the adjacency graph 𝒢 of 𝒳 has finite diameter, and that classical HMC allows the particles to move freely in each orthant by a single HMC step.

To show that Markov chains generated by PPHMC are ergodic, we also need to prove that the integrator can reach any subset of positive measure of the augmented state space with positive probability.
To enable such a result, we show:


For every sequence of topologies ω = {τ^(0), τ^(1), …, τ^(n_ω)} and every set with positive measure B ⊂𝒳, let B_ω be the set of all (τ', q') ∈ B such that (τ', q') can be reached from (τ^(0), q^(0)) in k PPHMC steps and such that the sequence of topologies crossed by the trajectory is ω.
We denote by I_B, ω the set of all sequences of initial momenta for each PPHMC step {p^(0), …, p^(k)} that make such a path possible.

Then, if μ(I_B, ω)=0, then μ(B_ω)=0.



We also need certain sets to be countable.



Given s ∈𝕋, we denote by R(s) the set of all augmented states s' such that there is a finite-size leap-prog step with path γ connecting s and s', and by K(s) the set of all such leap-prog paths γ connecting s and s' ∈ R(s).
Then R(s) and K(s) are countable.
Moreover, the probability P_∞( s, s') of moving from s to s' via paths with infinite number of topological changes is zero.






 §.§ Ergodicity of Probabilistic Path HMC


In this section, we establish that a Markov chain generated by PPHMC is ergodic with stationary distribution π(τ, q) ∝exp(-U(τ, q)).
To do so, we need to verify that the Markov chain generated by PPHMC is aperiodic, because we have shown k-accessibility of the integrator rather than 1-accessibility.
Throughout this section, we will use the notation P((τ, q, p), ·) to denote the one-step proposal distribution of PPHMC starting at augmented state (τ, q, p), and P((τ, q), ·) to denote the one-step proposal distribution of PPHMC starting at position (τ, q) and with a momentum vector drawn from a Gaussian as described above.

We first note that:

PPHMC preserves the target distribution π.


Given probabilistic volume preservation (<ref>), the proof is standard and is given in the Appendix.

[Ergodic]

The Markov chain generated by PPHMC is ergodic.



For every sequence of topologies ω =  {τ^(0), τ^(1), …, τ^(n_ω)} (finite by Lemma <ref>) and every set with positive measure B ⊂𝒳, we define B_ω and I_B, ω as in the previous section.
By Lemma <ref>, we have

    B = ⋃_ωB_ω.

Assume that μ(I_B, ω) = 0 for all ω.
From Lemma <ref>, we deduce that  μ(B_ω)=0 for all ω.
This makes μ(B) = 0, which is a contradiction.
Hence μ(I_B, ω) > 0 for some ω and P^n_ω((τ^(0),q^(0)), B ) is at least the positive quantity

    1/Z∫_p ∈ I_B, ωP^n_ω((τ^(0), q^(0), p),  B) exp (-K(p))dp

where Z is the normalizing constant.
This holds for all sets with positive measure B ⊂𝒳, so PPHMC is irreducible.

Now assume that a Markov chain generated by the leapfrog algorithm is periodic.
The reversibility of Hamiltonian dynamics implies that the period d must be equal to 2.
In other words, there exist two disjoint subsets X_1, X_2 of 𝒳 such that π(X_1)>0, and

    P(x, X_2)=1  ∀ x ∈ X_1,    and  P(x, X_1)=1  ∀ x ∈ X_2.

Consider x ∈ X_1 with all positive attributes.
There exists a neighborhood U_x around x such that any y ∈ U_x is reachable from x by Hamiltonian dynamics.
Since X_1, X_2 are disjoint, we deduce that μ(U_x ∩ X_1)=0.
Since the neighborhood U_x exists for almost every x ∈ X_1, this implies that μ(X_1)=0, and hence, that π(X_1)=0, which is a contradiction.
We conclude that any Markov chain generated by the leapfrog algorithm is aperiodic.

Lemma <ref> shows that PPHMC preserves the target distribution π.
This, along with π-irreducibility and aperiodicity, completes the proof <cit.>.




 §.§ An efficient surrogate smoothing strategy


One major advantage of HMC methods over traditional approaches is that HMC-proposed states may be distant from the current state but nevertheless have a high probability of acceptance.
This partially relies on the fact that the leapfrog algorithm with smooth energy functions has a local approximation error of order 𝒪(ϵ^3) (which leads to global error 𝒪(Tϵ^3), where T is the number of leapfrog steps in a Hamiltonian path).

However, when the potential energy function U(τ, q) is not differentiable on the whole space this low approximation error can break down.
Indeed, although PPHMC inherits many nice properties from vanilla HMC and RHMC <cit.>, this discontinuity of the derivatives of the potential energy across orthants may result in non-negligible loss of accuracy during numerical simulations of the Hamiltonian dynamics.
A careful analysis of the local approximation error of RHMC for potential energy functions with discontinuous first derivatives reveals that it only has an local error rate of order at least Ω(ϵ) (see proof in Appendix):


Given a potential function V, we denote by V^+ and V^- the restrictions of V on the half-spaces {x_1 ≥ 0} and {x_1 ≤ 0} and assume that V^+ and V^- are smooth up to the boundary of their domains.
If the first derivative with respect to the first component of the potential energy V(q) are discontinuous across the hyper-plane {x_1=0} (i.e., (∂ V^+)/(∂ q_1) and (∂ V^-)/(∂ q_1) are not identical on this set), then RHMC on this hyper-plane has a local error of order at least Ω(ϵ).



Since PPHMC uses RHMC, when the first derivatives of the potential energy are discontinuous, it also has a global error of order 𝒪(C ϵ + Tϵ^3), which depends on the number of critical events C along a Hamiltonian path (that is, the number of reflection/refraction events).
This makes it difficult to tune the step size ϵ for optimal acceptance rate, and requires small ϵ, limiting topology exploration.

To alleviate this issue, we propose the use of a surrogate induced Hamiltonian dynamics <cit.> with the Hamiltonian H̃(τ,q,p) = Ũ(τ,q) + K(p), where the surrogate potential energy is

    Ũ(τ,q)  = U(τ,G(q)),   G(q) = (g(q_1),…,g(q_n))

and g(x) is some positive and smooth approximation of |x| with vanishing gradient at x=0.
One simple example which will be used for the rest of this paper is

    g_δ(x) = {[             x,            x≥δ; 1/2δ(x^2+δ^2),        0≤ x <δ ].

where δ will be called the smoothing threshold.

Due to the vanishing gradient of g, Ũ now has continuous derivatives across orthants.
However, Ũ is no longer continuous across orthants since g(0)≠ 0 and we thus employ the refraction technique introduced in <cit.> (see Algorithm <ref> for more details).
The proposed state s_δ^∗ = (τ_δ^∗,q_δ^∗,p_δ^∗) at the end of the trajectory is accepted with probability according to the original Hamiltonian, that is, min(1,exp(H(s)-H(s_δ^∗))).

By following the same framework proposed in previous sections, we can prove that the resulting sampler  still samples from the exact posterior distribution 𝒫(τ,q).
A complete treatment, however, requires more technical adjustments and is beyond the scope of the paper.

We will leave this as a subject of future work.



As we will illustrate later, compared to the exact potential energy, the continuity of the derivative of the surrogate potential across orthants dramatically reduces the discretization error and allows for high acceptance probability with relatively large step size.

















§ EXPERIMENTS

In this section, we demonstrate the validity and efficiency of our PPHMC method by an application to Bayesian phylogenetic inference.
We compared our PPHMC implementations to industry-standard MrBayes 3.2.5, which uses MCMC to sample phylogenetic trees <cit.>.
We concentrated on the most challenging part: sampling jointly for the branch lengths and tree topologies, and assumed other parameters (e.g., substitution model, hyper-parameters for the priors) are fixed.
More specifically, for all of our experiments we continued to assume the Jukes-Cantor model of DNA substitution and placed a uniform prior on the tree topology τ∼ Z(𝒯_N) with branch lengths i.i.d. q_i ∼Exponential(λ = 10), as done by others when measuring the performance of MCMC algorithms for Bayesian phylogenetics <cit.>.
As mentioned earlier, although in the theoretical development we assumed that the lengths of the pendant edges are bounded from below by a positive constant e_0 to ensure that the likelihood stays positive on the whole tree space, this condition is not necessary in practice since the Hamiltonian dynamics guide the particles away from regions with zero likelihood (i.e., the region with U = ∞).
We validate the algorithm through two independent implementations in open-source software:


  * a Scala version available at <https://github.com/armanbilge/phyloHMC> that uses the Phylogenetic Likelihood Library[<https://github.com/xflouris/libpll>] <cit.>, and

  * a Python version available at <https://github.com/zcrabbit/PhyloInfer> that uses the ETE toolkit <cit.> and Biopython <cit.>.




 §.§ Simulated data

As a proof of concept, we first tested our PPHMC method on a simulated data set.
We used a random unrooted tree with N=50 leaves sampled from the aforementioned prior.
1000 nucleotide observations for each leaf were then generated by simulating the continuous-time Markov model along the tree.
This moderate data set provided enough information for model inference while allowing for a relatively rich posterior distribution to sample from.

We ran MrBayes for 10^7 iterations and sampled every 1000 iterations after a burn-in period of the first 25% iterations to establish a ground truth for the posterior distribution.
For PPHMC, we set the step size ϵ=0.0015 and smoothing threshold δ=0.003 to give an overall acceptance rate of about α = 0.68 and set the number of leap-prog steps T=200.
We then ran PPHMC for 10,000 iterations with a burn-in of 25%.
We saw that PPHMC indeed samples from the correct posterior distribution (see Figure <ref> in the Appendix).










 §.§ Empirical data


We also analyzed an empirical data set labeled DS4 by <cit.> that has become a standard benchmark for MCMC algorithms for Bayesian phylogenetics since <cit.>.
DS4 consists of 1137 nucleotide observations per leaf from N=41 leaves representing different species of fungi.
Notably, only 554 of these observations are complete; the remaining 583 are missing a character for one or more leaves so the likelihood is marginalized over all possible characters.
<cit.> observed that the posterior distribution for DS4 features high-probability trees separated by paths through low-probability trees and thus denoted it a “peaky” data set that was difficult to sample from using MrBayes.

To find the optimal choice of tuning parameters for DS4, we did a grid search on the space of step size ϵ and the smoothing threshold–step size ratio δ / ϵ.
The number of leap-prog steps T was adjusted to keep the total simulation time ϵ T fixed.
For each choice of parameters, we estimated the expected acceptance rate α by averaging over 20 proposals from the PPHMC transition kernel per state for 100 states sampled from the posterior in a previous, well-mixed run.
This strategy enabled us to obtain an accurate estimate of the average acceptance rate without needing to account for the different mixing rates in a full PPHMC run due to the various settings for the tuning parameters.

The results suggest that choosing δ≈ 2ϵ maximizes the acceptance probability (Figure <ref>b).
Furthermore, when aiming for the optimal acceptance rate of α = 0.65 <cit.>, the use of the surrogate function enables a choice of step size ϵ nearly 10 times greater than otherwise.
In practice, this means that an equivalent proposal requires less leap-prog steps and gives a more efficient sampling algorithm.

To see the difference this makes in practice, we ran long trajectories for exact and surrogate-smoothed PPHMC with a relatively large step size ϵ=0.0008.
Indeed, we found that the surrogate enables very long trajectories and large number of topology transformations (Figure <ref>a,c).


























§ CONCLUSION

Sophisticated techniques for sampling posteriors using HMC have thus far been restricted to manifolds with boundary.
To address this limitation, we have developed “PPHMC,” which is the first extension of HMC to a space with intricate combinatorial structure.
The corresponding integrator makes a random choice among alternatives when encountering a boundary.
To prove ergodicity, we extend familiar elements of HMC proofs to this probabilistic path setting.
We develop a smoothing surrogate function that enables long HMC paths with many boundary transitions across which the posterior is not differentiable.
Our surrogate method enables high acceptance probability for RHMC <cit.> in the case of potential functions with discontinuous derivatives; this aspect of our work is independent of the probabilistic nature of PPHMC.
Our implementation shows good performance on both simulated and real data.
There are many opportunities for future development, including extending the theory to other classes of combinatorially-described spaces and surrogate functions, developing adaptive path length algorithms, as well as extending our implementation for phylogenetics to sample other mutation model and demographic parameters along with topologies.










plainnat










§ APPENDIX




 §.§ Properties of the phylogenetic posterior distribution






Recall that L(τ, q) denotes the likelihood function of the tree T=(τ, q), we have

    U(τ, q) = - log L(τ, q) - logπ_0(τ, q)

Since -logπ_0(τ, q) is assumed to satisfy the Assumption <ref>, we just need to prove that the phylogenetic likelihood function is smooth while each orthant and is continuous on the whole space.

Without loss of generality, we consider the case when a single branch length of some edge e is contracted to zero.
To investigate the changes in the likelihood function and its derivatives, we first fix all other branches, partition the set of all extensions of ψ according to their labels at the end points of e, and split E(T) into two sets of edges E_left and E_right corresponding to the location of the edges with respect to e.
The likelihood function of the tree T=(τ, q) can be rewritten as

    L(T)    =   ∏_s=1^S∑_ij∑_a ∈𝒜_ij( ∏_(u,v) ∈ E_leftP^uv_a_ua_v( q_uv)) 
             × η(i) P^e_ij(t) ×( ∏_(u,v) ∈ E_rightP^uv_a_ua_v( q_uv))

where t is the branch length of e, η is the stationary distribution, 𝒜_ij denotes the set of all extensions of ψ for which the labels at the left end point and the right end point of e are i and j, respectively.
By grouping the products over E_left and E_right, the stationary frequency η(·), and the sum over a in a single term b_ij^s, we can define the one-dimensional log-likelihood function as a univariate function of the branch length of e

    L_T(t)=∏_s=1^S(∑_ijb^s_ijP^e_ij(t)).

Consider the tree T' obtained by collapsing edge e of the tree T to zero.
The likelihood of T' can be written as

    L(T') = ∏_s=1^S(∑_i=jb^s_ijP^e_ij(0)) =∏_s=1^S(∑_ib^s_ii)

since P_ij(0) = 1 if i=j and 0 otherwise.
Thus

    lim_t→ 0L_T(t) = L(T').

Since this is true for all (τ, q) and t ∈ E(τ, q), we deduce that the likelihood function is continuous up to the boundary of each orthant, and thus, is continuous on the whole tree space.
Moreover, using the same arguments, we can prove that likelihood function is smooth up to the boundary of each orthant.

Now fixing all but two branch lengths t_e, t_f, the likelihood can be rewritten as

    L_T(t_e, t_f)=∏_s=1^S(∑_ijb^s_ij(t_e)P^f_ij(t_f))

and the derivative of the log likelihood is

    1/L_T(t_e, t_f) ∂ L_T/∂ t_f(t_e, t_f)= ∑_s=1^S∑_ijb^s_ij(t_e)(P^f_ij)'(t_f)/∑_ijb^s_ij(t_e)P^f_ij(t_f).

By using the same argument as above, we have that b_ij^s (t_e) is continuous in t_e up to zero and so

    lim_t_e → 01/L_T(t_e, t_f) ∂ L_T/∂ t_f(t_e, t_f) = 1/L(T') ∂ L/∂ t_f(T').

Thus, when a Hamiltonian particle crosses a boundary between orthants, partial derivatives of the energy function with respect to positive branch lengths are continuous.





 §.§ Theoretical properties of the leap-prog integrator



Note that for PPHMC,  in a single leap-prog step γ of finite size ϵ, the algorithm only re-evaluates the gradient of the energy function at the end of the step when the final position q' has been fixed, and changes in topology on the path have no effect on the changes of position and momentum.
Thus, the projection γ̃ of γ to the (q,p) space is just a deterministic reflected Hamiltonian path.
As a result, for any s^(1) =(τ^(1), q^(1), p^(1)), s^(2) =(τ^(2), q^(2), p^(2)) ∈ R(s), we have (q^(1), p^(1))= (q^(2), p^(2)).
This, along with the fact that set of topologies is countable, implies that R(s) is countable.

Now denote by {t^(1) < t^(2) < … < t^(n) < …≤ϵ} the set of time points at which γ̃ hits the boundary.
Since this set is strictly increasing, it is countable.
Moreover, the τ-component of γ is only updated with finite choices at {t^(i)}.
This implies that K(s) is countable.

Finally, consider any leap-prog step γ that connects s and s' through infinite number of topological changes.
We note that at each t^(i), the next topology is chosen among x^(i)≥ 2 neighboring topologies.
Denote by P_γ( s, s') the probability of moving from s to s' via path γ, we have

    P_γ(s, s') ≤∏_i=1^∞1/x^(i)= 0.

Since K(s) is countable, we deduce that P_∞( s, s')=0.




Consider any possible path γ that connects s and s'. By definition, one can find a sequence of augmented states (s=s^(0), s^(1), s^(2), …, s^(k) = s' ) such that γ can be decomposed into segments on which the topology does not change.
From standard result about Hamiltonian dynamics, the Hamiltonian is constant on each segment.

For PPHMC, since the potential energy is continuous across the boundary and the magnitude of the momentum does not change when moving from one orthant to another one, we deduce that the Hamiltonian is constant along that path.

Similarly, for PPHMC with surrogates, the algorithm is designed in such a way that any changes in potential energy is balanced by a change in momentum, which conserves the total energy from one segment to another.
We also deduce that the Hamiltonian is constant along the whole path.




Define σ(τ, q, p):=(τ, q, -p).
Consider any possible leap-prog step γ that connects s and s'; say the sequences of augmented states (s=s^(0), s^(1), s^(2), …, s^(k) = s' ), topologies (τ=τ^(0), τ^(1), τ^(2), …, τ^(k) = τ' ) and times (t=t^(0), t^(1), t^(2), …,t^(k) = t' ) decompose γ into segments on which the topology is unchanged.
Denote by P_γ( s, s') the probability of moving from s to s' via path γ, we have

    P_γ( s, s')    = ∏_iℙ(s^(i+1)| s^(i), t^(i+1) -t^(i))
       × ∏_jℙ(τ^(j+1)|τ^(j)),

where each sub-step of the algorithm is a leapfrog update (ϕ^(i)) with some momentum reversing (σ^(i)), that is s^(i+1) = σ^(i)(ϕ^(i)(s^(i)))  and σ^(i) is a map that changes the sign of some momentum coordinates.

If we start the dynamics at σ(s^(i+1)), then since the particle is crossing the boundary, the momenta corresponding to the crossing coordinates are immediately negated and the system is instantly moved to the augmented state

    σ^(i)σ(s^(i+1)) = σσ^(i)(s^(i+1)) = σ(ϕ^(i)(s^(i))).

A standard result about reversibility of Reflective Hamiltonian dynamics implies that the system starting at σ(ϕ^(i)(s^(i))) will end at σ(s^(i)) after the same period of time t^(i+1)-t^(i).
We deduce that

    ℙ(s^(i+1)| s^(i), t^(i+1) -t^(i)))
    
    =    ℙ(σ(s^(i)) |σ(s^(i+1)), t^(i+1) -t^(i))).

On the other hand, at time t^(j), (τ^(j), q^(j)) and (τ^(j+1), q^(j)) are neighboring topologies, hence

    ℙ(τ^(j+1)|τ^(j)) = 1/|𝒩(τ^(j), q^(j))| =  
           1/|𝒩(τ^(j+1), q^(j))| = ℙ(τ^(j)|τ^(j+1)).

Therefore

    P_γ(s, s') = P_γ(σ(s'),  σ(s))

for any path γ.
This completes the proof.






We denote by C the set of pairs (s, s') ∈ A × B such that P(s, s')>0.
Let us consider any possible leap-prog step γ that connects s ∈ A and s' ∈ B crossing a finite number of boundaries and the sequences of augmented states (s=s^(0), s^(1), s^(2), …, s^(k) = s' ), topologies (τ=τ^(0), τ^(1), τ^(2), …, τ^(k) = τ' ), times (t=t^(0), t^(1), t^(2), …,t^(k) = t' ) and indices α = (α^(0), α^(1), …, α^(k)) (each α^(i) is a vector of  ± 1 entries characterizing the coordinates crossing zero in each sub-step) that decompose γ into segments on which the topology is unchanged.
By grouping the members of C by the value of α and ω, we have:

    C= ⋃_(α, ω)C_α, ω.

Because there will typically be many paths between s and s', the C_α, ω need not be disjoint.
However, we can modify the (countable number of) sets by picking one set for each (s, s') and dropping it from the rest, making a collection of disjoint sets { C_j } such that each C_j is a subset of some C_α, ω and

    C= ⋃_j ∈ JC_j.


We will write s ∈ A_j(s') and s' ∈ B_j(s) if (s, s') ∈ C_j and denote

    A_j=⋃_s' ∈ BA_j(s')and B_j=⋃_s ∈ AB_j(s).

We note that although the leap-prog algorithm is stochastic, if (α, ω) has been pre-specified, the whole path depends deterministically on the initial momentum.
Thus, by denoting the projection of C_α, ω to A by A_α, ω, we have that the transformation ϕ_α, ω that maps s to s' is well-defined on A_α, ω.
Since the projection of the particles (in a single leap-prog step) to the (q,p) space is exactly Reflective Hamiltonian Monte Carlo on ℝ^n_≥ 0.
Using Lemma 1, Lemma 2 and Theorem 1 in  <cit.>, we deduce that the determinant of the Jacobian of ϕ_α, ω is 1.

Now consider any j∈ J such that C_j ⊂ C_α, ω.
Because P(s,s') = P(s',s) for all s, s' ∈𝕋 and the determinant of the Jacobian of ϕ_α, ω is 1, we have

    ∫_B_j∑_s ∈ A_j(s')P(s',s) ds'   =  ∫_B_jP(s',ϕ_α, ω^-1(s')) ds'
       = ∫_A_jP(ϕ_α, ω(s),s)  ds
       = ∫_A_jP(s, ϕ_α, ω(s))  ds
       = ∫_A_j∑_s' ∈ B_j(s)P(s,s') ds.


Denote

    A^*=⋃_jA_jand B^*=⋃_jB_j.

Summing (<ref>) over all possible values of j gives

    ∫_B^*∑_s ∈ A(s')P(s',s) ds' = ∫_A^*∑_s' ∈ B(s)P(s,s') ds.

Moreover, we note that for s ∉ A^*, B(s)=∅.
Similarly, if  s' ∉ B^*, A(s')=∅.
We deduce that

    ∫_B∑_s ∈ A(s')P(s',s) ds' = ∫_A∑_s' ∈ B(s)P(s,s') ds.




By definition of k, for any state (τ', q') ∈ B, we can find a sequence of topologies (τ=τ^(0), τ^(1), τ^(2), …, τ^(k) = τ' ) for some l ≤ k such that τ^(i) and τ^(i+1) are adjacent topologies.
From the construction of the state space, let (τ^(i), q^(i)) denote a state on the boundary between the orthants for the two topologies τ^(i) and τ^(i+1).
Moreover, since (τ^(i), q^(i)) and (τ^(i+1), q^(i+1)) lie in the same orthant, we can find momentum values p^(i) and (p^(i))' such that

    P((τ^(i), q^(i), p^(i)) → (τ^(i+1), q^(i+1), (p^(i+1))')) >0

for all i.
That is, we can get from (τ^(i), q^(i), p^(i)) to (τ^(i+1), q^(i+1), (p^(i+1))') by a sequence of leapfrog steps Σ^(i) with length T.
By joining the Σ^(i)'s, we obtain a path Σ of k PPHMC steps that connects (τ^(0), q^(0)) and (τ',q').





For a path Σ of k PPHMC steps connecting (τ^(0), q^(0)) and (τ',q'), we define F_Σ = {(τ^(0), q^(0)), (τ^(1), q^(1)), …, (τ^(n_ω), q^(n_ω))}, where (τ^(i), q^(i)) denotes the state on Σ that joins the topologies τ^(i) and τ^(i+1).
We first note that although our leap-prog algorithm is stochastic, if the sequence of topologies crossed by a path Σ has been pre-specified, the whole path depends deterministically on the sequence of momenta p=(p^(0), …, p^(m)) along Σ.
Thus, the functions

    ϕ_i, ω(p):= q^(i)∀ p ∈ I_B, ω,

are well-defined.

We will prove that ϕ_n_ω, ω is Lipschitz by induction on n_ω. For the base case n_ω=0, the sequence ω is of length 1, which implies no topological changes along the path.
The leap-prog algorithm reduces to the baseline leapfrog algorithm and from standard results about HMC on Euclidean spaces <cit.>, we deduce that ϕ_1, ω is Lipschitz.

Now assume that the results holds true for n_ω = n.
Consider a sequence ω of length n+1.
For all (τ', q') ∈ B_ω, let Σ(τ', q') be a (k,T)-path connecting (τ^(0), q^(0)) and (τ', q').
We recall that

    F_Σ(τ', q')
       = {(τ^(0), q^(0)), (τ^(1), ϕ_1, ω(p)), …, (τ^(n_ω), ϕ_n_ω, ω(p))},

where ϕ_n_ω, ω(p))=(τ', q'), is the set of states that join the topologies on the path Σ(τ', q').

Define ω' = {τ^(0), τ^(1), …, τ^(n_ω-1)} and B'=ϕ_n_ω-1(I_B,ω), the induction hypothesis implies that the function ϕ_n_ω', ω' = ϕ_n_ω-1, ω is Lipschitz on I_B', ω'=I_B, ω.

On the other hand, since (τ^(n), q^(n)) and  (τ^(n+1), q^(n+1)) belong to the same topology, the base case implies that q^(n+1) is a Lipschitz function in p and q^(n) = ϕ_n_ω-1, ω(p).
Since compositions of Lipschitz functions are also Lipschitz, we deduce that ϕ_n_ω, ω is Lipschitz.

Since Lipschitz functions map zero measure sets to zero measure sets <cit.>, this implies μ(B_ω)=0 which completes the proof.






 §.§ Ergodicity of PPHMC





We denote

    ν(τ, q,p)    = 1/Zexp(-U(τ, q)) exp(-K(p)) 
       = 1/Zexp(-H(τ, q,p))

and refer to it as the canonical distribution.

It is straightforward to check that for all s, s' ∈𝕋, we have ν(s)r(s, s') = ν(s')r(s' , s).
Lemma <ref> implies that

    P(s, ds') ds = P(s', ds) ds'

in term of measures.
This gives the detailed balance condition

    ∫_A ∫_B ν(s)r(s, s') P(s, ds') ds
    
    =    ∫_B ∫_A ν(s')r(s' , s) P(s', ds) ds'

for all A, B ⊂𝕋.

We deduce that every update step of the second step of PPHMC satisfies detailed balance with respect to ν and hence, leaves ν invariant.
On the other hand, since ν is a function of |p|, the negation of the momentum p at the end of the second step also fixes ν.
Similarly, in the fist step, p is drawn from its correct conditional distribution given q and thus leaves ν invariant.

Since the target distribution π is the marginal distribution of ν on the position variables, PPHMC also leaves π invariant.




 §.§ Approximation error of reflective leapfrog algorithm


In this section, we investigate the local approximation error of the reflective leapfrog algorithm <cit.> without using surrogates.
Recall that V^+ and V^- are the restrictions of the potential function V on the sets {x_1 ≥ 0} and {x_1 ≤ 0},
and we assume that V^+ and V^- are smooth up to the boundary of their domains.


Consider a reflective leapfrog step with potential energy function V starting at (q^(0), p^(0)) (with q_1^(0)>0), ending at (q^(1), p^(1)) (with q_1^(1)<0) and hitting the boundary at x (with x_1=0, i.e., a refraction event happens on the hyper-plane of the first component).


Let p and p' denote the half-step momentum of a leapfrog step before and after the refraction events, respectively.
Recall that in a leapfrog approximation with refraction at x_1=0, we have

    p_i^(0) = p_i + ϵ/2∂ V/∂ q_i(q^(0)),  p_i^(1) = p_i' - ϵ/2∂ V/∂ q_i(q^(1)),

where

    p_1' = √(p_1^2 - 2dV(x)),

p_i' = p_i for i>1, and dV(x) = V^-(x) - V^+(x) denotes the change in potential energy across the hyper-plane.

The change in kinetic energy after this leapfrog step is

    Δ K =     -d V(x) - ϵ/2∑_i(p_i∂ V/∂ q_i(q^(0)) + p_i' ∂ V/∂ q_i(q^(1)) )
          +  ϵ^2/8∑_i((∂ V/∂ q_i(q^(1))  )^2- (∂ V/∂ q_i(q^(0)))^2 ).

We can bound the second-order term by

    (∂ V/∂ q_i(q^(1))  )^2- (∂ V/∂ q_i(q^(0)))^2 
            = 2 ∫_0^ϵ∂ V/∂ q_i(q^(0) + tp) ∂^2 V/∂ q_i^2(q^(0) + tp) p_i  dt
            = 𝒪(ϵ) ·sup_z, W=V^+, V^-∂ W/∂ q_i(z) ∂^2 W/∂ q_i^2(z).

On the other hand for the potential energy,

    Δ V    = V(q^(1)) - V(q^(0)) 
       = V(q^(1)) - V^-(x) + dV(x) +  V^+(x) - V(q^(0)) 
        = dV(x) +  ∫_ϵ_1^ϵ∇ V(q^(0) + tp) · p  dt
        +∫_0^ϵ_1∇ V(q^(0) + tp) · p  dt

where ϵ_1 and ϵ_2 := ϵ - ϵ_1  denote the integration times before and after refraction.
By the trapezoid rule for integration,

    Δ V = dV(x)    + ∑_i>1ϵ/2(p_i∂ V/∂ q_i(q^(0)) + p_i'∂ V/∂ q_i(q^(1))  ) 
       + p_1' ϵ_2/2∂ V/∂ q_1(q^(1)) +p_1' ϵ_2/2∂ V^-/∂ q_1(x) 
       + p_1ϵ_1/2∂ V/∂ q_1(q^(0)) +p_1 ϵ_1/2∂ V^+/∂ q_1(x)
       + 𝒪(ϵ^3) ·sup_z ∑_i, W=V^+, V^-(∂^3 W/∂ q_i^3(z) ) .

We recall that the error of the trapezoid rule on [a, b] with resolution h is a constant multiple of h^2(b-a), which is of order ϵ^3 in our case.
We deduce that

    Δ H= Δ V + Δ K =     -p_1' ϵ_1/2∂ V/∂ q_1(q^(1))  +p_1' ϵ_2/2∂ V^-/∂ q_1(x)
        -p_1 ϵ_2/2∂ V/∂ q_1(q^(0))  +p_1 ϵ_1/2∂ V^+/∂ q_1(x)  
       + 𝒪(ϵ^3).

Using Taylor expansion, we have

    ∂ V/∂ q_1(q^(1)) = ∂ V^-/∂ q_1(x) + 𝒪(ϵ),

and

    ∂ V/∂ q_1(q^(0))= ∂ V^-/∂ q_1(x)  +  𝒪(ϵ).

This implies

    Δ H =(ϵ_2 - ϵ_1) (p_1' ∂ V^-/∂ q_1(x) - p_1∂ V^+/∂ q_1(x) )  + 𝒪(ϵ^2).

In general, there is no dependency between ϵ_1 and ϵ_2, and the only cases where Δ H is not of order 𝒪(ϵ) are when

    √(p_1^2 - 2dV(x)) ∂ V^-/∂ q_1(x) - p_1∂ V^+/∂ q_1(x) =0.

In order for this to be true for all p, we need to have either

    dV(x)=0 and∂ V^-/∂ q_1(x) = ∂ V^+/∂ q_1(x),

or

    ∂ V^-/∂ q_1(x) = ∂ V^+/∂ q_1(x) = 0.

In both cases, the first derivative of V with respect to the first component must be continuous.
This completes the proof.





 §.§ Estimated posterior tree distributions for the simulated data




































































































































 §.§ Coordinate systems for branch lengths on trees


In this section we verify Assumption <ref> for phylogenetic trees.
Further explanation of the framework used here can be found in <cit.>.

Assume we are considering phylogenetic trees on N leaves, and that those leaves have labels [N] := {1, …, N}.
Every possible edge in such a phylogenetic tree can be described by its corresponding split, which is a partition of [N] into two non-empty sets, by removing that edge of the tree and observing the resulting partitioning of the leaf labels.
If a split can be obtained by deleting such an edge of a given phylogenetic tree, we say that the tree displays that split.
We use a vertical bar (|) to denote the division between the two sets of the bipartition.
For example, if we take the unrooted tree with four leaves such that 1 and 2 are sister to one another, the tree displays splits 1|234, 12|34, 134|2, 124|3, and 123|4.
Two splits A|B and C|D on the same leaf set are called compatible if one of A ∩ C, B ∩ C, A ∩ D, or B ∩ D is empty.
A set of splits that are pairwise compatible can be displayed on a phylogenetic tree <cit.>, and in fact the set of pairwise compatible sets of splits is in one-to-one correspondence with the set of (potentially multifurcating) unrooted phylogenetic trees.

When a single branch length goes to zero, 𝒩(τ, q) will have three elements: τ itself and its two NNI neighbors.
When multiple branch lengths go to zero, one can re-expand branch lengths for any set of splits that are compatible with each other and with the splits that did not originally go to zero.
This generalizes the NNI condition.
However, the correspondence between the branches that went to zero and the newly expanded branches is no longer obvious.

One can define such a correspondence using a global splits-based coordinate system.
Namely, such a coordinate system can be achieved by indexing branch length vectors by splits, with the proviso that for any two incompatible splits r and s, one of q_r or q_s is zero.
We could have used such a coordinate system for this paper, such that branch length vectors q would live in ℝ^2^N-1.

However, for simplicity of notation, we have indexed the branch lengths (e.g. in Algorithm <ref>) with integers [n] corresponding to the actual branches of a phylogenetic tree.
Thus our branch length vectors q live in 2N-3 dimensions.
One can use a total order on the splits to unambiguously define which branches map to which others when the HMC crosses a boundary.
We will describe how this works when two branch lengths, q_i and q_j, go to zero.
The extension to more branch lengths is clear.

Our branch indices i, j ∈ [2N-3] are always associated with a phylogenetic tree τ with numbered edges.
For any branch index i on τ, one can unambiguously take the split s_i.
Assume without loss of generality that s_i < s_j in the total order on splits.
Now, when q_i and q_j go to zero, one can transition to a new tree τ' which may differ from τ by up to two splits.
We assume without loss of generality that these are actually new splits (if not, we are in a previously defined setting) which we call s_1' and s_2' such that s_1' < s_2'.
We carry all of the branch indices for branches that aren't shrinking to zero across to τ'.
Then map branch i in τ to the branch in τ' corresponding to the split s_1', and branch j to the branch in τ' corresponding to the split s_2'.
Thus, for example, the momentum q_i in the τ orthant is carried over to this corresponding q_i in the τ' orthant.
