Parallel power posterior analyses for fast computation of marginal likelihoods in phylogenetics

Power posterior sampling

Both stepping-stone-sampling and path-sampling techniques construct and sample from a series of importance distributions. Lartillot & Philippe (2006) define the importance distributions as power posterior distributions, which are obtained by modifying the posterior probability density as (2)${\displaystyle f_{{\beta }_i}\left(\theta \right)=f{\left(Y|\theta ,M\right)}^{{\beta }_i}f\left(\theta |M\right)\text{.}}$

Here, β represent a vector of powers between 0 and 1. Then, for every value of β_i a draw from the power posterior distribution is needed and its likelihood score, l_i, is recorded (Lartillot & Philippe, 2006; Friel & Pettitt, 2008). In principle, one such likelihood sample per power posterior distribution is sufficient, although multiple samples improve the accuracy of the estimated marginal likelihood considerably (Baele et al., 2012b; Oaks et al., 2019; Fourment et al., 2020). We will use the notation l_ij to represent the jth likelihood sample from the ith power posterior distribution.

We illustrate the mean log-likelihood over different values of β in Fig. 1. Commonly, the values of the powers β are set to the ith quantile of a beta(0.3, 1.0) distribution (Xie et al., 2011; Baele et al., 2012a). The rationale is that more narrowly spaced intervals are needed for the range of β where the expected likelihood changes most rapidly, i.e., for β values close to 0 (Fig. 1).

Draws from the power posterior distribution are obtained by running a modified Markov chain Monte Carlo (MCMC Metropolis et al., 1953; Hastings, 1970) algorithm:

1. Let θ_j denote the current parameter values at the jth iteration, initialized at random at the start of the MCMC algorithm.

2. Propose a new values θ′ drawn from a proposal kernel with density q(θ′|θ_j).

3. The proposed state is accepted with probability (3)${\displaystyle \alpha =min\left(1,\frac{f{\left(D|{\theta }^{\prime }\right)}^{{\beta }_i}}{f{\left(D|{\theta }_j\right)}^{{\beta }_i}}\times \frac{f\left({\theta }^{\prime }\right)}{f\left({\theta }_j\right)}\times \frac{q\left({\theta }_j|{\theta }^{\prime }\right)}{q\left({\theta }^{\prime }|{\theta }_j\right)}\right)\text{.}}$

4. Set θ_j+1 = θ′ with probability α and to θ_j+1 = θ_j otherwise.

As can be seen from this brief description of the modified MCMC algorithm, only the likelihood values need to be raised to the power β_i. All remaining aspects of the MCMC algorithm stay the same as the standard implementations in Bayesian phylogenetics (Huelsenbeck & Ronquist, 2001; Drummond & Rambaut, 2007; Lakner et al., 2008; Lartillot, Lepage & Blanquart, 2009; Höhna & Drummond, 2012).

It is important to note that every MCMC simulation for each power β_j ∈ β necessarily includes its own burn-in period before the first sample can be taken. The power posterior analysis can be ordered to start from the full posterior (β_K−1 = 1.0) and then to use monotonically decreasing powers until the prior (β₀ = 0.0) has been reached. Thus, the last sample of the previous power posterior run can be used as the new starting state. This strategy has been shown to be more efficient because it is easier to disperse from the (concentrated) posterior distribution to the (vague) prior distribution thereby reducing the burn-in period significantly (Baele et al., 2012a).

Parallel power posterior analyses

The sequential algorithm of a power posterior analysis starts with a pre-burnin phase to converge to the posterior distribution. Then, consecutive power posterior simulations are performed sequentially, starting with β_K−1 = 1.0 (i.e., the posterior) to β₀ = 0.0 (i.e., the prior). Each power posterior simulation contains L iterations, with the likelihood of the current state recorded every Tth iteration. These ‘thinned’ samples are less correlated than the original draws from the MCMC simulation. The number of samples taken per power is n = L/T. At the beginning of each run a short burn-in phase is conducted, for example 10% or 25% of the run length.

