Paper deep dive
FLARE MCMC: Fidelity-based Layer-Adaptive REcursive proposals for MCMC
Harini Venkatesan, Christian Shelton, Ming-Feng Ho, Simeon Bird, Mengxuan Wu
Intelligence
Status: succeeded | Model: Gemma-4-26B-A4B | Prompt: intel-v1 | Confidence: 93%
Last extracted: 8/17/2026, 4:20:11 AM
Summary
The paper introduces FLARE MCMC, a multi-fidelity layered Markov chain Monte Carlo method that accelerates mixing by recursively using lower-fidelity approximations of the likelihood function as proposals for higher-fidelity chains. This approach reduces computational cost while maintaining effective sample sizes, demonstrated in hydrology and cosmology applications where gradient-based methods are infeasible.
Entities (10)
Relation Signals (8)
Christian Shelton → affiliatedwith → University of California, Riverside
confidence 95% · Christian Shelton ... Department of Computer Science and Engineering, University of California, Riverside
Harini Venkatesan → affiliatedwith → University of California, Riverside
confidence 95% · Harini Venkatesan ... Department of Computer Science and Engineering, University of California, Riverside
FLARE MCMC → improves → mixing rate
confidence 95% · exploits lower-fidelity approximations of the true likelihood calculation to improve mixing and leads to overall faster performance
FLARE MCMC → uses → lower-fidelity approximations
confidence 95% · FLARE MCMC is a multi-fidelity layered MCMC method that exploits lower-fidelity approximations of the true likelihood calculation to improve mixing
FLARE MCMC → appliedto → Hydrology
confidence 90% · We demonstrate experimentally that FLARE MCMC achieves larger effective sample sizes for the same computational time across different scientific domains including hydrology
FLARE MCMC → appliedto → Cosmology
confidence 90% · We demonstrate experimentally that FLARE MCMC achieves larger effective sample sizes for the same computational time across different scientific domains including ... cosmology.
FLARE MCMC → fundedby → National Science Foundation
confidence 90% · This work was funded by the U.S. National Science Foundation (NSF) under Grant No. IIS-2435579.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Markov chain Monte Carlo (MCMC) requires only the ability to evaluate the likelihood, making it a common technique for inference in complex models. However, it can have a slow mixing rate, requiring the generation of many samples to obtain good estimates and an overall high computational cost. FLARE MCMC is a multi-fidelity layered MCMC method that exploits lower-fidelity approximations of the true likelihood calculation to improve mixing and leads to overall faster performance. Such lower-fidelity likelihoods are commonly available in scientific and engineering applications where the model involves a simulation whose resolution or accuracy can be tuned. Our technique uses recursive, layered chains with simple layer tuning; it does not require the likelihood to take any form or have any particular internal mathematical structure. We demonstrate experimentally that FLARE MCMC achieves larger effective sample sizes for the same computational time across different scientific domains including hydrology and cosmology.
Tags
Links
- Source: https://arxiv.org/abs/2608.13774v1
- Canonical: https://arxiv.org/abs/2608.13774v1
Trouble viewing inline? Open PDF directly →
Full Text
98,985 characters extracted from source content.
Expand or collapse full text
FLARE MCMC: Fidelity-based Layer-Adaptive REcursive proposals for MCMC Thanks: This is the author’s accepted manuscript of an article published in SIAM Journal on Uncertainty Quantification. The final version is available at https://doi.org/10.1137/25M1795194. Harini Venkatesan Email: hvenk001@ucr.edu Email: cshelton@cs.ucr.edu Email: mwu171@ucr.edu Thanks: Department of Computer Science and Engineering, University of California, Riverside (, , ) Christian Shelton11footnotemark: 1 Ming-Feng Ho Email: mfho@umich.edu Thanks: Leinweber Center for Theoretical Physics, University of Michigan () Simeon Bird Email: sbird@ucr.edu Thanks: Department of Physics and Astronomy, University of California, Riverside () Mengxuan Wu11footnotemark: 1 Abstract Markov chain Monte Carlo (MCMC) requires only the ability to evaluate the likelihood, making it a common technique for inference in complex models. However, it can have a slow mixing rate, requiring the generation of many samples to obtain good estimates and an overall high computational cost. FLARE MCMC is a multi-fidelity layered MCMC method that exploits lower-fidelity approximations of the true likelihood calculation to improve mixing and leads to overall faster performance. Such lower-fidelity likelihoods are commonly available in scientific and engineering applications where the model involves a simulation whose resolution or accuracy can be tuned. Our technique uses recursive, layered chains with simple layer tuning; it does not require the likelihood to take any form or have any particular internal mathematical structure. We demonstrate experimentally that FLARE MCMC achieves larger effective sample sizes for the same computational time across different scientific domains including hydrology and cosmology. keywords Markov chain Monte Carlo, multi-layered models, Bayesian inference, simulation-based inference, hydrology, cosmology. Funding.This work was funded by the U.S. National Science Foundation (NSF) under Grant No. IIS-2435579. SB was supported by NASA-80NSSC21K1840. †runningheads: FLARE MCMC: Fidelity-based Layer-Adaptive REcursive proposals for MCMC / H. Venkatesan, C. Shelton, M. Ho, S. Bird, and M. Wu MSC 62F15, 62M05, 65C05, 65C40, 85A35 1 Introduction Markov chain Monte Carlo (MCMC) is a workhorse of scientific and engineering computation. Most frequently, it is employed to compute the posterior distribution of model parameters, based on observations. The calculated distributions (as represented by samples) give estimates that can be used in calibration and uncertainty quantification to aid in the generation of new scientific experiments, clarify the observability of the model, and resolve scientific theories. Among the many MCMC algorithms, Metropolis-Hastings MCMC (MH-MCMC) is popular because of its ability to sample from almost any distribution while requiring only the ability to evaluate the model’s likelihood given a parameter setting. Yet, this is also its weakness, as it has no additional knowledge of the problem setting to guide its sampling effectively. Therefore, its mixing time (speed of generating effectively new samples) can be slow and the overall algorithm computationally burdensome. Methods such as Hamiltonian Monte Carlo and its variants [17, 47, 36] speed up mixing by adding auxiliary momentum variables, allowing longer steps to reduce correlations between consecutive samples. Such methods require computing the gradient of the log target distribution with respect to the parameters, something that could be prohibitively expensive when the distribution is evaluated through lengthy simulation code. For instance, the cosmological simulation we use in our experimental results that aims to approximate the posterior density conditioned on the galaxy power spectrum from SDSS-I Baryon Oscillation Spectroscopic Survey (BOSS) Data [13, 4] cannot be modified to produce gradients. Due to the complexity and non-differentiability of the forward cosmological simulation, gradients with respect to the model parameters are not available. Therefore, methods like auto-differentiation cannot be applied, nor is there an analytic form for the gradients, prohibiting the use of gradient-based inference methods. FLARE MCMC speeds up the mixing time of MH-MCMC by exploiting lower-fidelity models of the same problem. Many engineering or scientific computational models can be run at multiple fidelities. FLARE MCMC exploits a set of computationally cheaper posterior calculations, each an approximation of the true posterior. Many posteriors involve solving a PDE, ODE, or integral. For these, coarsening the spatial or temporal grid leads to cheaper approximations. For those with constraint or optimization solvers, reducing the solvers’ tolerances or maximum number of iterations can similarly lead to cheaper approximations. We further show a physics example where the underlying simulation can be coarsened by reducing the number of representative particles. By recursively employing MCMC chains, we can use the coarser resolution models to guide the higher resolution MCMC chain. The result is a sampler for the target model that converges faster and generates more effective samples per computation time, even considering the extra time necessary to employ the lower-fidelity computations. 2 Background Markov chain Monte Carlo is a class of algorithms designed to sample from a complicated target distribution by constructing an easy-to-simulate Markov chain such that the stationary distribution of the Markov chain is the target distribution. Commonly, this target distribution is the posterior distribution of a set of parameters, conditioned on observations. Let D be the observations and θ∈Θ⊂ℝRθ∈ ^R be the parameters. Assuming a prior distribution on the parameters p(θ)p(θ), the target posterior distribution of interest, π(θ∣D)π(θ D), is obtained through Bayes’ theorem: π(θ∣D)=ℒ(D∣θ)p(θ)p(D)∝ℒ(D∣θ)p(θ)π(θ D)= L(D θ)p(θ)p(D) (D θ)p(θ) (1) where ℒ(D∣θ)L(D θ) is the likelihood of the data, which in many scientific applications requires a lengthy simulation to evaluate. We only require the ability to evaluate π(θ∣D)π(θ D) up to a constant of proportionality, and therefore the denominator of p(D)p(D) is safely ignored. That π(θ∣D)π(θ D) is a conditional distribution is largely irrelevant for MCMC, so we will just let π(θ)π(θ) denote the distribution of interest (equal to π(θ∣D)π(θ D) if the underlying distribution is a posterior, but it could be any distribution over θ). 2.1 Metropolis-Hastings MCMC We begin by focusing on the Metropolis-Hastings method for Markov Chain Monte Carlo (MH-MCMC) introduced by Hastings (1970) [30]. The (i+1)(i\!\!+\!\!1)th sample, θi+1θ^i+1, is generated based on the previous sample in the chain, θiθ^i, in a two-step process. First, a proposed next state, θ~i θ^i is generated from a proposal distribution, q(θ~i|θi)q( θ^i|θ^i). Then, θ~i θ^i is either accepted or rejected as θi+1θ^i+1 according to a carefully constructed acceptance probability. If accepted, θi+1=θ~iθ^i+1\!=\! θ^i, otherwise θi+1=θiθ^i+1\!=\!θ^i. Often, a normal distribution centered at θiθ^i is used as the proposal distribution q(θ~i|θi)q( θ^i|θ^i), but almost any proposal distribution can be used, subject to mild conditions (for instance, that q(θ~i|θi)q( θ^i|θ^i) is positive everywhere). With a chosen q(θ~i|θi)q( θ^i|θ^i), the acceptance probability, A, for the transition θi→θ~iθ^i\! \! θ^i is (θi→θ~i)=min(1,r(θi→θ~i))A(θ^i\! \! θ^i)= (1,r(θ^i\! \! θ^i)) (2) where r(θi→θ~i)=π(θ~i)π(θi)q(θi|θ~i)q(θ~i|θi).r(θ^i\! \! θ^i)= π( θ^i)π(θ^i) q(θ^i| θ^i)q( θ^i|θ^i)\,\,. (3) Although the standard Metropolis-Hastings MCMC algorithm can be an easy way to sample from a posterior distribution, it requires sufficient samples to be an effective approximation of the posterior distribution. When the chain is slow to mix (due to a less-than-optimal proposal distribution), consecutive samples are highly dependent and more samples must be taken to achieve a set representative of the true distribution. When the evaluation of π(θi)π(θ^i) (necessary for the calculation of Equation 3) is computationally expensive, this is particularly problematic. 2.2 Related Work Our goal of accelerating MCMC sampling is shared by a large body of work. These approaches involve methods that couple chains (like simulated tempering), methods that aim to reduce the variance of estimators for a target using cheap approximations from multiple fidelities (like MLMC), and methods that use cheap models to build MCMC proposals. Like FLARE MCMC, methods such as simulated tempering and coupled MCMC [57, 46, 5] use multiple chains. Samples are accepted or rejected by evaluating the energy of the process and adjusting the temperature of the model. Two chains are run in parallel at different temperatures, and the system swaps between different temperatures. Reversible jump MCMC [24, 3] also jumps between chains (of different dimensions). While FLARE MCMC shares the notion of multiple chains, because it solves a different problem (to take advantage of simulations that are orders of magnitude cheaper to evaluate), the resulting structure is very different. Methods such as sequential MCMC or particle filtering [43, 16] use the notion of approximations of the target by a large number of samples called particles that are propagated across time using importance sampling. However, those are filtering frameworks and do not converge to a stationary distribution. Thus, though appearing related in its structure, FLARE MCMC is quite different to these methods. A highly influential body of work focuses on reducing the variance of the final estimator for a target expectation in the Multilevel Monte Carlo (MLMC) framework. Taking inspiration from the multilevel Monte Carlo method [31] for high-dimensional, parameter-dependent integrals and Multilevel Monte Carlo Path Simulation [22], Hoang et al. [34] proposed a multilevel MCMC method that applied to Bayesian Inverse problems. The core idea is to decompose a high-fidelity expectation into a telescoping sum using a hierarchy of computational models with increasing model resolution. This method achieves computational speedup by estimating the low variance difference terms with small number of samples, while the bulk of the computational efforts is spent on the cheap, low-fidelity estimator. The authors also provide rigorous complexity proof, showing how quickly posterior expectation might converge when running iterative samplers on sparse grids using telescopic expansion on the discretization error. This MLMC framework has been extended in many directions. Multilevel sequential Monte Carlo samplers [7] and Multilevel Particle Filters [39] along with previous work [35, 25, 26] extended MLHC to sequential Monte Carlo. Jasra et al. [41] extended MLHC to the problem of static parameter estimation in partially observed diffusions. Problems with multiple ways of discretizing were addressed by Multi-Index MCMC [29, 40]. These MLMC methods use samples from all chains in a telescoping estimator. They target the MSE of a specific quantity of interest. In contrast, FLARE MCMC uses only samples from the finest chain (like methods discussed below) and targets the chain’s mixing time, rather than MSE. Our theoretical analyses in this paper focus on the ergodicity and convergence rates for FLARE MCMC and are not specific to any particular problem domain. More similar to FLARE MCMC, several previous methods have shown that replacing the proposal with an approximation with generally high acceptance probability reduces the computational cost of the standard Metropolis-Hastings algorithm significantly. This idea was first proposed by Christen and Fox [9, 20] as a two-stage MCMC method that tests the original proposal using a cheap approximation to find moves in the chain that are more likely to be accepted. In other words, a candidate is accepted with the likelihood of the approximate model before it is evaluated with the more expensive model. In preconditioned MCMC using coarse-scale simulation proposed by Efendiev, Hou and Luo [18], two-stages are used to reduce the computational cost incurred in the fine fidelity by testing the coarse model based on high-fidelity multiscale finite volume model. However, this only performs a single check with a cheap approximation, and does not exploit it to run a full MCMC subchain. Multilevel Markov chain Monte Carlo (MLMCMC) [14] achieves computationl efficiency on the finer levels. If the coarse proposal from the approximation is rejected by the fine level, the coarse chain continues independently of the fine chain instead of recursively starting the next coarse chain from the current sample of the fine chain. MLMCMC uses a user-specified variable that is internal to the likelihood computation and shared across the levels (for instance, the predicted observations to be compared with the true observations through a noise model). The samples drawn from the coarse approximation are used to reduce the variance of this internal variable achieving better proposals from the coarse fidelities. Lykkegaard et al. [44] proposed Adaptive Multilevel Delayed Acceptance (MLDA), which adapted a recursive version of MLMCMC over multiple levels. Here, the coarse inner subchain used to generate subsequent proposals for the current chain is initiated from the current sample from the outer chain again instead of independently continuing the fine chain even if coarse proposal is rejected. MLDA also applies an Adaptive Error Model (AEM) [42] to account for discrepancies between the different fidelities. It takes the two-level AEM from Adaptive Delayed Acceptance Metropolis Hastings [12, 11] and extends it by adding a telescoping sum of differences in the model output across multiple levels. Several multilevel MCMC methods based on delayed rejection, in contrast to delayed acceptance, have also been proposed and are summarized by Peherstorfer et al. [49]. Adaptive methods in multistage MCMC [58] proposed using an independence sampler that is a good approximation for the posterior distribution in the first stage and random walk in the second stage to help with poor approximation by the independence sampler. Delayed rejection in MCMC [23] suggested using a normal distribution as the proposal in the first level and a normal distribution with the same mean but higher variance in the second level. Higdon et al. [32] proposed using multiple MCMC chains from low and high fidelities and coupling them using a product chain and ”swapping” updates allowing information to move between the two fidelity scales. An accelerated MCMC method using local approximations was developed by Conrad et al. [10] that uses local approximations of either the log-likelihood function or the forward model of different simulations into the Metropolis-Hastings kernel. Although these methods use approximations as proposals, they do not exploit layered or recursive MCMC chains. Cai and Adams [8] proposed a multi-fidelity Monte Carlo method (MFMC) that uses randomized fidelities as the approximate for the target fidelity. The algorithm does not converge to the true posterior, but the resulting samples can be used to estimate expectations through a specific “sign-correction” formula. Our method follows a hierarchy of levels in its sampling while also sampling from the true posterior and provides a simpler alternative to previous multilevel methods. Our multi-fidelity layered MCMC algorithm, FLARE MCMC, has a similar structure to MLDA in terms of the recursive layers and achieves a similar amount of effective samples across multiple chains of MLDA. However, our method for mitigating the differences between approximations is simpler in construction and implementation than that of MLDA, does not require the identification of any internal variables of the distribution to be sampled, and generates more effective samples in a shorter amount of time and computational cost. We demonstrate this on real-world large scientific problems. We also show theoretical convergence rates, optimal value for the number of inner steps M and prove ergodicity of the adaptation in layer tuning. 3 FLARE MCMC: Fidelity-based Layer-Adaptive REcursive proposals for MCMC We consider a series of models, ordered by fidelity. For instance, we might have a model that evaluates a differential equation numerically as the main part of the likelihood calculation (simulating forward in time); the resolution of the spatial or temporal grid used to evaluate the model can be tuned to change its fidelity. The highest fidelity model is our “true” model, from whose posterior we wish to sample. FLARE MCMC draws samples from the true model. Its nested chains use the coarser fidelity models as cheap approximations of this finest fidelity model to speed up mixing. Where the standard Metropolis-Hastings algorithm uses a distribution q that proposes the next sample, FLARE MCMC uses nested Markov chains as the proposal distribution. In a recursive fashion, each layer uses the result of another MCMC chain with a coarser approximation as its proposal. The recent sample in the current chain is the starting sample in the nested chain. The coarser chain runs for M iterations, with each proposed sample evaluated by the likelihood of the cheaper layer. The last, MMth, sample of the coarser chain is proposed as the candidate for the next sample in the current chain. At the coarsest fidelity/layer, a standard proposal distribution is used, for instance a normal distribution centered on the current point. This avoids numerous expensive likelihood calculations in the fine fidelity that might end up rejected, and it allows the proposal to generate samples that are more likely to get accepted by the finest fidelity, since it was accepted by an approximation already. While the coarser chains have their own computational cost, they can often be orders of magnitude faster to evaluate, thus leading to an overall savings in the running time of the entire algorithm, as measured by the quality of the samples generated per computational time. Figure 1: One sampling step from the finest layer with two coarse fidelities and two iterations per nested chain. Refer to text in Section 3.1. 3.1 Algorithm Specification Let θ∈ℝRθ ^R be the set of parameters (over which we are sampling) and let j∈0,1,…,Jj∈\0,1,...,J\ be the fidelities ordered in a decreasingly complex fashion (00 is the “true” model and J is the coarsest fidelity). We let πj(θ) _j(θ) be the posterior distribution according to the jjth fidelity model and qj(θ~|θ)q_j( θ|θ) be the proposal distribution for layer j. Here, θji _j^i is the iith sample in the current chain at layer j. The goal is to sample from π0(θ) _0(θ). Algorithm 1 FLARE-Chain(θj0θ^0_j,n,j) for i=0,…,n−1i=0,…,n-1 do if j=Jj=J then ⊳ coarsest layer Sample θ~ji θ^i_j from qj(⋅|θji)q_j(·|θ^i_j) Accept θji+1=θ~jiθ^i+1_j= θ^i_j with probability A from Equation 2 Otherwise, reject and θji+1=θjiθ^i+1_j=θ^i_j else θj+11,…,θj+1M=θ^1_j+1,…,θ^M_j+1= FLARE-Chain(θjiθ^i_j, M, j+1j\!+\!1) θ~ji=θj+1M θ^i_j=θ^M_j+1 Accept θji+1=θ~jiθ^i+1_j= θ^i_j with probability jA_j from Equation 5 Otherwise, reject and θji+1=θjiθ^i+1_j=θ^i_j return θj1,…,θjnθ^1_j,…,θ^n_j In FLARE MCMC, the proposal distribution qj(θ~ji|θji)q_j( θ_j^i| _j^i) for iteration i of a chain at layer j is another MCMC chain of M steps targeting the (coarser) posterior πj+1(⋅) _j+1(·), starting this nested chain at θji _j^i. The result of M steps using a chain with stationary distribution πj+1(⋅) _j+1(·) is the proposal for θ~ji θ_j^i: qj(θ~ji|θji)q_j( θ_j^i| _j^i). More algorithmically, to generate θ~ji θ^i_j from θjiθ^i_j, we run the (coarser) MCMC algorithm at layer j+1j+1. We start with θj+10=θjiθ^0_j+1=θ^i_j and continue the coarser MCMC sampler until θj+1Mθ^M_j+1. We then set θ~ji=θj+1M θ^i_j=θ^M_j+1. At the coarsest layer, qJ(θ~i|θi)q_J( θ^i|θ^i) is a standard simple proposal distribution. Figure 1 pictorially demonstrates this for J=2J=2 inner layers, each with M=2M=2 steps. With the sampling scheme so defined, it remains to construct the acceptance probability for each layer: 0,1,…,JA_0,A_1,…,A_J. We follow a standard Metropolis-Hastings method for every layer and therefore j=min(1,rj(θji→θ~ji))A_j= (1,r_j( _j^i\! \! θ_j^i)). At the coarsest layer, the ratio rJ(θJi→θ~Ji)r_J( _J^i\! \! θ_J^i) is just as in Equation 3 because qJq_J is a standard proposal distribution. When j<Jj<J, the proposal distribution is from a Markov chain that obeys detailed balance. Therefore qj+1(θji|θ~ji)qj+1(θ~ji|θji)=πj+1(θji)πj+1(θ~ji)0≤j<J q_j+1( _j^i| θ_j^i)q_j+1( θ_j^i| _j^i)= _j+1( _j^i) _j+1( θ_j^i) 0≤ j<J (4) and thus j(θji→θ~ji)=min(1,πj(θ~ji)πj(θji)⋅πj+1(θji)πj+1(θ~ji)).A_j(θ^i_j\! \! θ^i_j)= (1, _j( θ^i_j) _j(θ^i_j)· _j+1(θ^i_j) _j+1( θ^i_j) )\,\,. (5) Note this equation does not depend on M (the number of steps for the coarser chain at layer j+1j+1). While this chain has almost certainly not mixed for small M, the ratio qj(θji|θ~ji)/qj(θ~ji|θji) q_j( _j^i| θ_j^i)q_j( θ_j^i| _j^i) is the same as if the chain had completely mixed and the proposed new state, θ~ji θ_j^i, were from the true posterior of the model at layer j+1j+1. The values πj+1(θ~ji) _j+1( θ^i_j) and πj+1(θji) _j+1(θ^i_j) were already calculated as part of the chain at layer j+1j+1 and therefore do not take any additional computation time. FLARE MCMC is summarized in Algorithm 1. To gather N samples from the true posterior, the algorithm is called with FLARE-Chain(θ0,N,J=0)(θ^0,N,J=0). 3.2 Convergence Rate We show a convergence rate for FLARE MCMC. We measure the distance to the stationary in terms of total variation distance as follows. Definition 1 (Strasser (1985) [55]). The total variation distance between two probability measures ν1 _1 and ν2 _2 is defined as ‖ν1−ν2‖=supA|ν1(A)−ν2(A)|.\| _1- _2\|= _A| _1(A)- _2(A)|. The minorization condition of Markov chains, used by Roberts and Rosenthal [51], provides a means of bounding the convergence rate. For any Markov chain with a one step transition probability of p(θi→θi+1)p(θ^i→θ^i+1), we let pn(θi→θi+n)p^n(θ^i→θ^i+n) denote the corresponding n step transition probability. Formally, the definition of the minorization condition is stated below. Definition 2 (Roberts and Rosenthal (2004) [51]). A Markov chain on Θ satisfies the minorization condition if there exists an ϵ>0ε>0, a positive integer n, and a probability measure ν(.)ν(.) such that pn(θ0→θn)≥ϵν(θn)∀θ0,θn∈Θ.p^n(θ^0→θ^n)≥εν(θ^n) ∀θ^0,θ^n∈ . (6) With respect to the stationary distribution, the probability of transitioning from θ to θ′θ can be minorized by a lower bound such that pj1(θ→θ′)≥ϵπj(θ′)p_j^1(θ→θ )≥ε _j(θ ). Lemma 3. Assume the minorization condition holds at the innermost level (j=Jj=J): pJ1(θJi→θJi+1)≥ξJ⋅πJ(θJi+1)p_ J^1(θ^i_ J→θ^i+1_ J)≥ _J· _ J(θ^i+1_J) for some ξJ>0 _ J>0. Then, there exists a minorized lower bound on levels j<Jj<J such that pj1(θji→θji+1)≥ξj⋅πj(θji+1)∀θji,θji+1p_j^1(θ^i_j→θ^i+1_j)≥ _j· _j(θ^i+1_j) ∀θ^i_j,θ^i+1_j (7) where ξj=(1−(1−ξj+1)M)⋅minθ(πj+1(θ)πj(θ)) _j=(1-(1- _j+1)^M)· _θ ( _j+1(θ) _j(θ) ). We call ξj _j the minorization constant for level j. This satisfies the necessary minorization condition of Theorem 8 in the original paper [51]. This allows us to get a quantitative bound on the distance to the stationary of every level as stated in Theorem 4. The proof for Lemma 3 can be found in Appendix A. Theorem 4. Let pjn(θj0→⋅)p_j^n(θ^0_j→·) be the distribution for layer j with an invariant target probability πj(⋅) _j(·). FLARE MCMC is uniformly ergodic and converges as ‖pjn(θj0→⋅)−πj(⋅)‖≤(1−ξj)n \|p_j^n(θ^0_j→·)- _j(·) \|≤(1- _j)^n where ξj=(1−(1−ξj+1)M)⋅minθ(πj+1(θ)πj(θ)) _j=(1-(1- _j+1)^M)· _θ ( _j+1(θ) _j(θ) ) as in Lemma 3. Most critically, it holds for layer j=0j=0. Proof. The results follow from the minorization condition established in Lemma 3. Although the theorem above establishes the convergence of the chain to the invariant target distribution, it illustrates several aspects about the inner chains. In particular, ergodicity requires that the chain be able to reach all regions of the target’s support. Therefore, the coarser approximations’ supports must be supersets of the finer ones. The coupling strength, ξj _j, in turn depends on the coupling strength of the coarser approximations’ chains, ξk,k>j _k,k>j. Thus, the effects of the mixing times of the inner chains on the outer chain are captured in this theorem. The similarities of the approximations to each other are captured in the minθ(πj+1(θ)πj(θ)) _θ ( _j+1(θ) _j(θ) ) terms. Thus, the theorem also quantifies the effects of similarities and dissimilarities between the approximations on the total convergence rate. If the modes of the target distribution are preserved across coarsening, then we would expect these terms to be larger and therefore the outer chain to mix faster. Ideally, the minimum in this term could be replaced with an expectation, thus turning it into the KL-divergence between adjacent layers. We have not yet determined whether or how this might be possible. 3.3 Optimal number of inner steps M To understand how to select the number of inner steps for each layer MjM_j, we derive a theoretical expression for optimal MjM_j that balances the decoupling rate of the chains and the cost of likelihood evaluations. The following lemma provides an analytical expression for the value of MjM_j that maximizes the cost-aware decoupling rate. The proof for Lemma 5 can be found in Appendix B. Lemma 5. For layers 0≤j<J0≤ j<J, suppose the minorization condition hold such that pj1(θji→θji+1)≥ξj⋅πj(θji+1)p_j^1(θ^i_j→θ^i+1_j)≥ _j· _j(θ^i+1_j), where ξj=(1−(1−ξj+1)Mj)⋅minθ(πj+1(θ)πj(θ)) _j= (1-(1- _j+1)^M_j )· _θ ( _j+1(θ) _j(θ) ) is the minorization constant, and Mj≥0M_j≥ 0 is the number of inner steps. Let the total cost per step at level j be Bj=bj+Mj⋅Bj+1B_j=b_j+M_j· B_j+1 where bj>0b_j>0 is the cost of evaluating the likelihood at level j and Bj+1>0B_j+1>0 is the cost of a single evaluation of the inner layer j+1j+1. Define the cost-aware computational decoupling rate as f(Mj)=(1−(1−ξj+1)Mj)⋅minθ(πj+1(θ)πj(θ))Bj.f(M_j)= (1-(1- _j+1)^M_j )· _θ ( _j+1(θ) _j(θ) )B_j. (8) Then the real-valued maximizer Mj∗M_j^* of f(Mj)f(M_j) is Mj∗=−1ΥW−1(−e−Υμ)−μM_j^*=- 1 W_-1 (-e^- μ )-μ (9) where Υ=−log(1−ξj+1) =- (1- _j+1), μ=bjBj+1+1Υμ= b_jB_j+1+ 1 and W−1W_-1 is −1-1 branch of the Lambert W function. However, this expression is not directly usable in practice, since it depends on the unknown coupling minorization constant ξj+1 _j+1 which is generally unknown and difficult to approximate in MCMC settings. The optimizer Mj∗M^*_j is a real-valued quantity, whereas in practice the number of inner steps must be an integer. Instead, we have found empirical values for M that offers the best trade-off between computation time and sampling efficiency in the experiment section. Nevertheless, this lemma provides a theoretical benchmark for the optimal trade-off between computational cost of likelihoods and effective mixing across layers. −2 -2−1 -10 01 12 2initial angle, α0 _0 (radians)0.5 0.51.0 1.01.5 1.52.0 2.0pendulum length, L LFine Model−2 -2−1 -10 01 12 2initial angle, α0 _0 (radians)0.5 0.51.0 1.01.5 1.52.0 2.0pendulum length, L LCoarse Model Figure 2: Posterior of different fidelities; coarse model is the small angle approximation. See Section 4.1. 3.4 Layer Tuning The algorithm above uses the coarser fidelities to guide the finer ones. Early in the chain, this is useful for quickly driving the samples toward high-probability regions. However, this mismatch between the fidelities can cause problems later because it can steer the chain away from high-probability regions in the fine fidelity model that do not overlap with high-probability regions of the coarse fidelity model. Figure 2 demonstrates an example of such partial, but not complete, overlap in one of our examples. To combat this, we present a simple modification that does not require estimation of any internal variables of the probability models, nor estimation of means or variances from multiple chains. Recall πj(θj) _j( _j) is known up to a normalizing constant: πj(θj)=π~j(θj)/Zj _j( _j)= π_j( _j)Z_j where π~j(θj) π_j( _j) is the unnormalized distribution and ZjZ_j is the normalzing constant. We modify the target distributions for coarser chains (and thus the proposal distributions for all j>0j>0) as ψj(θj) _j( _j) =(π~j(θj)+ωj)/ζj(ωj) = ( π_j( _j)+ _j ) _j( _j) ∀<j≤J ∀\ \ 0\!<\!j\!≤\!J (10) where ζj(ωj) _j( _j) is the normalizing constant of this new distribution which depends on ωj _j.11 1 We assume the domain of θ, Θ , is of finite volume. We now use ψj(θj) _j( _j) in place of πj(θj) _j( _j) in Equation 5, therefore modifying the acceptance ratio for all <j≤J\ 0\!<\!j\!≤\!J\,\, as j(θji→θ~ji) _j(θ^i_j\! \! θ^i_j) =min(1,ψj(θ~ji)ψj(θji)⋅ψj+1(θj+1i)ψj+1(θ~j+1i)). = (1, _j( θ^i_j) _j(θ^i_j)· _j+1(θ^i_j+1) _j+1( θ^i_j+1) )\,\,. (11) For the finest layer, things remain the same (or alternatively, ω0=0 _0=0), because we do not want to change the distribution of the overall sampler. This effectively mixes the stationary distribution of the jjth layer with a uniform distribution (we have added a constant to the posterior and then renormalized), encouraging the proposal to explore more widely than the coarser layer would normally. While unsophisticated, we found it simpler to implement and compute than other options and just as effective. Instead of leaving (ω1,ω2,…,ωJ)( _1, _2,…, _J) as hyper-parameters, we use gradient descent to adapt them over the course of the sampling. We adjust ωj+1 _j+1 to minimize the Kullback-Leibler divergence between layers ψj _j and ψj+1 _j+1: KL(ψj∥ψj+1)=θ∼ψj[ln(ψj)]−θ∼ψj[ln(ψj+1)].KL( _j\| _j+1)= _θ _j[ ( _j)]- _θ _j[ ( _j+1)]\,\,. (12) This tries to make the coarser (approximating) distribution ψj+1 _j+1 more similar to the distribution ψj _j. Because the first term does not depend on ωj+1 _j+1, the objective function is to maximize Hj+1=θ∼ψj[ln(ψj+1)].H_j+1= _θ _j [ ( _j+1) ]\,\,. (13) Using Equation 10, ∂ωj+1Hj+1 ∂ _j+1H_j+1 =∂ωj+1(θ∼ψj[ln(π~j+1(θ)+ωj+1)]−lnζj+1(ωj+1)) = ∂ _j+1 ( _θ _j [ ( π_j+1(θ)+ _j+1 ) ]- _j+1( _j+1) ) =θ∼ψj[∂ωj+1ln(π~j+1(θ)+ωj+1)]−θ∼ψj+1[∂ωj+1ln(π~j+1(θ)+ωj+1)] = _θ _j [ ∂ _j+1 ( π_j+1(θ)+ _j+1 ) ]- _θ _j+1 [ ∂ _j+1 ( π_j+1(θ)+ _j+1 ) ] (14) where the second step replaces the derivative of the log-partition function with the expected derivative of the log-probability. The first term is an expectation with respect to the distribution at the lower layer j. We assume that the lower layer has mixed and therefore, the starting state for the chain at layer j+1j+1 is a sample drawn from ψj _j. The second term is an expectation with respect to the distribution at this layer, j+1j+1. We let the sample at the end of this chain after M steps approximate a sample from this distribution. This is similar to the approximation employed by M-step contrastive divergence [33]. Although this is not guaranteed to converge [56], in practice we have found it to work well. Thus, the total derivative for the gradient ascent update is ∂ωj+1Hj+1≈1π~j+1(θj+10)+ωj+1−1π~j+1(θj+1M)+ωj+1. ∂ _j+1H_j+1≈ 1 π_j+1( _j+1^0)+ _j+1- 1 π_j+1( _j+1^M)+ _j+1\,\,. (15) Note that these denominators are calculated during the MCMC chain and therefore the derivative requires very little extra computation. A single ωj _j is kept for each layer and is maintained across subchains at that layer. We use a learning rate of 10−310^-3 to adjust ωj _j for all experiments. An update is made on layer j once after each M-step subchain. Algorithm 2 FLARE-with-layer-tuning(θj0θ^0_j,n,j) for i=0,…,n−1i=0,…,n-1 do if j=Jj=J then ⊳ coarsest layer Sample θ~ji θ^i_j from qj(⋅|θji)q_j(·|θ^i_j) Accept θji+1=θ~jiθ^i+1_j= θ^i_j with probability A from Equation 2 Otherwise, reject and θji+1=θjiθ^i+1_j=θ^i_j else θj+11,…,θj+1M=θ^1_j+1,…,θ^M_j+1= FLARE-with-layer-tuning(θjiθ^i_j, M, j+1j\!+\!1) θ~ji=θj+1M θ^i_j=θ^M_j+1 Update gradient of ωj+1 _j+1 using ∂ωj+1Hj+1 ∂ _j+1H_j+1 from Equation 15 Accept θji+1=θ~jiθ^i+1_j= θ^i_j with probability jA_j from Equation 11 Otherwise, reject and θji+1=θjiθ^i+1_j=θ^i_j return θj1,…,θjnθ^1_j,…,θ^n_j To make the innermost Gaussian proposal more robust, we adaptively update the covariance of the proposal distribution as initially proposed in the AM algorithm [27]. We use the history of chains from the coarsest layer θJ0,θJ1,…,θJtθ^0_J,θ^1_J,…,θ^t_J to update the covariance for the inner most proposal distribution. By using all previous states of the coarsest layer, the proposal distribution quickly adapts using the accepted samples. This rapid start of adaptation ensures good mixing in the inner most layer which gives higher quality candidate samples for the finer chains. We show the recursive algorithm with layer tuning adaptation added in Algorithm 2. Here we update the gradient after M steps of each layer and use it in the acceptance probability with ω mixed in as a uniform distribution to the target distribution. 3.5 Ergodicity of Layer Tuning We show FLARE MCMC with adaptive tuning of the proposals at each layer is ergodic. This can be shown with diminishing adaptation and simultaneous uniform ergodicity. Lemma 6. For layers 0≤j<J0≤ j<J, let γj∈Γj _j∈ _j be the adaptations for the proposal at layer j or the chain at level j+1j+1, i.e, γj=ωj+1↔ψj+1(θ)↔pj+1(θ→⋅) _j= _j+1 _j+1(θ) p_j+1(θ→·) where Γj∈ℝ _j and ψj+1 _j+1 is the target at layer j+1j+1 with the layer tuning adaptation added. Let pj,γj(θ→⋅)p_j, _j(θ→·) denote the transition distribution of chain at level j using adaptation γj _j, starting in state θ. Assume ∀j,ωj∈[ω¯,ω¯]∀ j, _j∈ [ ω, ω ] for some 0<ω¯<ω¯0< ω< ω, and, at the inner most layer, there exists a minorization constant ξJ>0 _ J>0 such that ‖pJ,γJM(θ→⋅)−ψJ(⋅)‖≤(1−ξJ)M \|p_ J, _ J^M(θ→·)- _ J(·) \|≤(1- _ J)^M. Then, (a) Simultaneous uniform ergodicity: For all τ>0τ>0. there exists n=n(τ)∈ℕn=n(τ) such that ‖pJ,γJn(θ→⋅)−ψj(⋅)‖≤τ \|p_ J, _ J^n(θ→·)- _j(·) \|≤τ (16) for all θ∈Θθ∈ and γj∈Γj _j∈ _j. (b) Diminishing adaptation: The amount of adaptation diminishes in probability with the number of steps t in the adaptation as limt→∞supθ‖pj,γjt(θ→⋅)−pj,γjt+1(θ→⋅)‖=0. _t→∞ _θ \|p_j, _j^t(θ→·)-p_j, _j^t+1(θ→·) \|=0. (17) Proof for Lemma 6 can be found in Appendix C. Theorem 7. FLARE MCMC with an adaptive layer tuning parameter is ergodic. Proof. We use Lemma 6 to show the conditions necessary in Theorem 1 of Roberts and Rosenthal (2007) [52]. This shows that the adaptive algorithm is ergodic. 4 Experiments We measure the efficiency of the MCMC methods tested using the effective sample size (ESS) [50] estimated across multiple chains as NESS=(N⋅K)/(1+2∑k=12m+1ρ(k))N_ESS= (N· K ) (1+2 _k=1^2m+1ρ(k) ) (18) where N is the number of samples, K is the number of chains, ρ(k)ρ(k) is the lag-k correlation, and m is the largest value such that ρ(2m)+ρ(2m+1)>0ρ(2m)+ρ(2m+1)>0. We compute ESS for each parameter for the “bulk” (entire distribution) and “tail” (largest and smallest 5%5\% of the samples) of the distributions. We compare FLARE MCMC with standard Metropolis-Hastings (with proposal adaptation introduced by Haario et al. [27]) and other multi-fidelity MCMC methods: Multilevel Delayed Acceptance MCMC (MLDA) [44], MLDA with Adaptive Error Model (AEM) [12, 44], Multi Level MCMC (MLMCMC) [14, 45], and Multi-fidelity Monte Carlo [8]. We give more detail about the different methods used for comparison: 1. MLDA without any adaptation: Introduced by Lykkegaard et al. [44], this method uses recursive chains of approximations as proposals. However, there is no adaptation being done to “correct” the approximations. Even the authors note that without adaptation, the chains do not mix well, and have poor effective sample sizes. 2. MLDA with AEM: Extended by Lykkegaard et al. [44] in the same paper, this method uses a similar structure to the above. They also use Adaptive Error Model (AEM) as a way to deal with the discrepancies between the different layers which uses a telescoping sum of differences in the mean of the approximations. They demonstrate with the subsurface flow model that MLDA with AEM leads to good mixing and high ESS. We demonstrate similar results for two of our experiments. Our subsurface flow model experiment uses the same fidelities as set up by the original authors; however, we run it to collect more samples using a higher number of chains. 3. Multi Level MCMC (MLMCMC): This method was proposed by Dodwell et al. [14] and was then applied to MLDA. A quantity of interest, Q, is proposed that is related to the parameters of the model. The samples drawn from the posterior are used to reduce the variance of Q. Since in MLDA, samples are not only drawn from a “true” posterior, but also approximations, the samples from the approximate levels are used to reduce the variance of Q. They state that it thus requires fewer samples to achieve the same variance. Using a telescopic sum, the difference of Q estimates between levels are used to correct Q with respect to the next coarser level. For the pendulum model, Q is the mean of the outputs. For the subsurface flow experiment used by the original authors of MLDA, Q is the hydraulic head at some fixed point (x,y)=(0.5,0.45)(x,y)=(0.5,0.45); that is, the model PDE is solved at these points at each level using samples from the coarser approximate level. 4. Multi-Fidelity Monte Carlo (MFMC): This method was proposed by Cai and Adams [8]. It uses a continuum of models with increasing fidelity and has a single Markov chain with a random choice of the fidelity, K, at each step. The fidelity K is part of the sampled state-space and therefore also part of the proposal distribution and acceptance probability. We map K to a reasonable range of fidelities for each experiment. For the pendulum model, we map K to the error tolerance of the integrator, ϵε, as ϵ=eK/10+10−6ε=e K10+10^-6. For the subsurface flow experiment, we let the grid resolution be equal to 10K10K (K is the sampled fidelity of this method) in order to map the resolution to the fidelity range expected by the algorithm’s implementation. The samples from this method are not from the true posterior, but rather can be corrected to estimate an expectation (like the mean). Therefore, we do not plot the evolution of effective samples with respect to time in the results, as they are not samples from the true posterior. Table 1: Mean ESS for bulk and tail distributions across 50 runs for 10 chains each, and mean value of each parameter with std error across all 500 total runs are listed. Average acceptance rates are listed per layer and the total samples are based on the average cost of likelihood evaluations per sample per method. Each chain is run for a total of 1000 seconds. α L acceptance rate total samples bulk ESS/s tail ESS/s mean sd bulk ESS/s tail ESS/s mean sd j=0j=0 j=1j=1 j=2j=2 MCMC 21.43 29.94 0.953 0.492 17.94 29.94 1.265 0.355 0.29 100000 FLARE MCMC(single) 52.47 58.42 1.081 0.019 45.85 54.95 1.372 0.016 0.98 0.24 60000 AEM MLDA (single) 39.41 53.91 1.064 0.062 42.20 52.24 1.371 0.051 0.98 0.27 50000 MLMCMC (single) 25.03 34.43 1.044 0.291 21.54 32.99 1.351 0.167 0.92 0.32 44500 MLDA (single) 11.64 1.49 1.059 0.041 15.36 7.73 1.361 0.022 0.91 0.31 45000 FLARE MCMC(double) 64.26 72.10 1.086 0.001 56.68 68.20 1.374 0.009 0.99 0.86 0.29 35000 AEM MLDA (double) 56.92 65.53 1.085 0.006 50.95 58.51 1.375 0.003 0.98 0.89 0.3 25000 MLMCMC (double) 33.82 37.05 1.049 0.015 33.09 39.94 1.361 0.015 0.90 0.81 0.28 25500 MLDA (double) 10.01 7.54 1.053 0.035 4.96 12.99 1.362 0.024 0.86 0.72 0.24 30000 MFMC 1.075 0.057 1.348 0.132 0.46 370000 (a) Simple Pendulum θ1 _1 θ2 _2 θ3 _3 acceptance rates total samples bulk ESS/s tail ESS/s mean sd bulk ESS/s tail ESS/s mean sd bulk ESS/s tail ESS/s mean sd j=0j=0 j=1j=1 j=2j=2 MCMC 5.35 7.81 -0.457 0.0037 5.49 8.04 0.466 0.0036 5.63 8.19 0.076 0.0034 0.27 10000 FLARE MCMC (single) 8.74 11.99 -0.460 0.0030 8.58 11.81 0.467 0.0036 8.63 11.58 0.076 0.0028 0.98 0.26 8365 AEM MLDA (single) 7.77 4.67 -0.459 0.0032 6.96 4.29 0.466 0.0033 8.16 6.24 0.076 0.0030 0.99 0.29 4100 MLMCMC (single) 4.35 5.86 -0.460 0.0028 4.37 3.27 0.465 0.0031 5.29 6.49 0.077 0.0021 0.93 0.32 8000 MLDA (single) 3.65 1.52 -0.463 0.0035 4.78 1.56 0.490 0.0036 3.33 3.02 0.075 0.0031 0.87 0.34 8250 FLARE MCMC (double) 15.49 21.35 -0.460 0.0026 15.141 19.69 0.469 0.0026 15.07 19.76 0.077 0.0023 0.99 0.95 0.3 6500 AEM MLDA (double) 13.29 14.70 -0.460 0.0035 12.32 13.22 0.468 0.0035 12.70 14.45 0.077 0.0033 0.99 0.93 0.24 2285 MLMCMC (double) 8.46 10.49 -0.459 0.0027 8.04 9.13 0.468 0.0030 7.47 5.81 0.077 0.0024 0.92 0.89 0.31 6000 MLDA (double) 4.48 5.05 -0.461 0.0036 5.33 6.60 0.469 0.0033 5.28 7.35 0.076 0.0030 0.93 0.91 0.27 6350 MFMC -0.475 0.026 0.418 0.041 0.102 0.003 0.39 5700 (b) Subsurface Flow model For the MLDA-based methods, we use the authors’ implementations in the open-source probabilistic programming package PyMC3 [53] by Lykkegaard et al. [44] For MFMC, we use the author-provided implementation. These implementations have significant computational overhead compared with our implementation of FLARE MCMC. Therefore, we only measure the time taken in likelihood computation (which is the same code for all methods). We present three different experimental posterior sampling problems across different scientific domains: a simple pendulum, a hydrology simulation that was used by prior methods as a benchmark, and a cosmology simulation that stresses computational limits. Because of the computational expense of the cosmology simulation, we are not able to collect a sufficient number of samples to get a reliable estimate of the effective sample sizes. Instead, we compare our estimates to those in the cosmology literature. For each experimental domain, we construct three fidelities by adjusting the relevant simulation parameter. In all cases, we measure our abilities to sample from the highest fidelity (j=0j\!\!=\!\!0). For methods labeled “(single),” there is a single higher fidelity layer (J=1J=1). For methods labeled “(double),” there are two higher fidelity layers (J=2J=2): the one from the “(single)” experiments, plus one more that is even more coarse. For the coarsest fidelity, a Gaussian proposal distribution is used with an adaptive covariance matrix. In the pendulum and cosmology experiments, this normal distribution is reflected to keep parameters within their respective ranges. FLARE MCMC can be extended beyond J=2J\!=\!2 layers. However, just two layers improves over the standard MCMC and other multi-level methods significantly. Layers’ costs should be roughly orders-of-magnitude different in computational costs. For these examples, J=2J\!=\!2 is the limit of how many layers can practically be constructed with orders-of-magnitude different computational costs. 4.1 Simple Pendulum The equation of motion for a pendulum of length L, mass M, and initial angle α0 _0 is α¨=−(g/L)sinα α=- ( gL ) α. Our goal is to sample from the posterior of the distribution the two parameters θ=(L,α)θ=(L,α) conditioned on the observations of α at three irregularly spaced times during the motion: α(1)=−0.85α(1)=-0.85, α(2.3)=0.9α(2.3)=0.9, α(5.0)=0.95α(5.0)=0.95. Observations of these angles are assumed to be corrupted by Gaussian noise with known standard deviation: σ=0.1σ=0.1. Different fidelities correspond to adjusting the error tolerance of an adaptive Runge-Kutta 4(5) ODE integrator [15] (10−310^-3 or 10−610^-6 in our experiments) with stepsize control and dense output [28]. As a separate coarsest layer of approximation, we use the small angle approximation (which does not hold for the observations), sin(α)≈α (α)≈α, reducing the equation of motion to a simple harmonic motion which can be solved analytically as α(t)=α0cos(tg/L)α(t)= _0 (t gL). The difference between the finest fidelity posterior and this small angle approximation is shown in Figure 2. 50001000015000200002500030000Bulk ESSinitial angle (α0 _0)length (L L)5101520M values50001000015000200002500030000Tail ESS5101520M valuesFLARE MCMC (single)MLDA AEM (single) Figure 3: ESS as a function of M for two single layered (J=1J=1) methods with fixed computation time of 500 seconds each ran across 50 different chains. Results To judge the importance of setting M, we evaluated our method across different values of M with fixed total computation time. The effective sample size (ESS) is plotted as a function of M for the two competing multilevel methods in Figure 3. We can see that beyond M=5M=5, the increased computation time from running longer inner subchains leads to a decrease in overall sampling efficiency, indicating that M=5M=5 offers the best trade-off between computation time and sampling efficiency. Therefore, the layered subchains were run for M=5M=5 steps. We ran 10 chains of the finest fidelity for all methods.We replicated this experiment (of 10 chains) 50 times. Table 1a summarizes the mean effective sample size per second (ESS/s) and the average mean of the parameters across all 500 chains with the standard deviation along with the acceptance rates for every layer. Figure 4 shows ESS (across all 10 chains) as a function of computation time for each method, with the total number of samples (N) generated in 1000 seconds. The standard deviations are plotted as (barely visible) vertical bars. For the sake of readability, we have separated our plots to show how each method performs with one level of nesting (single) and two levels of nesting (double). We note that the MFMC method obtains significantly more samples in the same time budget. This result arises from its randomized fidelity selection, which collects samples at low-fidelity evaluations more frequently than other methods. Since the samples from this method are not from the true posterior, we do not list the effective samples in the table. However, for comparison, the average bulk and tail ESS/s for 50 runs of 10 chains each measured for parameters [α,l][α,l] are [24.80,20.46][24.80,20.46] and [26.52,30.61][26.52,30.61] for MFMC. All methods are able to improve by using more fidelities. Our FLARE MCMC method is consistently and significantly better than the other methods (including the best one, MLDA with AEM as shown in the table) in terms of ESS/s in both the bulk and tail of the distribution. The acceptance rates for the multilevel methods indicate that while the coarsest level accepts about a third of the samples (consistent with Gelman et al. [21]), the proposed sample from that level is accepted by the finer levels frequently since it was already accepted by an approximate coarse level. The mean of the parameter across different runs of the standard MCMC has a higher standard deviation compared to the other multilevel methods suggesting that some runs of MCMC do a poor job at finding the modes in the posterior. FLARE MCMC produces more samples in the same amount of time compared to its competing method MLDA reflecting that our sampler requires fewer likelihood evaluations per step. 0 020000 2000040000 4000060000 60000Bulk ESSinitial angle (α0 _0)length (L L)250 250500 500750 7501000 1000time (seconds)0 020000 2000040000 4000060000 60000Tail ESS250 250500 500750 7501000 1000time (seconds)MCMCFLARE MCMC (single)MLDA AEM (single)MLMCMC (single)MLDA (single) (a) 0 025000 2500050000 5000075000 75000Bulk ESSinitial angle (α0 _0)length (L L)250 250500 500750 7501000 1000time (seconds)0 025000 2500050000 5000075000 75000Tail ESS250 250500 500750 7501000 1000time (seconds)MCMCFLARE MCMC (double)MLDA AEM (double)MLMCMC (double)MLDA (double) (b) Figure 4: Pendulum Model: ESS for bulk and tail across 50 runs (mean and std. dev.) for single and double layers of nesting. 10−1 10^-1100 10^0101 10^1102 10^2103 10^3time (seconds)10−2 10^-210−1 10^-1100 10^0Wasserstein distanceinitial angle (α0 _0)10−1 10^-1100 10^0101 10^1102 10^2103 10^3time (seconds)10−2 10^-210−1 10^-1100 10^0Wasserstein distancelength (L L)MCMCFLARE MCMC (single)FLARE MCMC (double)MLDA AEM (single)MLDA AEM (double)MLMCMC (single)MLMCMC (double)MLDA (single)MLDA (double) Figure 5: Wasserstein distance to true distributions for the pendulum model. Figure 5 shows the distributional distance to the true distribution as a function of number of samples, averaged over 500 chains. The true mean and standard deviation are unknown for all the experiments in the paper, and measuring the distance between a multi-dimensional distribution which can be evaluated (only up to a normalizing constant) and a distribution represented by samples is non-trivial. However, we have analyzed the 11-dimensional marginals of the pendulum model in the following way. We evaluate the true unnormalized distribution on a grid, normalize it, project it to the marginal of interest and then treat it as (weighted) samples for a sample-to-sample Wasserstein distance between it and the samples from the MCMC methods. As we refine the grid, the distances become smaller (for almost all methods). We refine the grid until these distances stabilize, resulting in about a 1000-by-1000 grid (1 million llh evaluations). From the figure, it is clear that for both the parameters, FLARE MCMC ends up with the smallest Wasserstein distance to the true distribution. 4.2 Estimation of Soil Permeability in Subsurface Flow We consider a simple problem in subsurface flow modeling [14]. This model was also used to evaluate the MLDA methods by the original authors [44], and we did not modify the code used by MLDA (except to increase the number of chains and measure time). The classical equations governing (steady state) single–phase subsurface flow consist of Darcy’s law coupled with an incompressibility condition: w+k∇p=g and ∇⋅w=0w+k∇ p=g and ∇· w=0 (19) subject to suitable boundary conditions. All quantities are fields over =[0,1]2⊂ℝ2D=[0,1]^2 ^2 for these experiments. Here p denotes the hydraulic head of the fluid, k is the permeability tensor, w is filtration velocity (or Darcy flux) and g is the (known) source term. 0 05000 500010000 10000Bulk ESSθ1 _1θ2 _2θ3 _3500 5001000 1000time (seconds)0 05000 500010000 10000Tail ESS500 5001000 1000time (seconds)500 5001000 1000time (seconds)MCMCFLARE MCMC (single)MLDA AEM (single)MLMCMC (single)MLDA (single) (a) 0 010000 1000020000 20000Bulk ESSθ1 _1θ2 _2θ3 _3500 5001000 1000time (seconds)0 010000 1000020000 20000Tail ESS500 5001000 1000time (seconds)500 5001000 1000time (seconds)MCMCFLARE MCMC (double)MLDA AEM (double)MLMCMC (double)MLDA (double) (b) Figure 6: Subsurface Flow Model: ESS for bulk and tail across 50 runs (mean and std. dev.) for single and double layers of nesting. We are interested in the permeability given observations (with known-variance Gaussian noise) of the hydraulic head at 16 regularly spaced points in D. k is simplified to be the gradient of a random scalar field. The log-Gaussian scalar field is parameterized with a truncated Karhunen-Loéve (KL) expansion (to three terms, following MLMCMC [14]). These three parameters (θ) have a standard normal prior and we sample from their posterior. Computing the likelihood involves solving a partial differential equation (PDE) with known boundary conditions for a given θ and comparing the results for p at the observation points. The fidelities correspond to different grid resolutions for the PDE solver: 120×120120× 120 (highest), 30×3030× 30, and 10×1010× 10 (coarsest). Results Previous work reports that M=5M=5 achieves the best trade-off between effective sample size and computation time for this experimental setup [44]. Therefore, we adopt the same value to ensure a fair comparison with our method. Table 1b summarizes the same statistics for this model with the same set-up as the pendulum experiments. Figure 6 shows ESS (across all 10 chains) as a function of computation time for each method. In terms of ESS/s, our method improves over the standard MCMC and outperforms the multilevel methods for the same amount of likelihood computational budget, especially in the tail of the distribution. All methods converge to similar means of the parameters with low standard deviation among chains. We note that MFMC collects fewer number of samples compared to other methods since each sample requires multiple log likelihood calculations in the same high fidelity to update K, leading to significant add up of computational costs. For MFMC, the estimated mean ESS/s for bulk and tail for parameters [θ1,θ2,θ3][ _1, _2, _3] are [0.128,0.081,0.710][0.128,0.081,0.710] and [0.227,0.161,1.107][0.227,0.161,1.107]. But, again, the samples from MFMC were never intended to be interpreted as from the true distribution. 4.3 Structure Formation in the Universe with N-body Gravitational Simulation An important problem in modern-day cosmology is to generate theoretical models of the Universe on very large scales (tens of MpcMpc across) that can be compared to observations. Bayesian inference allows cosmologists to measure quantities of fundamental physics significance, such as the nature of dark energy and dark matter [48]. The theoretical models needed for next generation telescopes, such as euclid [6] and the Roman Space Telescope (wfirst) [54], are based on expensive numerical simulations, some of which require many days of computer time for each evaluation. For such a computationally expensive model, we show the efficacy of FLARE MCMC as compared with the standard Metropolis-Hastings algorithm and MLDA. One of the most frequently used summary statistics is the galaxy power spectrum, g P_g (a bold P is a power spectrum, not a distribution): the two-point clustering of galaxies in Fourier space as a function of the wavenumber scale, k. We use a slightly simplified model for the galaxy power spectrum for (relative) ease of computation. We perform a forward simulation which starts from a given set of cosmological parameters and predicts the galaxy power spectrum. It works by following the evolution of the Universe under the influence of gravity, from its beginnings in an almost uniform density state to the diverse collection of galaxies sitting in dark matter potentials observed today. We sample from the posterior density of four cosmological parameters: θ1 _1 The dimensionless Hubble constant, h, which characterizes the Universal expansion rate and thus the recession velocity of distant galaxies. A redshift zero galaxy at distance d MpcMpc recedes at a speed v=H0dv=H_0d, where H0=h×100kms−1Mpc−1H_0=h× 100\,km\,s^-1Mpc^-1. Measuring h is of importance to understand dark energy. θ2 _2 The dimensionless total matter density, 0<Ω0<10< _0<1. Ω0 _0 is the energy density of matter as a function of the critical density. Ω0 _0 is important because it can be used to infer the density of dark matter. θ3 _3 The dimensionless scalar perturbation amplitude, AsA_s, of the primordial fluctuations at the wavenumber k=0.05Mpc−1k=0.05\,Mpc^-1. AsA_s is of interest because it connects to the uncertain high energy physics of the Early Universe. Larger values of AsA_s correspond to a clumpier early Universe and so lead to larger g P_g. θ4 _4 The dimensionless linear bias, b, which is used to shift the amplitude of our simulated matter power spectrum to match the amplitude of the galaxy power spectrum. This is to account for the difference between observed galaxies and dark matter (which is used by the forward model): model(θ)=b2⋅dm(h,Ω0,As), P_model(θ)=b^2· P_dm(h, _0,A_s), (20) where b is the scale-independent linear bias and dm P_dm is the simulated dark-matter power spectrum directly computed from the output density field of FastPM. The posterior density is conditioned on the galaxy power spectrum from SDSS-I Baryon Oscillation Spectroscopic Survey (BOSS) Data Release 12 (DR12) as our observational data source [13, 4]. We have used a subset of the BOSS data from the North Galactic Cap (NGC) at z=0.38z=0.38, which includes ∼106 10^6 galaxies, from Ivanov et al. [38]. The likelihood function is a multivariate Gaussian between the galaxy power spectrum from BOSS, g P_g, and the galaxy power spectrum from the forward model, model(θ) P_model(θ): lnℒ(θ)=−12(model(θ)−g)⊺−1(model(θ)−g)+k. (θ)=- 12( P_model(θ)- P_g) C^-1( P_model(θ)- P_g)+k. (21) C is the covariance matrix of the galaxy power spectrum, also estimated observationally. The most expensive part of the forward model, evolution under gravitational force, is computed using FastPM [19]. FastPM has a couple of tunable fidelity parameters. The size of the region simulated controls the amount of data available and may have a non-linear effect on the accuracy of the result. We thus fix the size of this region to 1024 Mpc/hMpc/h and instead change the number of particles. More particles in the simulation mean higher resolution, more accurate power spectrum at higher wavenumber k, and thus the likelihood function is higher fidelity. For N-body simulations, the compute time usually scales as NlogN N, where N is the number of particles. Thus a simulation with a 5123512^3 number of particles is ≃80 80 times more expensive than a 1283128^3 simulation. We therefore set the fidelities by only adjusting the number of particles used in the simulation to be 512512 (highest), 384384, and 256256 (coarsest). Note this calculation is distributed across 20 cores (using MPI) and therefore a saving of 1 hour corresponds to 20 core-hours. 102030400.000.050.10average absolute mean errorh h102030400.100.20Ω0 _010203040time (hours)0.100.150.20average absolute mean errorAs A_s10203040time (hours)0.100.200.30b bMCMCFLARE MCMC (single)FLARE MCMC (double)AEM MLDA (single)AEM MLDA (double) Figure 7: Average (across 10 chains) absolute error in mean estimates for the four parameters, as a function of total simulation time. In the allotted 48 hours, the chains each sampled 800 samples for (plain) MCMC, 600 samples for FLARE MCMC (single), 500 samples for FLARE MCMC (double), 290 samples for MLDA (single) and 200 samples for MLDA (double) We compare our estimated distributional means to previous computations on the same data. For the parameters h and Ω0 _0, we compare to the means reported by Ivanov et al. [37] on the same data using their own MCMC simulation (h=0.661h=0.661 and Ωm=0.290 _m=0.290). We have only a single linear bias term, compared with the multiple such terms of Ivanov et al. [37]. Therefore, we can compare neither b nor AsA_s (which is heavily related to b) to their results. Instead, we measure AsA_s against the best fit value from the Planck Satellite [2], As=2.09A_s=2.09 and b=2b=2, consistent with comparable BOSS measurements [37]. While these are modes (and not means), they are the best independent estimates we can obtain. 0500.70.8Valueh h0500.20.4Ω0 _0050231e−-9As A_s050123b b050time (hours)0.70.8Value050time (hours)0.20.4050time (hours)231e−-9050time (hours)123MCMCFLARE MCMC (double) Figure 8: Cosmology model; Trace plot for a random single run of each method Results Extreme running time dictated smaller values for M for this experiment. We reduced them by a factor of 2 (approximately) and used M=2M=2 for inner substeps. Figure 7 shows that the FLARE MCMC methods converge to the mean values from previous literature better than the standard Metropolis Hastings method using fewer samples and less time. The AsA_s parameter has slightly strange behavior. We can still see better convergence of our methods. However, note that the best-fit value of AsA_s we are taking as “ground truth” is measured (with error) from a different dataset, and thus is likely not the true mean of our posterior. Many large scale structure experiments prefer a lower value of this parameter than Planck, a feature known as the S8 tension [1]. The pairwise plots of the posterior can be found in the Appendix. From the posteriors, it is clear that FLARE MCMC is better at approximating the modes of the distribution as compared to MCMC. Figure 8 shows the trace plot for a random run of MCMC and FLARE MCMC. Our method shows better mixing than MCMC and is less likely to reject proposed samples. 5 Summary Many scientific and engineering problems involve simulations or solving differential equations. In this paper, we present an efficient multi-fidelity layered MCMC that exploits the ability to reduce the accuracy of models leading to approximations of the posterior. In a recursive, nested fashion, these approximations act as proposals for MCMC-based inference. We add layer tuning that successfully encourages the approximate proposals to explore the distribution well. We demonstrate with experimental results using models from three different scientific domains with varying costs that out method, FLARE MCMC, is simple, and yet produces more efficient samples than existing adaptive multilevel MCMC methods with the same computational budget. Appendix A Convergence Rate Proofs We show proofs for convergence rates in the main paper here. We first use Lemma 8 to show that after M steps of a coarse chain, we can obtain a minorized lower bound that can be recursively used in its finer layer. Lemma 8. Let pj(⋅→⋅)p_j(·→·) be the transition distribution of the Markov chain at level j with an invariant target distribution πj(⋅) _j(·). For any level j, if there exists a ξj>0 _j>0 such that pj1(θj0→θj1)≥ξjπj(θj1)p_j^1( _j^0→ _j^1)≥ _j _j( _j^1) for all θj0,θj1∈Θ _j^0, _j^1∈ , then pjM(θj0→θjM)≥(1−(1−ξj)M)πj(θjM)∀θj0,θjM∈Θp_j^M(θ^0_j→θ^M_j)≥ (1-(1- _j)^M ) _j(θ^M_j) ∀θ^0_j,θ^M_j∈ (22) Proof. We prove this using induction. We can verify the base case for k=1k=1 such that pj1(θj0→θj1)≥(1−(1−ξj)1)πj(θjk)=ξjπj(θj1)p_j^1(θ^0_j→θ^1_j)≥(1-(1- _j)^1) _j(θ^k_j)= _j _j(θ^1_j). This is held by the assumption made in the lemma. Assume using the induction hypothesis that, pjk(θj0→θjk)≥(1−(1−ξj)k)πj(θjk)p_j^k(θ^0_j→θ^k_j)≥(1-(1- _j)^k) _j(θ^k_j). We need to show that pjk+1(θj0→θjk+1)≥(1−(1−ξj)k+1)πj(θjk+1)p_j^k+1(θ^0_j→θ^k+1_j)≥(1-(1- _j)^k+1) _j(θ^k+1_j). Note pjk(θj0→θjk)p_j^k(θ^0_j→θ^k_j) can be written as pjk(θj0→θjk)=ϕkπj(θjk)+(1−ϕk)rj(θjk∣θj0)p_j^k(θ^0_j→θ^k_j)=φ^k _j(θ^k_j)+(1-φ^k)r_j(θ^k_j θ^0_j), where ϕkφ^k is the probability that the chain couples to the stationary distribution in k steps, and rj(θjk∣θj0)r_j(θ^k_j θ^0_j) is the remaining distribution that depends on θj0θ^0_j. pjk+1(θj0→θjk+1) p_j^k+1(θ^0_j→θ^k+1_j) =∫[pjk(θj0→θjk)⋅pj1(θjk→θjk+1)]dθjk = [p_j^k(θ^0_j→θ^k_j)· p_j^1(θ^k_j→θ^k+1_j) ]dθ^k_j =∫pj1(θjk→θjk+1)[ϕkπj(θjk)+(1−ϕk)rj(θjk∣θj0)]dθjk = p_j^1(θ^k_j→θ^k+1_j) [φ^k _j(θ^k_j)+(1-φ^k)r_j(θ^k_j θ^0_j) ]dθ^k_j =∫ϕkπj(θjk)pj1(θjk→θjk+1)dθjk+∫(1−ϕk)rj(θjk∣θj0)pj1(θjk→θjk+1)dθjk = φ^k _j(θ^k_j)p_j^1(θ^k_j→θ^k+1_j)dθ^k_j+ (1-φ^k)r_j(θ^k_j θ^0_j)p_j^1(θ^k_j→θ^k+1_j)dθ^k_j We know ϕk≥1−(1−ξj)kφ^k≥ 1-(1- _j)^k. Thus, ≥(1−(1−ξj)k)πj(θjk+1)+(1−ξj)k∫rj(θjk∣θj0)pj1(θjk→θjk+1)dθjk ≥(1-(1- _j)^k) _j(θ^k+1_j)+(1- _j)^k r_j(θ^k_j θ^0_j)p_j^1(θ^k_j→θ^k+1_j)dθ^k_j Replacing pj1(θjk→θjk+1)p_j^1(θ^k_j→θ^k+1_j) between two consecutive samples with the base assumption, ≥(1−(1−ξj)k)πj(θjk+1)+ξj(1−ξj)k∫rj(θjk∣θj0)πj(θjk+1)dθjk ≥(1-(1- _j)^k) _j(θ^k+1_j)+ _j(1- _j)^k r_j(θ^k_j θ^0_j) _j(θ^k+1_j)dθ^k_j =(1−(1−ξj)k)πj(θjk+1)+ξj(1−ξj)kπj(θjk+1) =(1-(1- _j)^k) _j(θ^k+1_j)+ _j(1- _j)^k _j(θ^k+1_j) =(1−(1−ξj)k+1)πj(θjk+1) =(1-(1- _j)^k+1) _j(θ^k+1_j) Therefore using proof by induction, we have that pjM(θj0→θjM)≥(1−(1−ξj)M)πj(θjM)p_j^M(θ^0_j→θ^M_j)≥(1-(1- _j)^M) _j(θ^M_j). Proof of Lemma 3. The transition kernel is given by pj1(θji→θji+1) p_j^1(θ^i_j→θ^i+1_j) =j(θji→θji+1)⋅qj(θji+1|θji)+δ(θji+1−θji)∫(1−j(θji→θj′))qj(θj′|θji)dθj′ =A_j(θ^i_j→θ^i+1_j)· q_j(θ^i+1_j|θ^i_j)+δ(θ^i+1_j-θ^i_j) (1-A_j(θ^i_j→θ _j) )q_j(θ _j|θ^i_j)\,dθ _j ≥j(θji→θji+1)⋅qj(θji+1|θji) _j(θ^i_j→θ^i+1_j)· q_j(θ^i+1_j|θ^i_j) =j(θji→θji+1)⋅pj+1M(θj+10→θj+1M) =A_j(θ^i_j→θ^i+1_j)· p_j+1^M(θ^0_j+1→θ^M_j+1) The sample, θji+1 _j^i+1 is proposed using the Mth sample from the j+1j+1 chain, therefore is the same as θj+1Mθ^M_j+1. Thus, from Lemma 8, ≥j(θji→θji+1)⋅(1−(1−ξj+1)M)⋅πj+1(θji+1) _j(θ^i_j→θ^i+1_j)·(1-(1- _j+1)^M)· _j+1(θ^i+1_j) =min(1,πj(θji+1)πj(θji)⋅πj+1(θji)πj+1(θji+1))⋅(1−(1−ξj+1)M)⋅πj+1(θji+1) = (1, _j(θ^i+1_j) _j(θ^i_j)· _j+1(θ^i_j) _j+1(θ^i+1_j) )·(1-(1- _j+1)^M)· _j+1(θ^i+1_j) Let r(θ)=πj+1(θ)πj(θ)r(θ)= _j+1(θ) _j(θ). Then, =min(1,r(θji)r(θji+1))⋅(1−(1−ξj+1)M)⋅r(θji+1)⋅πj(θji+1) = (1, r(θ^i_j)r(θ^i+1_j) )·(1-(1- _j+1)^M)· r(θ^i+1_j)· _j(θ^i+1_j) =(1−(1−ξj+1)M)⋅min(r(θji+1),r(θji))⋅πj(θji+1) =(1-(1- _j+1)^M)· (r(θ^i+1_j),r(θ^i_j) )· _j(θ^i+1_j) ≥(1−(1−ξj+1)M)⋅minθ(r(θ))⋅πj(θji+1) ≥(1-(1- _j+1)^M)· _θ (r(θ) )· _j(θ^i+1_j) =ξj⋅πj(θji+1) = _j· _j(θ^i+1_j) where ξj=(1−(1−ξj+1)M)⋅minθ(πj+1(θ)πj(θ)) _j=(1-(1- _j+1)^M)· _θ ( _j+1(θ) _j(θ) ). Appendix B Optimal M proof Proof of Lemma 5. Since the ratio minθ(πj+1(θ)πj(θ)) _θ ( _j+1(θ) _j(θ) ) is constant with respect to MjM_j, we simplify the objective function to be maximized as: f(Mj)=1−(1−ξj+1)Mjbj+MjBj+1,Mj≥0,ξj+1∈(0,1),bj,Bj+1>0.f(M_j)= 1-(1- _j+1)^M_jb_j+M_jB_j+1, M_j≥ 0, _j+1∈(0,1), b_j,B_j+1>0. We define c:=1−ξj+1∈(0,1),β:=bjBj+1,Υ:=−logc>0andμ:=β+1Υ.c:=1- _j+1∈(0,1), β:= b_jB_j+1, :=- c>0 μ:=β+ 1 . Then, f(Mj)=1−cMjβBj+1+MjBj+1=1−cMjBj+1(β+Mj).f(M_j)= 1-c^M_jβ B_j+1+M_jB_j+1= 1-c^M_jB_j+1(β+M_j). Differentiating the objective function, f′(Mj)=−cMjlogc(β+Mj)−(1−cMj)(β+Mj)2∗Bj+1.f (M_j)= -c^M_j c\,(β+M_j)-(1-c^M_j)(β+M_j)^2*B_j+1. Setting the derivative to zero and using Υ=−logc =- c, −cMjlogc(β+Mj)=1−cMj⟹cMj(1+Υ(β+Mj))=1.-c^M_j c\,(β+M_j)=1-c^M_j c^M_j (1+ (β+M_j) )=1. Since cMj=elog(cMj)=eMjlogc=eMj(−Υ)c^M_j=e (c^M_j)=e^M_j c=e^M_j(- ), we replace cMj=e−ΥMjc^M_j=e^- M_j to obtain (1+Υ(β+Mj))e−ΥMj=1. (1+ (β+M_j) )e^- M_j=1. By the definition of μ, we have 1+Υ(β+Mj)=Υμ+ΥMj1+ (β+M_j)= μ+ M_j. Therefore, (Υμ+ΥMj)e−ΥMj=1.( μ+ M_j)e^- M_j=1. Let y=Υμ+ΥMjy= μ+ M_j. Then, ΥMj=y−Υμ M_j=y- μ, substituting this gives ye−(y−Υμ) ye^-(y- μ) =1 =1 ye−y ye^-y =e−Υμ =e^- μ −ye−y -ye^-y =−e−Υμ. =-e^- μ. Using the Lambert W function for branch k=−1k=-1, y=−W−1(−e−Υμ).y=-W_-1 (-e^- μ ). (23) Since ΥMj=y−Υμ M_j=y- μ, we get Mj=yΥ−μ=−1ΥW−1(−e−Υμ)−μ.M_j= y -μ=- 1 W_-1 (-e^- μ )-μ. Therefore, the maximizer for our objective function is Mj∗=−1ΥW−1(−e−Υμ)−μM_j^*=- 1 W_-1 (-e^- μ )-μ (24) where Υ=−log(1−ξj+1) =- (1- _j+1), μ=bjBj+1+1Υμ= b_jB_j+1+ 1 and W−1W_-1 is −1-1 branch of the Lambert W function. Appendix C Layer Tuning Ergodicity Proof Proof of Lemma 6. (a) Consider layer J−1J-1. Using Theorem 4, ‖pJ−1,γJ−1M(θ→⋅)−ψJ−1(⋅)‖ \|p_ J-1, _ J-1^M(θ→·)- _ J-1(·) \| ≤(1−ξJ−1)M,ξJ−1=(1−(1−ξJ)M)minθψJ(θ)ψJ−1(θ) ≤(1- _ J-1)^M, _ J-1=(1-(1- _ J)^M) _θ _ J(θ) _ J-1(θ) From the definition of layer tuning, =(1−(1−ξJ)M)minθ(π~J(θ)+ωJ)⋅ζJ−1(ωJ−1)(π~J−1(θ)+ωJ−1)⋅ζJ(ωJ) =(1-(1- _ J)^M) _θ ( π_ J(θ)+ _ J)· _ J-1( _ J-1)( π_ J-1(θ)+ _ J-1)· _ J( _ J) =(1−(1−ξJ)M)minθπ~J(θ)+ωJπ~J−1(θ)+ωJ−1 =(1-(1- _ J)^M) _θ π_ J(θ)+ _ J π_ J-1(θ)+ _ J-1 ×Z+ωJ−1⋅VZ+ωJ⋅V × Z+ _J-1· VZ+ _J· V where ζ is the normalizing constant of the new distribution that depends on ω, Z is the normalizing constant of the original distribution, and V is the volume of Θ . With the bounds for ω from the assumption, ≥(1−(1−ξJ)M)minθ(π~J(θ)+ω¯π~J−1(θ)+ω¯)(Z+ω¯⋅VZ+ω¯⋅V) ≥(1-(1- _ J)^M) _θ ( π_ J(θ)+ ω π_ J-1(θ)+ ω ) ( Z+ ω· VZ+ ω· V ) ≜ξ¯J−1 ξ_ J-1 By induction with base case at layer J, ∀j,‖pj,γjn(θ→⋅)−ψj(⋅)‖≤(1−ξj)n where ξj≥ξ¯j.∀ j, \|p_j, _j^n(θ→·)- _j(·) \|≤(1- _j)^n where _j≥ ξ_j\,\,. (25) We need to show that for all τ>0τ>0. there exists n=n(τ)∈ℕn=n(τ) such that ‖pj,γjn(θ→⋅)−ψj(⋅)‖≤τ \|p_j, _j^n(θ→·)- _j(·) \|≤τ for all θ∈jθ _j and γj∈Γj _j∈ _j. From Equation 25, we want (1−ξj)n (1- _j)^n ≤τ ≤τ n n ≥lnτln(1−ξj) ≥ τ (1- _j) Since ln(1−ξj)≤ln(1−ξj¯) (1- _j)≤ (1- _j), n n ≥lnτln(1−ξj¯). ≥ τ (1- _j)\,\,. Thus, for all τ>0τ>0, there exists n=maxτlnτln(1−ξj¯)n= _τ τ (1- _j) such that ‖pJ,γJn(θ→⋅)−ψj(⋅)‖≤τ \|p_ J, _ J^n(θ→·)- _j(·) \|≤τ. (b) At every step t, the change in γj _j maps to change in ωj+1 _j+1. Diminishing adaptation is guaranteed by a gradient descent algorithm with diminishing stepsize that updates ωj+1 _j+1 at each layer to minimize the Kullback-Leibler divergence between layers ψj _j and ψj+1 _j+1. At each step of the GD algorithm, ωj+1 _j+1 is updated as ωj+1t+1=ωj+1t−ηi∂ωj+1Hj+1 _j+1^t+1= _j+1^t- _i ∂ _j+1H_j+1. To get diminishing adaptation, the update needs to converge as limt→∞‖ηt∂ωj+1Hj+1‖≈0. _t→∞ \| _t ∂ _j+1H_j+1 \|≈ 0. Since we bound γj↔ωj+1 _j _j+1 away from zero, ∂ωj+1Hj+1 ∂ _j+1H_j+1 =1π~j+1(θj+10)+ωj+1−1π~j+1(θj+1M)+ωj+1 = 1 π_j+1( _j+1^0)+ _j+1- 1 π_j+1( _j+1^M)+ _j+1 ≤1π~j+1(θj+10)+ωj+1 ≤ 1 π_j+1( _j+1^0)+ _j+1 ≤1π~j+1(θj+10)+ω¯ ≤ 1 π_j+1( _j+1^0)+ ω ≤1ω¯ ≤ 1 ω Therefore, the update is limt→∞‖ηt∂ωj+1Hj+1‖ _t→∞ \| _t ∂ _j+1H_j+1 \| ≤limt→∞‖ηt1ω¯‖ ≤ _t→∞ \| _t 1 ω \| ≤1ω¯limt→∞ηt ≤ 1 ω _t→∞ _t If the stepsize, ηt _t asymptotes to 00, adapation decreases to 00 as t→∞t→∞. Therefore limt→∞‖ωj+1t−ωj+1t+1‖=0 _t→∞ \| _j+1^t- _j+1^t+1 \|=0, and thus limt→∞supθ‖pj,γjt(θ→⋅)−pj,γjt+1(θ→⋅)‖=0 _t→∞ _θ \|p_j, _j^t(θ→·)-p_j, _j^t+1(θ→·) \|=0. Appendix D Uniform Smoothing Parameter Plotted in Figure 9 is the evolution of our tuning parameter ωj _j as a function of samples collected for one example chain of the doubly nested method for the pendulum model. Each sample at j=1j=1 starts a chain of length M=5M=5 at inner layer j=2j=2. The inner layer j=2j=2 uses the small angle approximation of the pendulum as the fidelity. Since it is a poor approximation of the posterior as shown in Fig 2, we start with a relatively high value of ω. This helps the coarsest layer better explore the high-probability regions. As shown, the tuning parameter ω converges close to zero after a few samples in both the layers. We use a learning rate of 10−310^-3 for both the layers. 0 05000 500010000 10000samples0.00 0.000.02 0.020.04 0.040.06 0.060.08 0.080.10 0.10ω1 _1j = 10 020000 2000040000 40000samples0.00 0.000.25 0.250.50 0.500.75 0.751.00 1.001.25 1.251.50 1.501.75 1.752.00 2.00ω2 _2j = 2 Figure 9: Evolution of ωj _j for doubly nested layers j=1j=1 and j=2j=2 Appendix E Computational Infrastructure Our experiments were performed on a machine with 4 Intel® Xeon® Silver 4214 CPUs running at 2.20GHz for our experiments (a total of 48 cores). The machine has 250GB of memory, but memory was never a restriction during our experiments. All methods use multiproccessing, that is, each chain is run in parallel using a different core. The time listed is across one run of a single chain; however, the effective sample size is calculated across 10 different chains. For the pendulum and hydrology models, likelihood calculations were carried out on a single core. For the cosmology model, the likelihood calculations were carried out in parallel across 20 cores. Therefore, for the cosmology experiments, saving a day’s worth of computation time on the graphs corresponds to saving 20 days worth of core-hours. Appendix F Pairwise Plot for Cosmology Model Plotted in Figure 10 is the pairwise plot for all four parameters of the cosmology model. MCMC generates more samples in the same period of time, yet these samples have not yet converged to the distribution and are still scattered across the space, compared with the relatively compact FLARE MCMC samples. Figure 10: Cosmology model: Pairwise plots for all four parameters. Plotted for a random run of samples collected for 48 hours of MCMC (black) and FLARE MCMC (green). References [1] E. Abdalla, G. F. Abellán, A. Aboubrahim, A. Agnello, Ö. Akarsu, et al. (2022) Cosmology intertwined: a review of the particle physics, astrophysics, and cosmology associated with the cosmological tensions and anomalies. Journal of High Energy Astrophysics 34, p. 49–211. External Links: ISSN 2214-4048, Link, Document Cited by: §4.3. [2] N. Aghanim, Y. Akrami, M. Ashdown, J. Aumont, C. Baccigalupi, et al. (2020) Planck 2018 results. Astronomy & Astrophysics 641, p. A6. External Links: Document, Link Cited by: §4.3. [3] F. Al-Awadhi, M. Hurn, and C. Jennison (2004) Improving the acceptance rate of reversible jump mcmc proposals. Statistics & Probability Letters 69, p. 189–198. External Links: Document Cited by: §2.2. [4] S. Alam, M. Ata, S. Bailey, F. Beutler, Bizyaev, et al. (2017) The clustering of galaxies in the completed SDSS-I Baryon Oscillation Spectroscopic Survey: cosmological analysis of the DR12 galaxy sample. MNRAS 470 (3), p. 2617–2652. External Links: Document, 1607.03155 Cited by: §1, §4.3. [5] G. Altekar, S. Dwarkadas, J. P. Huelsenbeck, and F. Ronquist (2004) Parallel metropolis coupled markov chain monte carlo for bayesian phylogenetic inference. Bioinformatics 20 (3), p. 407–415. External Links: Document Cited by: §2.2. [6] L. Amendola, S. Appleby, A. Avgoustidis, D. Bacon, T. Baker, et al. (2018) Cosmology and fundamental physics with the Euclid satellite. Living Reviews in Relativity 21 (1), p. 2. External Links: Document, 1606.00180 Cited by: §4.3. [7] A. Beskos, A. Jasra, K. Law, R. Tempone, and Y. Zhou (2017) Multilevel sequential monte carlo samplers. Stochastic Processes and their Applications 127 (5), p. 1417–1440. Cited by: §2.2. [8] D. Cai and R. P. Adams (2022) Multi-fidelity Monte Carlo: a pseudo-marginal approach. In Advances in Neural Information Processing Systems, S. Koyejo, S. Mohamed, A. Agarwal, D. Belgrave, K. Cho, and A. Oh (Eds.), Vol. 35, p. 21654–21667. External Links: Link Cited by: §2.2, item 4, §4. [9] J. A. Christen and C. Fox (2005) Markov chain monte carlo using an approximation. Journal of Computational and Graphical Statistics 14 (4), p. 795–810. External Links: Document, Link, https://doi.org/10.1198/106186005X76983 Cited by: §2.2. [10] P. R. Conrad, Y. M. Marzouk, N. S. Pillai, and A. Smith (2016) Accelerating asymptotically exact mcmc for computationally intensive models via local approximations. Journal of the American Statistical Association 111 (516), p. 1591–1607. External Links: ISSN 1537-274X, Link, Document Cited by: §2.2. [11] T. Cui, C. Fox, and M. J. O’Sullivan (2019) A posteriori stochastic correction of reduced models in delayed-acceptance MCMC, with application to multiphase subsurface inverse problems. International Journal for Numerical Methods in Engineering 118 (10), p. 578–605 (English). External Links: Document, ISSN 0029-5981 Cited by: §2.2. [12] T. Cui, C. Fox, and M. O’Sullivan (2012) Adaptive error modelling in MCMC sampling for large scale inverse problems. Cited by: §2.2, §4. [13] K. S. Dawson, D. J. Schlegel, C. P. Ahn, S. F. Anderson, Aubourg, et al. (2013) The Baryon Oscillation Spectroscopic Survey of SDSS-I. AJ 145 (1), p. 10. External Links: Document, 1208.0022 Cited by: §1, §4.3. [14] T. J. Dodwell, C. Ketelsen, R. Scheichl, and A. L. Teckentrup (2015) A hierarchical multilevel Markov chain Monte Carlo algorithm with applications to uncertainty quantification in subsurface flow. SIAM/ASA Journal on Uncertainty Quantification 3 (1), p. 1075–1108. External Links: Document, Link, https://doi.org/10.1137/130915005 Cited by: §2.2, item 3, §4.2, §4.2, §4. [15] J.R. Dormand and P.J. Prince (1980) A family of embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics 6 (1), p. 19–26. External Links: ISSN 0377-0427, Document, Link Cited by: §4.1. [16] A. Doucet, N. de Freitas, and N. Gordon (2001) An introduction to sequential monte carlo methods. In Sequential Monte Carlo Methods in Practice, A. Doucet, N. de Freitas, and N. Gordon (Eds.), p. 3–14. External Links: ISBN 978-1-4757-3437-9, Document, Link Cited by: §2.2. [17] S. Duane, A.D. Kennedy, B. J. Pendleton, and D. Roweth (1987) Hybrid Monte Carlo. Physics Letters B 195 (2), p. 216–222. Cited by: §1. [18] Y. Efendiev, T. Hou, and W. Luo (2006) Preconditioning markov chain monte carlo simulations using coarse-scale models. SIAM Journal on Scientific Computing 28 (2), p. 776–803. External Links: Document, Link, https://doi.org/10.1137/050628568 Cited by: §2.2. [19] Y. Feng, M. Chu, U. Seljak, and P. McDonald (2016) FASTPM: a new scheme for fast simulations of dark matter and haloes. MNRAS 463 (3), p. 2273–2286. External Links: Document, 1603.00476 Cited by: §4.3. [20] C. Fox and G. Nicholls (1997) Sampling conductivity images via MCMC. In The Art and Science of Bayesian Image Analysis, p. 91–100. Cited by: §2.2. [21] A. Gelman, G. O. Roberts, and W. R. Gilks (1996) Efficient Metropolis jumping rules. In Bayesian Statistics, J. M. Bernardo, J. O. Berger, A. P. Dawid, and A. F. M. Smith (Eds.), p. 599–608. Cited by: §4.1. [22] M. B. Giles (2008) Multilevel Monte Carlo path simulation. Operations Research 56 (3), p. 607–617. External Links: ISSN 0030364X, 15265463, Link Cited by: §2.2. [23] P. J. Green and A. Mira (2001) Delayed rejection in reversible jump Metropolis-Hastings. Biometrika 88 (4), p. 1035–1053. External Links: ISSN 00063444, Link Cited by: §2.2. [24] P. J. Green (1995) Reversible jump Markov chain Monte Carlo computation and Bayesian model determination. Biometrika 82 (4), p. 711–732. External Links: ISSN 0006-3444, Document, Link, https://academic.oup.com/biomet/article-pdf/82/4/711/699533/82-4-711.pdf Cited by: §2.2. [25] A. Gregory, C. J. Cotter, and S. Reich (2016) Multilevel ensemble transform particle filtering. SIAM Journal on Scientific Computing 38 (3), p. A1317–A1338. Cited by: §2.2. [26] A. Gregory and C. J. Cotter (2017) A seamless multilevel ensemble transform particle filter. SIAM Journal on Scientific Computing 39 (6), p. A2684–A2701. Cited by: §2.2. [27] H. Haario, E. Saksman, and J. Tamminen (2001) An adaptive Metropolis algorithm. Bernoulli 7 (2), p. 223 – 242. Cited by: §3.4, §4. [28] E. Hairer, S.P. Nørsett, and G. Wanner (2000) Solving ordinary differential equations I nonstiff problems. Second edition, Springer, Berlin. Cited by: §4.1. [29] A. Haji-Ali, F. Nobile, L. Tamellini, and R. Tempone (2016) Multi-index stochastic collocation convergence rates for random pdes with parametric regularity. Foundations of Computational Mathematics 16 (6), p. 1555–1605. Cited by: §2.2. [30] W. K. Hastings (1970) Monte Carlo sampling methods using Markov chains and their applications. Biometrika 57 (1), p. 97–109. External Links: ISSN 00063444, Link Cited by: §2.1. [31] S. Heinrich (2001) Multilevel Monte Carlo methods. In Proceedings of the Third International Conference on Large-Scale Scientific Computing-Revised Papers, LSSC ’01, Berlin, Heidelberg, p. 58–67. External Links: ISBN 3540430431 Cited by: §2.2. [32] D. Higdon, H. Lee, and Z. Bi (2002) A bayesian approach to characterizing uncertainty in inverse problems using coarse and fine-scale information. IEEE Transactions on Signal Processing 50 (2), p. 389–399. External Links: Document Cited by: §2.2. [33] G. E. Hinton (2002) Training products of experts by minimizing contrastive divergence. Neural computation 14 (8), p. 1771–1800. Cited by: §3.4. [34] V. H. Hoang, C. Schwab, and A. M. Stuart (2013) Complexity analysis of accelerated mcmc methods for bayesian inversion. Inverse Problems 29 (8), p. 085010. Cited by: §2.2. [35] H. Hoel, K. J. Law, and R. Tempone (2016) Multilevel ensemble kalman filtering. SIAM Journal on Numerical Analysis 54 (3), p. 1813–1839. Cited by: §2.2. [36] M. D. Hoffman and A. Gelman (2014) The no-U-turn sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research 15, p. 1593–1623. Cited by: §1. [37] M. M. Ivanov, M. Simonović , and M. Zaldarriaga (2020) Cosmological parameters from the BOSS galaxy power spectrum. Journal of Cosmology and Astroparticle Physics 2020 (05), p. 042–042. External Links: Document, Link Cited by: §4.3. [38] M. M. Ivanov, M. Simonović, and M. Zaldarriaga (2020) Cosmological parameters from the BOSS galaxy power spectrum. J. Cosmology Astropart. Phys 2020 (5), p. 042. External Links: Document, 1909.05277 Cited by: §4.3. [39] A. Jasra, K. Kamatani, K. J. Law, and Y. Zhou (2017) Multilevel particle filters. SIAM Journal on Numerical Analysis 55 (6), p. 3068–3096. Cited by: §2.2. [40] A. Jasra, K. Kamatani, K. J. Law, and Y. Zhou (2018) A multi-index markov chain monte carlo method. International Journal for Uncertainty Quantification 8 (1). Cited by: §2.2. [41] A. Jasra, K. Kamatani, K. Law, and Y. Zhou (2018) Bayesian static parameter estimation for partially observed diffusions via multilevel monte carlo. SIAM Journal on Scientific Computing 40 (2), p. A887–A902. Cited by: §2.2. [42] J. Kaipio and E. Somersalo (2007) Statistical inverse problems: discretization, model reduction and inverse crimes. Journal of Computational and Applied Mathematics 198 (2), p. 493–504. Note: Special Issue: Applied Computational Inverse Problems External Links: ISSN 0377-0427, Document, Link Cited by: §2.2. [43] J. Liu and R. Chen (1998) Sequential monte carlo methods for dynamic systems. Journal of the American Statistical Association 93, p. . External Links: Document Cited by: §2.2. [44] M. B. Lykkegaard, T. J. Dodwell, C. Fox, G. Mingas, and R. Scheichl (2023) Multilevel delayed acceptance mcmc. SIAM/ASA Journal on Uncertainty Quantification 11 (1), p. 1–30. External Links: Document, Link, https://doi.org/10.1137/22M1476770 Cited by: §2.2, item 1, item 2, §4.2, §4.2, §4, §4. [45] M. B. Lykkegaard, G. Mingas, R. Scheichl, C. Fox, and T. J. Dodwell (2020) Multilevel delayed acceptance MCMC with an adaptive error model in PyMC3. arXiv. External Links: Document, Link Cited by: §4. [46] E. Marinari and G. Parisi (1992) Simulated tempering: a new Monte Carlo scheme. Europhysics Letters (EPL) 19 (6), p. 451–458. External Links: ISSN 0295-5075, Document, Link Cited by: §2.2. [47] R. M. Neal (1996) Bayesian learning for neural networks. Springer New York, NY. Cited by: §1. [48] P. J. E. Peebles (1980) The large-scale structure of the universe. Cited by: §4.3. [49] B. Peherstorfer, K. Willcox, and M. Gunzburger (2018) Survey of multifidelity methods in uncertainty propagation, inference, and optimization. SIAM Review 60 (3), p. 550–591 (English (US)). External Links: Document, ISSN 0036-1445 Cited by: §2.2. [50] B. D. Ripley (1987) Stochastic simulation. John Wiley & Sons, Inc., New York, NY, USA. External Links: ISBN 0-471-81884-4 Cited by: §4. [51] G. O. Roberts and J. S. Rosenthal (2004) General state space Markov chains and MCMC algorithms. Probability Surveys 1 (none). External Links: Document, Link Cited by: §3.2, §3.2, Definition 2. [52] G. O. Roberts and J. S. Rosenthal (2007) Coupling and ergodicity of adaptive markov chain monte carlo algorithms. Journal of Applied Probability 44 (2), p. 458–475. External Links: ISSN 00219002, Link Cited by: §3.5. [53] J. Salvatier, T. V. Wiecki, and C. Fonnesbeck (2016) Probabilistic programming in python using PyMC3. PeerJ Computer Science 2, p. e55. External Links: Document, Link Cited by: §4. [54] D. Spergel, N. Gehrels, J. Breckinridge, M. Donahue, A. Dressler, B. S. Gaudi, et al. (2013) Wide-Field InfraRed Survey Telescope-Astrophysics Focused Telescope Assets WFIRST-AFTA Final Report. arXiv e-prints, p. arXiv:1305.5422. External Links: 1305.5422 Cited by: §4.3. [55] H. Strasser (1985) Mathematical theory of statisticsMathematical theory of statistics. De Gruyter, Berlin, New York. External Links: Link, Document, ISBN 9783110850826 Cited by: Definition 1. [56] I. Sutskever and T. Tieleman (2010) On the convergence properties of contrastive divergence. In Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, Y. W. Teh and M. Titterington (Eds.), Proceedings of Machine Learning Research, Vol. 9, Chia Laguna Resort, Sardinia, Italy, p. 789–795. External Links: Link Cited by: §3.4. [57] R. H. Swendsen and J. Wang (1986) Replica monte carlo simulation of spin-glasses. Phys. Rev. Lett. 57, p. 2607–2609. External Links: Document, Link Cited by: §2.2. [58] L. Tierney and A. Mira (1999) Some adaptive Monte Carlo methods for Bayesian inference.. Statistics in medicine 18 17-18, p. 2507–15. Cited by: §2.2.