The parallel algorithm for a power posterior analysis is set up almost identically to the sequential algorithm (see Fig. 2). Let us assume we have M CPUs available. Then, we split the set of powers into M consecutive blocks; the mth block containing the powers from ${\beta }_{\lfloor \ast K-\frac{\left(m-1\right)}{M}K-1\rfloor }$ to β_{⌊∗K−(mK/M)⌋}, e.g., the first out of four blocks for 128 analyses contains {β₁₂₇, …, β₉₆}, the second block contains {β₉₅, …, β₆₄}, etc. If the set of β cannot be split evenly into blocks then some blocks have one additional simulation, which is enforced by using only the integer part of the index. This block-strategy ensures that each CPU works on a set of consecutive powers which has the advantage of a shorter burn-in between simulations because the importance distributions are more similar to one another.

Regardless, each parallel sampler needs to start with an independent pre-burnin phase which creates an additional overhead. Thus, instead of running only one pre-burnin phase, as under the sequential power posterior analysis, we need to run M pre-burnin phases. This overhead could be removed only if it would be possible to draw initial values directly from the power posterior distribution.

Figure 2 shows a schematic of our parallelization algorithm. After the initial pre-burnin phase, the workload is divided into blocks and equally distributed over the available CPUs. Note that CPUs can be combined for distributed likelihood computation. No synchronization or communication between samplers is necessary because each power posterior simulation is independent. The only parallelization barrier occurs at the end when all power posterior simulations have finished. Finally, the master CPU collects all likelihood samples, combines the results, and computes the marginal likelihood using one of equations given below. These equations are computationally cheap compared with obtaining the likelihood samples. We thus expect that the performance gain is close to linear with the number of available cores. The algorithm described here is implemented in the open-source software RevBayes (Höhna et al., 2014; Höhna et al., 2016), available at http://www.RevBayes.com.

Our implementation in RevBayes uses the Message Passing Interface (MPI). That is, a RevBayes instance that was compiled using MPI can be used on any standard Unix based computer or high performance cluster (HPC) and executed in parallel. From the user perspective, no additional commands between using the standard version of RevBayes and using the MPI version are needed. For example, the command powerPosterior(mymodel, moves, monitors, “output/powers.out”, cats=100, sampleFreq=10) will automatically perform the power posterior analysis—as used in the study—in parallel. Hence, our implementation in RevBayes takes care of the parallelization for the user without the need of further specifications, assuming RevBayes was executed using MPI and several processes (e.g., mpirun -np 16 rb myscript.Rev). Such an analyses can be run on any standard HPC with submission systems such as SLURM and TORQUE.

Path-sampling

Path-sampling was the first numerical approximation method developed for marginal likelihood computation in Bayesian phylogenetic inference (Lartillot & Philippe, 2006). Path-sampling uses the trapezoidal rule to compute the integral of the log-likelihood samples between the prior and the posterior (see Fig. 1), which equals the marginal likelihood (Lartillot & Philippe, 2006). The equation of the trapezoidal rule for a single likelihood sample from each power posterior simulation is (4)${\displaystyle lnf\left(D|M\right)=\sum _{k=0}^{K-1}\frac{\left(ln\left(l_k\right)+ln\left(l_{k+1}\right)\right)\ast \left({\beta }_{k+1}-{\beta }_k\right)}{2}\text{.}}$ Samples of the log-likelihood have a large variance. Hence, it is more robust to take many log-likelihood samples and use the mean instead. This yields the equation to estimate the marginal log-likelihood, (5)${\displaystyle lnf\left(D|M\right)=\sum _{k=0}^{K-1}\frac{\left(\frac{\sum _{i=1}^{n}ln\left(l_{k,i}\right)}{n}+\frac{\sum _{i=1}^{n}ln\left(l_{k+1,i}\right)}{n}\right)\ast \left({\beta }_{k+1}-{\beta }_k\right)}{2}}$ which was proposed by Baele et al. (2012a).

Stepping-stone-sampling

Stepping-stone-sampling approximates the marginal likelihood by computing the ratio between the likelihood sampled from the posterior and the likelihood sampled from the prior. However, this ratio is unstable to compute and thus a series of intermediate ratios is computed: the stepping-stones (Xie et al., 2011; Fan et al., 2011). The stepping-stones can be chosen to be exactly the same powers as those used for path-sampling. The equation to approximate the marginal likelihood using stepping stone sampling is (6)${\displaystyle f\left(D|M\right)=\prod _{k=0}^{K-1}\left(\frac{1}{n}\sum _{i=1}^n\frac{l_{k,i}^{{\beta }_{k+1}}}{l_{k,i}^{{\beta }_k}}\right)=\prod _{k=0}^{K-1}\left(\frac{1}{n}\sum _{i=1}^n{l}_{k,i}^{{\beta }_{k+1}-{\beta }_k}\right)}$

Numerical stability of the computed marginal likelihood can be improved by retrieving first the highest log-likelihood sample, denoted by max_k, for the kth power. Re-arranging Eq. (6) accordingly yields (7)${\displaystyle ln\left(f\left(D|M\right)\right)=\sum _{k=0}^{K-1}\left[\right.ln\left(\sum _{i=1}^n\frac{exp\left(\right.\left(ln\left(l_{k,i}\right)-\underset{k}{max}\right)\ast \left({\beta }_{k+1}-{\beta }_k\right)\left)\right.}{n}\right)+\left({\beta }_{k+1}-{\beta }_k\right)\ast \underset{k}{max}\left]\right..}$

As seen in Eqs. (5) and (7), only the set of likelihood, or log-likelihood, samples is needed to approximate the marginal likelihood. Both marginal likelihood estimates approach the true marginal likelihood when the number of samples and powers increases. Since both computations are comparably fast, they can be applied jointly and, for example, be used to test for accuracy without additional time requirements.

Simulation design

Several previous studies have evaluated and compared the accuracy of different approaches to estimate the marginal likelihood of phylogenetic models (Baele et al., 2012a; Baele, Lemey & Suchard, 2016; Oaks et al., 2019). Overall, the findings suggest that power posterior based estimators (PS, SS and GSS) outperform other estimates for the cost of being computationally more expensive. The objective of this simulation study was to test the performance gain when using multiple CPUs. Thus, we tested the performance of the parallel power posterior analyses using two phylogenetic examples; a smaller and a larger dataset. As the small example dataset we chose 23 primate species representing the majority of primate genera. We used only a single gene sequence, the cytochrome b subunit, containing 1141 base pairs. For the larger example data set we chose an alignment with 4 genes (with a total of 6720 base pairs) from 305 taxa of the superfamily Muroidea (Schenk, Rowe & Steppan, 2013). For both examples we used the same model with the only difference that the larger dataset was partitioned into four subsets of sites (see protocols 1 and 2 from Höhna, Landis & Heath, 2017). We assumed that molecular evolution can be modeled by a general time reversible (GTR) substitution process (Tavaré, 1986) with four gamma-distributed rate categories (Yang, 1994). Furthermore, we assumed a strict, global clock (Zuckerkandl & Pauling, 1962) and calibrated the age of the root. As a prior distribution on the tree we used a constant-rate birth-death process with diversified taxon sampling (Höhna et al., 2011; Höhna, 2014) motivated by the fact that one representative species per genus was sampled, which is clearly a non-random sampling approach. The specific models correspond to the protocols described in Höhna, Landis & Heath (2017) and can also be found as tutorials at https://revbayes.github.io/tutorials/.

Each analysis consisted of a set of K = 100 power posterior simulations (see Fig. 2 for a schematic overview). The analyses started with a pre-burnin period of 10,000 iterations to converge to the posterior distribution. Note that in RevBayes each MCMC iterations consists of several moves (compared with MrBayes (Ronquist et al., 2012) and BEAST (Suchard et al., 2018) where each iteration consists of only a single move); the primates analysis included 38 moves and the Murdoidea analysis included 73 moves per iteration. Then, each power posterior analysis was run for 10,000 iterations and samples of the likelihood were taken every 10 iterations. The 25% initial samples of each power posterior distribution were discarded as additional burnin. The marginal likelihood was estimated using both path-sampling and stepping-stone-sampling once all power posterior simulations had finished as they contribute to performance overhead in practice. We ran each analysis 10 times and measured the computation time on the San Diego Supercomputer (SDSC) Gordon. Each compute node on Gordon contains two 8-core 2.6 GHz Intel EM64T Xeon E5 (Sandy Bridge) processors. The experiment was executed using 1, 2, 4, 8, 16, 32 and 64 CPUs, respectively. For each number k of CPUs used, we repeated the analyses by assigning 1, 2, 4,…64 CPUs to parallelizing the likelihood computation instead of distributing the stones. Thus, we additionally tested if parallelization over stones, the likelihood computation, or a mixture is most efficient.