Stratification as a general variance reduction method for Markov chain Monte Carlo

TL;DR

This paper presents a theoretical analysis of EMUS, a stratified sampling method, demonstrating its efficiency in computing averages over complex distributions, especially in multimodal or tail probability scenarios.

Contribution

The paper provides a detailed theoretical foundation for EMUS, establishing it as a versatile and more efficient alternative to traditional MCMC in challenging sampling problems.

Findings

EMUS can outperform MCMC in multimodal distributions.

EMUS is effective for tail probability estimation.

Theoretical analysis confirms EMUS's efficiency.

Abstract

The Eigenvector Method for Umbrella Sampling (EMUS) belongs to a popular class of methods in statistical mechanics which adapt the principle of stratified survey sampling to the computation of free energies. We develop a detailed theoretical analysis of EMUS. Based on this analysis, we show that EMUS is an efficient general method for computing averages over arbitrary target distributions. In particular, we show that EMUS can be dramatically more efficient than direct MCMC when the target distribution is multimodal or when the goal is to compute tail probabilities. To illustrate these theoretical results, we present a tutorial application of the method to a problem from Bayesian statistics.

Equations

$$p_M:=\mathbf{P}[X\ge M].$$

$$\pi [g]:={\int }_{\Omega }g(x)\pi (dx),$$

$${\pi }_i(dx):=\frac{{\psi }_i(x)\pi (dx)}{\pi [{\psi }_i]}$$

$$\sum _{i=1}^L{\psi }_i(x)>0\text{ for all }x\in \Omega .$$

$$h^{\ast }(x):=\frac{h(x)}{{\sum }_{k=1}^L{\psi }_k(x)/u_k}.$$

$$\pi [g]=\frac{{\sum }_{i=1}^L{z}_i{\pi }_i[g^{\ast }]/u_i}{{\sum }_{i=1}^L{z}_i{\pi }_i[1^{\ast }]/u_i},$$

$$z_i=\frac{\pi [{\psi }_i]}{{\sum }_{k=1}^L\pi [{\psi }_k]}.$$

$$w^{\mathrm{t}}F=w^{\mathrm{t}},\text{ where }{F}_{ij}:={\pi }_i[{\psi }_j^{\ast }]/u_j.$$

$$w(G{)}^{\mathrm{t}}G=w(G{)}^{\mathrm{t}}\text{ with }\sum _{i=1}^L{w}_i(G)=1.$$

$$\pi [g]=\frac{{\sum }_{i=1}^L{w}_i(F){\pi }_i[g^{\ast }]}{{\sum }_{i=1}^L{w}_i(F){\pi }_i[1^{\ast }]}.$$

$${\bar{g}}_i^{\ast }:=\frac{1}{N_i}\sum _{t=1}^{N_i}{g}^{\ast }(X_t^i),{\bar{1}}_i^{\ast }:=\frac{1}{N_i}\sum _{t=1}^{N_i}{1}^{\ast }(X_t^i),\text{ and }{\bar{F}}_{ij}:=\frac{1}{N_i}\sum _{t=1}^{N_i}{\psi }_j^{\ast }(X_t^i)/u_j.$$

$${\pi }_{\text{US}}[g]:=\frac{{\sum }_{i=1}^L{w}_i(\bar{F}){\bar{g}}_i^{\ast }}{{\sum }_{i=1}^L{w}_i(\bar{F}){\bar{1}}_i^{\ast }}$$

$$\pi \left[\left(\sum _{i\in A}{\psi }_i\right)\left(\sum _{j\notin A}{\psi }_j\right)\right]>0.$$

$$N=\sum _{i=1}^L{N}_i$$

$$\underset{N\to \infty }{\lim }{N}_i/N={\kappa }_i>0.$$

$$\bar{F}\stackrel{\mathrm{as}}{\to }F,{\bar{g}}_i^{\ast }\stackrel{\mathrm{as}}{\to }{g}_i^{\ast },\text{ and }{\bar{1}}_i^{\ast }\stackrel{\mathrm{as}}{\to }{1}_i^{\ast }\text{ as }N\to \infty .$$

$$N_j=\sum _{i=1}^L{N}_i{\bar{F}}_{ij}(u),\sum _{i=1}^L{u}_i{N}_i=1$$

$$w_j=\sum _{i=1}^L{w}_i{\bar{F}}_{ij}(u),\sum _{i=1}^L{w}_i=1$$

$$z_i^{m+1}=\frac{u_i{w}_i}{{\sum }_{k=1}^L{u}_k{w}_k}.$$

$${\bar{G}}_{i:}=({\bar{g}}_i^{\ast },{\bar{1}}_i^{\ast }).$$

$$\sqrt{N_i}\left(\left({\bar{F}}_{i:},{\bar{G}}_{i:}\right)-(F_{i:},{\pi }_i[g^{\ast }],{\pi }_i[1^{\ast }])\right)\stackrel{\mathrm{d}}{\to }\mathrm{N}(0,{\mathrm{\Sigma }}_{\mathrm{i}})$$

$${\Sigma }_i=\left(\begin{array}{cc}{\sigma }^i& {\rho }_i\\ {\rho }_i^{\mathrm{t}}& {\tau }_i\end{array}\right),$$

$$A{A}^{\#}A=A\text{, }{A}^{\#}A{A}^{\#}=A^{\#}\text{, and }A{A}^{\#}=A^{\#}A.$$

$$\Psi =\frac{1}{{\sum }_{i=1}^L\pi [{\psi }_i]}\text{and}\ell :=\Psi {\left(1,-\pi [g]\right)}^{\mathrm{t}}\in {\mathbb{R}}^2.$$

$$\sqrt{N}\left({\pi }_{\text{US}}[g]-\pi [g]\right)\stackrel{\mathrm{d}}{\to }\mathrm{N}(0,{\sigma }_{\text{US}}^2(\mathrm{g})),$$

$$\displaystyle\sigma^{2}_{\text{US}}(g)=\sum_{i=1}^{L}\frac{z_{i}^{2}}{\kappa_{i}}\big{\{}$$

$$h=g^{\ast }-\pi [g]1^{\ast },$$

$$\displaystyle\sigma^{2}_{\text{US}}(g)\leq 2\sum_{i=1}^{L}\frac{1}{\kappa_{i}}\Bigg{\{}$$

$$R_{jk}^i:=\frac{{\sigma }_{jk}^i}{\sqrt{{\mathrm{var}}_{{\pi }_i}({\psi }_j^{\ast })}\sqrt{{\mathrm{var}}_{{\pi }_i}({\psi }_k^{\ast })}}.$$

$$\sum _{\begin{array}{c}j\ne i\\ F_{ij}>0\end{array}}\frac{{\mathrm{var}}_{{\pi }_i}({\psi }_j^{\ast })}{{\mathbf{P}}_i[t_j<t_i{]}^2}$$

Full text

1 Introduction

Markov chain Monte Carlo (MCMC) methods have been widely used with great success throughout statistics, engineering, and the natural sciences. However, when estimating tail probabilities or when sampling from multimodal distributions, accurate MCMC estimates often require a prohibitively large number of samples. In this article, we analyze the Eigenvector Method for Umbrella Sampling (EMUS) [48]. We first proposed EMUS as a method for computing free energies, and we demonstrated that it was useful for treating the multimodality that typically arises in that context. Here, we demonstrate that EMUS is an effective general means of addressing the challenges posed not just by multimodality but also tail events, with potential applications to a broad range of problems in statistics, engineering, and the natural sciences.

EMUS was inspired by Umbrella Sampling [49] and other methods such as the Weighted Histogram Analysis Method (WHAM) [27] and the Multistate Bennett Acceptance Ratio (MBAR) [44] for computing potentials of mean force and free energies in statistical mechanics. 111A potential of mean force is the logarithm of a marginal density. A free energy is the logarithm of a normalization constant. Both quantities play fundamental roles in statistical mechanics, e.g., in theories of rates of chemical reactions. We call these stratified MCMC methods since they each adapt the principle of stratified survey sampling to MCMC simulation. That is, they each draw samples concentrated in a collection of subregions of state space and combine the resulting averages in a consistent manner. Stratified MCMC methods are among the most powerful, most successful, and most widely used tools in molecular simulation. (However, in contrast to our presentation here, they are not typically used in molecular simulation to compute averages of general observables.) WHAM, for example, has been instrumental for treating biomolecular processes ranging from protein folding [7] to conductance by ion channels [3].

While the practical utility of stratification has been established in many applications, the advantages and disadvantages of the method have remained poorly understood; cf. [48]. Motivated by the substantial gap between theory and application of stratified MCMC within statistical mechanics, and also by the general challenges posed by multimodality and tail probabilities, the goal of this paper is to develop a clear theoretical explanation of the advantages of EMUS. Our theory suggests new applications of stratified MCMC (and EMUS in particular) to broad classes of sampling problems arising in statistics and statistical mechanics. For example, very recently EMUS was successfully applied to a parameter estimation problem in cosmology [34].

Summary of Main Results

Our most general results are a central limit theorem (CLT) for the EMUS method and a convenient upper bound on the asymptotic variance, cf. Theorem 3.1 and Theorem 3.4. We note that the proof of the upper bound relies on a new class of perturbation estimates for Markov chains which we derived in [47]. These estimates are substantially more detailed than previous results [9]. After proving the CLT, we address the dependence of the sampling error on the choice of strata. In particular, for a representative MCMC method, we estimate the asymptotic variances of trajectory averages sampling the biased distributions, cf. Theorem 3.6. Our estimate shows how factors such as the diameters of the strata influence the asymptotic variances.

In Section 4, we apply the general theory developed in Section 3 to case studies involving tail probabilities and multimodality. Our results concern two limits: a small probability limit and a low-temperature limit. In the small probability limit, we consider estimation of probabilities of the form

[TABLE]

For a broad class of random variables $X$, we show that while the cost of computing $p_{M}$ with relative precision by direct MCMC increases exponentially with $M$, the cost by EMUS increases only polynomially; cf. Section 4.2. In the low-temperature limit, a parameter of the target distribution decreases, intensifying the effects of multimodality on the efficiency of MCMC sampling. We show that the cost of computing an average to fixed precision by direct MCMC increases exponentially in this limit, whereas the cost by EMUS increases only polynomially; cf. Section 4.1. We conclude that EMUS may be dramatically more efficient than direct MCMC sampling when the target distribution is multimodal or when the goal is to compute a small tail probability.

To illustrate our theoretical results, we present a tutorial numerical study applying EMUS to a problem in Bayesian statistics in Section 5. In addition to illustrating the theory, our numerical study demonstrates the problems that may occur when EMUS and other similar stratified MCMC methods are used carelessly. It also addresses practical issues such as the choice of strata and the computation of error bars for averages estimated by EMUS.

The results in this article significantly extend and generalize the ideas in [48]. We first proposed the EMUS method with the goal of analyzing and improving umbrella sampling approaches in free energy calculations. Here, our goal is to establish EMUS as a general variance reduction technique, and we present entirely new results, including an upper bound on the asymptotic variance of EMUS (Theorem 3.4), a condition to guide some aspects of the choice of strata (Remark 3.12), a theoretical argument demonstrating the benefits of EMUS for computing tail probabilities (Section 4.2), numerical results applying EMUS to Bayesian inference (Section 5), a method of correcting problems related to poorly chosen strata (Section 5.3), and a greatly improved numerical method for estimating the standard deviations of quantities computed by EMUS (Appendix G). In addition, we give complete justifications of some results that were stated without proof in [48], including Theorem 3.6 concerning the dependence of the sampling error on the choice of strata. Finally, we note that our results concerning multimodal distributions and the low-temperature limit generalize and clarify the results given in [48]; in particular, our Theorem 4.1 covers periodic boundary conditions and stratification in more than one variable.

2 The Eigenvector Method for Umbrella Sampling

In this section, we present the Eigenvector Method for Umbrella Sampling (EMUS), and we prove that it is consistent. A detailed derivation of the estimator can be found in [48]. We also review a related method, iterative EMUS, and we compare iterative EMUS with the MBAR method from statistical mechanics [44].

2.1 The EMUS estimator

The objective of EMUS is to compute the average

[TABLE]

of a function $g$ with respect to a measure $\pi$ defined on a set $\Omega$. In EMUS, instead of sampling directly from $\pi$, we sample from biased distributions analogous to the strata in stratified survey sampling methods. We then weight the samples from the biased distributions to estimate $\pi[g]$.

We assume that the biased distributions take the form

[TABLE]

for some set $\{\psi_{i}\}_{i=1}^{L}$ of non-negative bias functions defined on $\Omega$. We assume that

[TABLE]

We call the support of $\psi_{i}$ the $i$’th stratum to make an analogy between the biased distributions of EMUS and the strata of stratified survey sampling.

The EMUS estimator is based on two observations. First, one may write $\pi[g]$ in terms of weighted sums of averages over the biased distributions. Let $u\in(0,\infty)^{L}$ be arbitrary, and for any function $h:\Omega\rightarrow\mathbb{R}$, define

[TABLE]

By [48, Equation 16], we have

[TABLE]

where $\mathrm{\mathbf{1}}$ denotes the function equal to one everywhere and

[TABLE]

(For the reader’s convenience, we present complete derivations of equations (1) and (3) in Appendix A.) The parameter $u$ above can be thought of as an initial guess of the unknown weight vector $z\in\mathbb{R}^{L}$ in (2). Its explicit presence here will simplify the description of an iterative version of EMUS in Section 2.2. Outside of this section and Section 2.2, the choice of $u$ is absorbed in the definition of the $\psi_{i}$, i.e. $u_{i}=1$.

Second, the weight vector $z$ can be found by solving a certain eigenproblem. Let $w\in\mathbb{R}^{L}$ be defined by $z_{i}=u_{i}w_{i}$. By [48, Equation 17], we have

[TABLE]

We call $F$ the overlap matrix.

We observe that equations (1) and (3) combine to express $\pi[g]$ as a function of averages over the biased distributions, namely $\pi_{i}[g^{\ast}]$, $\pi_{i}[\mathrm{\mathbf{1}}^{\ast}]$, and $F$. To see this, it suffices to recognize that $F$ is stochastic and $w$ is a probability vector. Therefore, whenever $F$ is irreducible, the solution of the eigenproblem is unique by the Perron–Frobenius theorem. We will assume throughout the remainder of this work that $F$ is irreducible, which amounts to requiring some overlap between neighboring strata; see Lemma 2.1. In general, for any irreducible, stochastic matrix $G\in\mathbb{R}^{L\times L}$, we will let $w(G)\in\mathbb{R}^{L}$ denote the unique solution of

[TABLE]

With this notation, by (1) and (3) we have

[TABLE]

In EMUS, we substitute MCMC estimates for the averages over biased distributions on the right hand side of (4) to estimate $\pi[g]$. To be precise, let $X^{i}_{t}$ be a Markov process ergodic for $\pi_{i}$. We call $X^{i}_{t}$ the biased process sampling the biased distribution $\pi_{i}$. The EMUS algorithm proceeds as follows:

For each $i=1,\dots,L$, compute $N_{i}$ steps of the process $X^{i}_{t}$. 2. 2.

Compute the averages

[TABLE] 3. 3.

Compute $w(\bar{F})$ numerically, for example from the QR factorization of $I-\bar{F}$ [19]. 4. 4.

Compute the estimate

[TABLE]

of $\pi[g]$.

Recall that $w(\bar{F})$ is defined only if $\bar{F}$ is irreducible. In the following lemma, whose proof is in Appendix B, we state simple criteria for the irreducibility of $F$ and $\bar{F}$:

Lemma 2.1**.**

The overlap matrix $F$ is irreducible if and only if for every $A\subset\{1,2,\dots,L\}$, we have

[TABLE]

*The approximate overlap matrix $\bar{F}$ is irreducible if and only if for every $A\subset\{1,2,\dots,L\}$, the set $\cup_{i\in A}\{x:\psi_{i}(x)>0\}$ contains at least one sample point generated from one of the biased processes $X^{j}_{t}$ with $j\notin A$. *

We claim that the EMUS estimator is consistent; that is, $\pi_{\text{US}}[g]$ converges almost surely to $\pi[g]$ as the total number of samples tends to infinity. To make this precise, we require the following assumption on the growth of $N_{i}$ with the total number of samples: {assumption} Let

[TABLE]

be the total number of samples from all biased distributions. Assume that for each $i$,

[TABLE]

That is, assume that when $N$ is large, the proportion of samples drawn from the $i$’th biased distribution is fixed and greater than zero.

We now prove that EMUS is consistent:

Lemma 2.2**.**

*Under Assumption 2.1 and the irreducibility condition (6), $\pi_{\text{US}}[g]$ converges almost surely to $\pi[g]$ as the total number of samples $N$ tends to infinity. *

Proof 2.3**.**

Since the processes $X^{i}_{t}$ are ergodic,

[TABLE]

*Moreover, by Lemma D.1 in Appendix D, $w(G)$ is continuous at $F$. (Technically, $w(G)$ admits an extension to the set of all $L\times L$ matrices, which is continuous at $F$.) Therefore, as a function of $\bar{F}$, $\bar{g}^{\ast}_{i}$, and $\bar{\mathrm{\mathbf{1}}}^{\ast}_{i}$, $\pi_{\text{US}}[g]$ is continuous at $F$, $\pi_{i}[g^{\ast}]$, and $\pi_{i}[\mathrm{\mathbf{1}}^{\ast}]$. It follows by the continuous mapping theorem and equation (4) that ${\pi_{\text{US}}[g]\xrightarrow[]{\rm{as}}\pi[g]}$. *

2.2 Iterative EMUS and Vardi’s Estimator for Selection Bias Models

EMUS resembles certain methods for computing normalization constants of families of probability densities [17, 51, 35, 26]. The resemblance arises because the weights in the third step of EMUS are the normalization constants of the biased distributions. These methods have been used, for example, to compute Bayes factors in model selection problems [17] and for computations related to selection bias models [51]. In this section, we explain how EMUS relates to Vardi’s estimator for selection bias models [51] and its descendants such as the popular Multistate Bennett Acceptance Ratio (MBAR) method [44]. In addition to comparing EMUS with these methods, we review a method, iterative EMUS, for solving the nonlinear system of equations defining Vardi’s estimator [48]. The first iterate of this method is exactly the EMUS estimator described in Section 2.1. We note that our analysis and calculations in Sections 3, 4, and 5 pertain only to the EMUS method and not iterative EMUS or MBAR. However, the similarity between the methods suggests that our results may generalize.

Given the notation developed in Section 2.1, we can express Vardi’s estimate $z^{\rm V}$ of the weight vector as the vector with entries $z^{\rm V}_{i}=u_{i}\,N_{i}$ where $u$ solves the equation

[TABLE]

and where we have now made the dependence of the matrix $\bar{F}$ in (5) on the choice of $u\in\mathbb{R}^{L}$ explicit. By [51, Theorem 1], this nonlinear equation determines $z^{V}$ uniquely whenever the irreducibility criterion of Lemma 2.1 holds.

Vardi’s estimator was originally derived assuming that the samples $X^{i}_{t}$ from the biased distributions were i.i.d. In that case, it is the nonparametric maximum likelihood estimator of the target distribution $\pi$ given samples from the biased distributions $\pi_{i}$ [51], and it has certain optimality properties [18]. Several adjustments to the estimator have been proposed for the case of samples from Markov processes. In the Multistate Bennett Acceptance Ratio (MBAR) method, one replaces the factors $N_{i}$ appearing in the summand in (5) with effective sample sizes $n_{i}$, which are computed from estimates of the integrated autocovariance of a family of functions [44]. (In some versions of MBAR, the sample average over all $N_{i}$ points is replaced with a sample average over the $n_{i}$ points obtained by including only every $N_{i}/n_{i}$’th point along the trajectory $X^{i}_{t}$.) Another recent work proposes different effective sample sizes computed by minimizing an estimate of the asymptotic variance of the estimator [12]. In fact, the estimator is consistent with $N_{i}$ replaced by any fixed positive number [12]. We have found that our numerical results do not depend sensitively on the choice of effective sample size, so we use $N_{i}$ for simplicity.

We now review iterative EMUS, which we introduced in [48]. Iterative EMUS may be understood as a fixed point iteration for solving equation (8). The iteration proceeds as follows:

As an initial guess for $z^{V}$, choose a positive vector $z^{0}\in\mathbb{R}^{L}$. Set $m=0$. Choose a tolerance $\tau>0$. 2. 2.

Compute $\bar{F}_{ij}(u)$ for $u_{i}=z^{m}_{i}/N_{i}$. Solve the eigenvector equation

[TABLE]

for $w\in\mathbb{R}^{L}$ to obtain the updated estimate $z^{m+1}$ of $z^{\rm V}$ with entries

[TABLE] 3. 3.

If $\max_{i}\left\lvert w_{i}-\frac{N_{i}}{N}\right\rvert>\tau$, then increment $m$ and repeat step 2.

In [48] we show that the eigenvector equation (9) has a unique solution for every $m$, and we suggest a numerical method for finding the solution. We also discuss the convergence of iterative EMUS, and we show that for every fixed $m$, $z^{m}$ is a consistent estimator of the weight vector $z$. If one chooses $z^{0}_{i}=N_{i}/N$, then $z^{1}$ is the EMUS estimate of $z$, $w(\bar{F})$.

2.3 Related MCMC methods

EMUS belongs to a large class of MCMC methods that by various mechanisms promote a more uniform sampling of space. For example, in parallel tempering [46, 16], one uses MCMC samples drawn from a distribution or sequence of distributions close to the uniform distribution to speed sampling of the target distribution. The bias introduced by the choice of distributions is corrected either by reweighting the samples or by a replica exchange strategy [16]. The Wang–Landau [52] and Metadynamics methods [28] adaptively construct a biased distribution to achieve uniform sampling in certain coordinates. The temperature accelerated molecular dynamics method [33] is also designed to achieve uniform sampling in a given coordinate, but it works by entirely different means. In EMUS and other stratified MCMC methods, one achieves more uniform sampling by ensuring that each stratum contains points from at least one MCMC simulation. EMUS is perhaps most similar in spirit to the parallel Markov chain Monte Carlo method [50].

3 Error Analysis of EMUS

Here, we develop tools for analyzing the error of EMUS. First, in Section 3.1, we prove a CLT for EMUS, and we derive a convenient upper bound on the asymptotic variance. Then, in Section 3.2, we analyze the dependence of the asymptotic variance of EMUS on the choice of biased distributions. We use these tools in Section 4 to prove limiting results demonstrating the advantages of EMUS for multimodal distributions and tail probabilities. In addition, our CLT for EMUS is the basis for both practical error estimates (Section 5 and Appendix G) and also a method of optimizing the allocation $\kappa_{i}$ of the samples among the different biased distributions [48].

3.1 A CLT for EMUS and an Estimate of the Asymptotic Variance

In this section, we prove a Central Limit Theorem (CLT) for EMUS, and we derive an upper bound on the asymptotic variance $\sigma^{2}_{\text{US}}(g)$ of $\pi_{\text{US}}[g]$. To prove the CLT for EMUS, we must assume that a CLT holds for trajectory averages over the biased processes: {assumption} For any matrix $H$, let $H_{i:}$ denote the $i$’th row of $H$. Define $\bar{G}\in\mathbb{R}^{L\times 2}$ by

[TABLE]

Assume that

[TABLE]

for some asymptotic covariance matrix $\Sigma_{i}\in\mathbb{R}^{(L+2)\times(L+2)}$ of the form

[TABLE]

where $\sigma^{i}\in\mathbb{R}^{L\times L}$ denotes the asymptotic covariance of $\bar{F}_{i:}$ with itself, $\rho_{i}\in\mathbb{R}^{L\times 2}$ denotes the asymptotic covariance of $\bar{F}_{i:}$ with $\bar{G}_{i:}$, and $\tau_{i}\in\mathbb{R}^{2\times 2}$ denotes the asymptotic covariance of $\bar{G}_{i:}$ with itself.

We expect a CLT to hold for most MCMC methods, target distributions, and target functions of interest in statistics and statistical mechanics. We refer to [41] for a comprehensive review of conditions guaranteeing a CLT. In Theorem 3.6 of Section 3.2, we prove a CLT and an estimate of the asymptotic variance for a simple family of processes which one might use to sample the biased distributions in an application of EMUS.

We now prove a CLT for $\pi_{\text{US}}[g]$, and we give a formula expressing the asymptotic variance $\sigma^{2}_{\text{US}}(g)$ of $\pi_{\text{US}}[g]$ in terms of the asymptotic variances $\Sigma_{i}$ of the trajectory averages. In this formula, $(I-F)^{\#}$ denotes the group generalized inverse of $I-F$; the group inverse $A^{\#}$ of a matrix $A$ is characterized by the properties

[TABLE]

We refer to [19] for a detailed explanation of the properties of the group inverse, a proof that $(I-F)^{\#}$ exists whenever $F$ is stochastic and irreducible, and an algorithm for computing $(I-F)^{\#}$.

In Theorem 3.1 and below we impose the following assumption: {assumption} The processes $X^{i}_{t}$ sampling the biased distributions are independent.

This assumption does not hold for all stratified MCMC methods. For example, in replica exchange umbrella sampling one periodically allows configuration exchanges between neighboring processes; see [32] for a general discussion of replica exchange strategies and [45] for an application of replica exchange in a method similar to EMUS. The result is a single process taking values in $\mathbb{R}^{L\times d}$ and sampling the product distribution $\Pi(x_{1},x_{2},\dots,x_{L})=\pi_{1}(x_{1})\pi_{2}(x_{2})\dots\pi_{1}(x_{L})$. In this case, a CLT would still hold for EMUS, but the asymptotic variance would take a different form. We also assume that $u_{i}=1$ for all $i$. As already mentioned, we can equivalently absorb the choice of $u_{i}$ into the choice of the bias function $\psi_{i}$.

Theorem 3.1**.**

Let Assumptions 2.1, 3.1, and 3.1 and the irreducibility condition (6) hold. Let $g$ be square integrable over $\pi$, so $\pi[g^{2}]<\infty$. Defining

[TABLE]

Let $\mathfrak{g}\in\mathbb{R}^{L}$ be the vector with $\mathfrak{g}_{i}:=\ell\cdot(\pi_{i}[g^{\ast}],\pi_{i}[\mathrm{\mathbf{1}}^{\ast}])$. We have

[TABLE]

where

[TABLE]

Proof 3.2**.**

*The result follows using the delta method and a formula expressing $w^{\prime}(F)$ in terms of $(I-F)^{\#}$; we give the details in Appendix D. *

We now derive a convenient upper bound on the asymptotic variance $\sigma^{2}_{\text{US}}(g)$. In Section 4, we use this bound to analyze the efficiency of EMUS in the low-temperature limit and in the limit of small tail probabilities. Our bound is based on the probability $\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]$ defined below:

Definition 3.3**.**

*Let $Y_{n}$ be the Markov chain with state space $\{1,2,\dots,L\}$ and transition matrix $F$. Let ${\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]}$ denote the probability that $Y_{n}$ hits $j$ before returning to $i$, conditioned on $Y_{0}=i$. *

Theorem 3.4**.**

Let Assumptions 2.1, 3.1, and 3.1 and the irreducibility condition (6) hold. Let $g$ be square integrable over $\pi$, so $\pi[g^{2}]<\infty$. Let $\sigma^{2}(g)$ be the asymptotic variance of $\pi_{\text{US}}[g]$, and for any measure $\nu$ and function $f$ let $\operatorname{var}_{\nu}(f)$ be the variance of $f$ over $\nu$. Define the function

[TABLE]

and let $\operatorname{C}(\bar{h}_{i})$ be the asymptotic variance of the trajectory average of $h$ over the biased process $X^{i}_{t}$. We have

[TABLE]

where $R^{i}\in\mathbb{R}^{L\times L}$ with

[TABLE]

Proof 3.5**.**

*The result follows from Theorem 3.1, using the perturbation bounds which we derived in [47]. Details appear in Appendix D. *

3.2 Dependence of the Asymptotic Variance on the Choice of Strata

In this section, we consider how the choice of strata influences the factors in the upper bound (14) on $\sigma^{2}_{\text{US}}(g)$. Roughly, the asymptotic variances $C(\bar{h}_{i})$ and $\operatorname*{tr}(R^{i})$ characterize the sampling error, and for each $i$ the factor

[TABLE]

measures the sensitivity of the EMUS estimator to sampling errors associated with $\pi_{i}$.

We show in Section 3.2.1 that the factors $C(\bar{h}_{i})$ and $\operatorname*{tr}(R^{i})$ characterizing the sampling error may be controlled by decreasing the diameters of the strata. We show in Section 3.2.2 that that $\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]$ may be controlled by ensuring sufficient overlap between neighboring strata. This last observation leads to a practical condition guiding the choice of strata; see Remark 3.12 and (38).

Our theorems in this section apply to a specific class of strata and Markov processes that are broadly representative of those employed in practical applications. Thus, the assumptions made here are much stronger than those made in proving the CLT for EMUS, for example. We discuss how our results might extend to more general implementations of stratified MCMC after the statement of Theorem 3.6 in Section 3.2.1 and in Remark 3.12.

3.2.1 Asymptotic Variances of MCMC Averages

Here, we consider the effect of the choice of strata on the asymptotic variances $C(\bar{h}_{i})$ and $\operatorname*{tr}(R^{i})$. Because such a diverse variety of biased processes and distributions could in principle be used, it is futile in our opinion to try for a completely general result. Instead, motivated by the efficiency analysis undertaken in Section 4, we introduce a simple parametric family of bias functions, and for this family we state Assumption 3.2.1 relating the diameters of the strata with the asymptotic variances. In Theorem 3.6, we verify Assumption 3.2.1 for one representative class of biased processes. Finally, at the end of this section, we explain why we expect the assumption to hold for other choices of biased processes and distributions.

Consider the following representative class of bias functions: Given a family of sets ${\{U_{i}:i=1,\dots,L\}}$ with ${\cup_{i=1}^{L}U_{i}=\Omega}$, define

[TABLE]

where $\mathrm{\mathbf{1}}_{U_{i}}$ denotes the characteristic function of $U_{i}$. Assume that the sets $U_{i}$ are chosen so that the irreducibility criterion of Lemma 2.1 holds. For example, suppose that $\Omega=[0,1]^{d}$ is the $d$-dimensional unit cube. One might choose $K\in\mathbb{N}$, set $h:=1/K$, and define

[TABLE]

covering $\Omega$ uniformly by a grid of strata having diameters proportional to $h$. We use this uniform grid as a device when analyzing the effect of the stratum size on the efficiency of EMUS in Section 4. However, while such a naïve choice may suffice for small $d$, it is not practical for large $d$. We discuss appropriate bias functions for high-dimensional problems later in this section and again in Section 5.1.

Since we wish to study grids like (17) as $h=1/K$ varies, we state our assumption on asymptotic variances in terms of the following parametric family of strata: Let $x_{0}\in\Omega$, and let $Z\subset\mathbb{R}^{d}$ be a bounded set containing [math]. For each $h>0$, define a stratum and a biased distribution by

[TABLE]

To make the connection with (17), one should imagine that $Z=[-1,1]^{d}$ and that $x_{0}$ is one of the grid points $h\mathbf{i}$.

Assumption 3.2.1 characterizes the dependence of the asymptotic variance of MCMC averages over $\pi_{h}$ on the parameter $h$:

{assumption}

Assume that $f:\Omega\rightarrow\mathbb{R}$ has finite variance $\operatorname{var}_{h}(f)$ over $\pi_{h}$, and define $\sigma^{2}_{h}(f)$ to be the asymptotic variance of an MCMC trajectory average approximating $\pi_{h}[f]$. Write

[TABLE]

for some potential $V:\Omega\rightarrow\mathbb{R}$ and inverse temperature $\beta>0$. We assume

[TABLE]

for some $C,a,b\geq 0$ independent of $h$,$Z$, and $f$.

To motivate Assumption 3.2.1, we prove that a special case holds for a representative class of processes sampling the biased distributions, cf. Theorem 3.6. Assume that the potential $V$ appearing in the assumption is continuously differentiable. Let $Z\subset\mathbb{R}^{d}$ be either a convex polyhedron or a set with $C^{3}$ boundary.222The boundary of a set is $C^{3}$ if in a neighborhood of each point on the boundary, the boundary is the graph of a three times continuously differentiable function. Now let $X^{h}_{t}$ be the overdamped Langevin process with reflecting boundary conditions on $Z_{h}$. This process is defined by the Fokker–Planck equation

[TABLE]

where $\beta$ and $V$ are the inverse temperature and potential defined in Assumption 3.2.1 and $\mathbf{n}(x)$ denotes the inward unit normal to $\partial Z_{h}$ at $x$. That is, $X^{h}_{t}$ is the unique Markov process so that if $X^{h}_{0}$ has density $p(x)$, then $X^{h}_{t}$ has density $u(x,t)$. The existence of the reflected process is established in [53, 2] when $Z$ is a convex polyhedron and in [13, Chapter 8] when $Z$ has $C^{3}$ boundary. A simple introduction to the reflected process and its properties appears in [36, Chapter 4]. We show in Theorem 3.6 that $X^{h}_{t}$ is ergodic for $\pi_{h}$, at least when $Z$ is bounded. The reflected process $X^{h}_{t}$ shares many features with the processes used in practical stratified MCMC methods. In particular, it is closely related to the (unreflected) overdamped Langevin process $Y_{t}$ [43, 40].

We now verify Assumption 3.2.1 for the reflected process:

Theorem 3.6**.**

Assume that $f:\Omega\rightarrow\mathbb{R}$ has finite variance $\operatorname{var}_{h}(f)$ over $\pi_{h}$. Let $Z$ either have $C^{3}$ boundary or be convex. Assume that $V$ is continuously differentiable. Suppose that $X^{h}_{t}$ is stationary; that is, $X^{h}_{0}$ has distribution $\pi_{h}$. Let $\bar{f}^{h}:=\frac{1}{T}\int_{t=0}^{T}f(X^{h}_{t})\,dt$ be the continuous time trajectory average of $f$. We have

[TABLE]

where

[TABLE]

*The constant $\Lambda$ depends only on $Z$, not on $h$, $\beta$, $V$, or $f$. *

Proof 3.7**.**

*See Appendix E. *

There are three major differences between the choice of strata and sampling scheme specified in this subsection and those typical of practical applications: First, in molecular simulations, one typically chooses Gaussian bias functions instead of piecewise constant. Second, practical methods must be discrete in time, e.g., one might use a discretization of the continuous time process $X^{h}_{t}$. Third, for high-dimensional problems, one typically stratifies only a certain low-dimensional reaction coordinate or collective variable.

In the first case, for Gaussian bias functions, a version of Theorem 3.6 holds with minor adjustments; we omit the exact statement and proof for simplicity. In the second case, for discretizations of Langevin dynamics, the asymptotic variances of trajectory averages are closely related to the corresponding averages for the continuous time dynamics: In fact, under some conditions on the potential $V$,

[TABLE]

where $C_{\Delta t}(f)$ is the asymptotic variance of the trajectory average of $f$ for the discretization with time step $\Delta t$ and $C(f)$ is the asymptotic variance for the continuous time process  [31, Section 3.2]. For other discrete time processes, we expect Assumption 3.2.1 to hold with different exponents $a$ and $b$. For example, the affine invariance property of the affine invariant ensemble sampler [20] suggests $a=0$.

The third case is subtle. When $d$ is large, one typically stratifies only in a function ${\theta:\Omega\subset\mathbb{R}^{d}\rightarrow\mathbb{R}^{\ell}}$ with $\ell$ much smaller than $d$. To be precise, one might choose a uniform grid of nonnegative functions $\eta_{i}:\mathbb{R}^{\ell}\rightarrow\mathbb{R}$ defined as in (17), but with supports covering $\theta(\Omega)\subset\mathbb{R}^{\ell}$ instead of $\Omega\subset\mathbb{R}^{d}$. One would then define the bias functions

[TABLE]

(We make a similar choice in our calculations in Section 5, cf. the natural stratification (31).) For a clever choice of $\theta$, these biased distribution may be much easier to sample than the target distribution. For example, suppose that the marginal $\pi_{\theta}$ of $\pi$ in $\theta$ were multimodal, but that the conditional distributions $\pi(\cdot\mid\theta=\theta_{0})$ were unimodal or otherwise easy to sample for each fixed $\theta_{0}$. In that case, for $h$ sufficiently small, each biased distribution would be unimodal, hence easy to sample. (Recall that $h$ sets the diameters of the strata for the grid of bias functions defined in (17), so $h$ small means that the diameter of the support of $\eta_{i}$ is small.) In free energy calculations, it is often possible to choose such a $\theta$ based on intuition or scientific principles; see [48, 30] for discussion. When computing tails or marginals, the problem itself typically suggests a particular $\theta$; cf. the natural stratification in Section 5.1.

The reader will notice that bias functions of the form (22) will typically have infinite support, rendering the bound in Assumption 3.2.1 useless. In this case, one might hope for a similar bound with the potential function $V$ replaced by the free energy

[TABLE]

where $\pi_{\theta}$ is the marginal density of $\pi$ in $\theta$. Roughly, this replacement will be valid when, under the dynamics of the MCMC processes sampling $\pi$, the distribution of any variable (any function of the process) converges rapidly to its conditional distribution under $\pi$ given the current value of the $\theta$ variable. This will occur, for example, when the marginal in $\theta$ is multimodal or otherwise difficult to sample, but the conditional distributions are easy to sample. In general an effective choice of $\theta$ will be one for which conditional equilibration given $\theta$ occurs much more rapidly than the overall time to convergence of the process. More on the effective dynamics of low-dimensional variables can be found in [37] or [29].

3.2.2 Controlling the Probabilities $\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]$

Here, we examine the effect of the choice of strata on the factors in display (15) that appear in our upper bound (14) on $\sigma^{2}_{\text{US}}(g)$.

We begin with a lemma estimating ${\operatorname{var}_{\pi_{i}}(\psi_{j}^{\ast})}/{\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]^{2}}$ in terms of $F_{ij}$:

Lemma 3.8**.**

We have

[TABLE]

Proof 3.9**.**

We have

[TABLE]

where $X_{t}$ denotes the Markov chain with transition matrix $F$. Therefore, since $\psi_{j}^{\ast}(x)\in[0,1]$,

[TABLE]

We now estimate the size of $F_{ij}$ for piecewise constant bias functions such as the uniform grid (17):

Lemma 3.10**.**

Assume as in (16) that the bias functions are piecewise constant, and write ${\pi(x)\propto\exp(-\beta V(x))}$. We have

[TABLE]

In particular, for the uniform grid of strata (17), we have

[TABLE]

*for any $\mathbf{i},\mathbf{j}\in\mathbb{Z}^{d}$ so that $F_{\mathbf{i}\mathbf{j}}>0$. *

Proof 3.11**.**

We have

[TABLE]

which proves the first claim made in the statement of the lemma.

Now, for the uniform grid of strata (17), the minimum nonzero value of $\lvert U_{\mathbf{i}}\cap U_{\mathbf{j}}\rvert/\lvert U_{\mathbf{i}}\rvert$ is $1/2^{d}$, attained when $\mathbf{j}=(1,1,\dots,1)+\mathbf{i}$. Moreover, except for a set of measure zero, each $x\in\mathbb{R}^{d}$ lies within $2^{d}$ strata, so $\left\lVert\sum_{k=1}^{L}\mathrm{\mathbf{1}}_{U_{k}}\right\rVert_{\infty}=2^{d}$. Finally, we have

[TABLE]

*and the result follows. *

Remark 3.12** (A Condition to Guide the Choice of Strata).**

*Lemmas 3.8 and 3.10 suggest a practical constraint on the choice of strata: To ensure that the calculation of the weights is not too sensitive to sampling errors, it will suffice to choose strata so that nonzero entries of $F$ are not too small. We let this condition guide the choice of strata in Section 5, cf. (38). However, the condition is only sufficient, not necessary. For example, consider a uniform grid of Gaussian bias functions similar to (17), but with Gaussian densities having mean $\mu=h[-1,1]^{d}+h\mathbf{i}$ and variance $\sigma^{2}=h^{2}$ replacing the characteristic functions $\mathrm{\mathbf{1}}_{U_{\mathbf{i}}}$. In that case, even though $F$ will be dense and may have some extremely small nonzero entries, one can still control (15) by decreasing $h$, under some conditions on $\pi$. We omit the exact statement and proof for simplicity. *

Despite the exponential dependence on $d$ in (23), EMUS and other stratified MCMC methods are advantageous for high-dimensional problems because it often suffices to stratify only a low-dimensional collective variable. In such cases, the dimension of the grid of strata is much smaller than dimension of the state space $\Omega$; see our discussion of collective variables in Section 3.2.1 and our computations in Section 5. It is important to keep this in mind when reading our results below. Also, one may define a uniform grid of strata so that (15) increases only as $d^{2}$ with dimension, not exponentially. We construct such a grid in Appendix C.

4 Limiting Results as a Rationale for EMUS

In this section, we analyze the efficiency of EMUS in two limits: First, we consider a low temperature limit, where we write $\pi(x)\propto\exp(-\beta V(X))$ and let the inverse temperature $\beta$ increase, concentrating the target distribution at its modes and intensifying the effects of multimodality on the efficiency of MCMC sampling. Second, we consider the estimation of increasingly small tail probabilities. Our goal in each case is to elucidate the advantages and disadvantages of EMUS for a broad class of problems, providing a rationale for the use of the method. We hope that others will use the tools of Section 3 in similar fashion to develop their own novel applications of EMUS.

4.1 Limit of Low Temperature

Let the target distribution take the form

[TABLE]

for some potential $V$ and inverse temperature $\beta>0$, as in Section 3.2. In this section, we analyze the efficiency of EMUS in the low temperature limit as $\beta$ tends to infinity with $V$ fixed. We observe that $\pi_{\beta}$ concentrates at its modes (the minima of $V$) in this limit. As a consequence, MCMC methods for sampling $\pi_{\beta}$ undergo transitions between modes only rarely, which makes direct MCMC sampling increasingly inefficient. To be precise, we show that the asymptotic variance of a trajectory average of the overdamped Langevin dynamics increases exponentially with $\beta$ in the worst case. On the other hand, we show that the asymptotic variance of the EMUS estimate of the same average increases only polynomially. Therefore, EMUS is dramatically more efficient than direct sampling in the low temperature limit.

We consider the low temperature limit because it provides a convenient sequence of increasingly difficult to sample multimodal distributions: By analyzing EMUS in the low temperature limit, we hope to clarify its advantages for multimodal problems in general. We have no other interest in low temperature.

We now examine the overdamped Langevin dynamics

[TABLE]

in the low temperature limit. (The overdamped Langevin dynamics is ergodic for $\pi_{\beta}$ under certain conditions on $V$; see [42] for example.) For typical potentials $V$, the generator

[TABLE]

of (24) has a spectral gap that shrinks exponentially with $\beta$; that is, for some $c>0$,

[TABLE]

where $\lambda_{1}$ is the greatest nonzero eigenvalue of $\mathscr{L}$. We refer to [31, Section 2.5] for a review of results on the spectrum of $\mathscr{L}$, and we refer to [21] for precise conditions on $V$ which guarantee (25). Now let $v_{1}$ be an eigenfunction corresponding to $\lambda_{1}$ normalized so that $\pi_{h}[v_{1}^{2}]=1$. By formula (52), the asymptotic variance $\sigma^{2}_{\beta}(v_{1})$ of the trajectory average of $v_{1}$ satisfies

[TABLE]

indicating that the cost of estimating $\pi[v_{1}]$ by direct MCMC grows exponentially with $\beta$.

Having analyzed the overdamped Langevin dynamics, we now examine EMUS in the low temperature limit. For convenience, we assume that $\Omega$ is the unit cube $[0,1]^{d}\subset\mathbb{R}^{d}$ with periodic boundary conditions; to be more precise, we let ${\Omega=\mathbb{R}^{d}/\mathbb{Z}^{d}}$ be the set of all points in $\mathbb{R}^{d}$ with $x$ and $y$ identified if and only if $x-y\in\mathbb{Z}^{d}$. Periodic boundary conditions are typical of problems in chemistry and computational statistical mechanics. We do not see any difficulties in generalizing our results to other types of domains.

As $\beta$ increases, we must make the supports of the bias functions smaller. We accomplish this by adjusting the parameter $h$ in a uniform grid of bias functions similar to those defined in (17). To be precise, we fix $K\in\mathbb{N}$, set $h:=1/K$, and define

[TABLE]

This family of $K^{d}$ bias functions is a partition of unity over $\Omega$, and the support of the $\mathbf{i}$’th bias function is

[TABLE]

For convenience, we treat the index $\mathbf{i}$ as an element of $\mathbb{Z}^{d}/K\mathbb{Z}^{d}$; that is, we let $\mathbf{i}$ be periodic with period $K$ in each of its components, identifying $(0,i_{2},\dots,i_{d})$ with $(K,i_{2},\dots,i_{d})$, for example. Figure 1 illustrates such a family of bias functions, and it demonstrates the appropriate relationship between $\beta$ and $h$.

We now show that the asymptotic variance of EMUS increases at most polynomially with $\beta$ when $K$ is chosen appropriately. In light of the above discussion, this means that EMUS may be dramatically more efficient than direct sampling for multimodal problems. We note that despite the exponential dependence on $d$ in (27) below, EMUS and other stratified MCMC methods are often advantageous for high-dimensional multimodal problems; see our discussion of low-dimensional collective variables in Section 3.2 and also our computations in Section 5.

Theorem 4.1**.**

For any bounded continuous function $g$, let $\sigma^{2}_{\beta,\text{US}}(g)$ denote the asymptotic variance of $\pi_{\beta,\text{US}}[g]$. Let the bias functions be defined by (26) with $K$ equal to the least integer greater than $\beta$; that is,

[TABLE]

Take $\kappa_{\mathbf{i}}=1/K^{d}$. Let Assumption 3.2.1 hold. We have

[TABLE]

*for constants $C,q>0$ independent of $g$ and $\beta$, but depending on $V$ and the constants in Assumption 3.2.1. *

Proof 4.2**.**

*The proof is a straightforward application of the theory developed in Section 3; we present the details in Appendix F. *

Our proof of Theorem 4.1 relies on the perturbation bounds which we derived in [47]. These bounds allow one to estimate the sensitivity of $w(F)$ to small perturbations of $F$. Most perturbation bounds in the literature predict that $w(F)$ is highly sensitive when the spectral gap of $F$ is small, but ours show that this is not always the case. (The spectral gap is $1-\lvert\lambda_{2}\rvert$, where $\lambda_{2}$ is the eigenvalue of $F$ with second largest absolute value.) In the low-temperature limit, the spectral gap of $F$ decreases exponentially with $\beta$; see [47] for a simple example of this phenomenon. Nonetheless, using our bounds, we show that the cost to compute averages by EMUS increases only polynomially in $\beta$.

4.2 Limit of Small Probability

In this section, we assess the performance of EMUS for computing tail probabilities. To be precise, we let $\Omega=[0,\infty)$, and we consider estimation of probabilities of the form

[TABLE]

We show that for a broad class of distributions $\pi$, the cost of computing $p_{M}$ with relative precision by direct MCMC increases exponentially with $M$, whereas the cost by EMUS increases only polynomially. Thus, EMUS is dramatically more efficient than direct sampling for computing the probabilities of tail events.

In Assumption 4.2 below, we state the conditions which we will impose on $\pi$ in our analysis. These conditions specify a simple class of problems for which strong conclusions may be drawn. Similar results hold more generally. For example, in Section 5, we report the results of a computational experiment demonstrating the advantages of EMUS for computing tails of a marginal density.

{assumption}

Write

[TABLE]

for some potential function $V:[0,\infty)\rightarrow\mathbb{R}$. Assume that for some $M_{0}\geq 0$:

Whenever $x\geq M_{0}$,

[TABLE] 2. 2.

For some $\alpha\in(0,1)$ and $c>0$, whenever $x\geq M_{0}$,

[TABLE]

For example, we might have

[TABLE]

Condition (29) in Assumption 4.2 implies geometric ergodicity of the overdamped Langevin dynamics with potential $V$ [41]. We rely on this fact to motivate Assumption 4.2 concerning the convergence of MCMC processes sampling biased distributions with unbounded support. Interestingly, we use the same condition to prove lower bounds on some of the entries of the overlap matrix; cf. Lemma F.3.

Condition (28) in Assumption 4.2 implies

[TABLE]

whenever $M\geq M_{0}$ for some $D,\gamma>0$. Therefore, the relative variance $\rho_{M}^{2}$ of $\mathrm{\mathbf{1}}_{[M,\infty)}$ over $\pi$ satisfies

[TABLE]

We conclude that estimating $p_{M}$ with relative accuracy by a direct MCMC method (or even Monte Carlo with independent samples) requires a number of samples increasing exponentially with $M$.

By contrast, we show that for an appropriate choice of bias functions, the cost to estimate $p_{M}$ by EMUS increases only polynomially in $M$. For each $M>0$ and $K\in\mathbb{N}$, let

[TABLE]

and define the family of $K+2$ bias functions

[TABLE]

As in Section 4.1, let $U_{i}$ denote the support of $\psi_{i}$. This family of bias functions is a partition of unity on $[0,\infty)$; see Figure 2.

We now address the cost of estimating $p_{M}$ by EMUS. First, we observe that Assumption 3.2.1 on the asymptotic variances of MCMC averages does not cover the sampling of $\pi_{K}$ and $\pi_{K+1}$, since the supports of these distributions are unbounded. Thus, we require the following assumption.

{assumption}

Let $f:[0,\infty)\rightarrow\mathbb{R}$, and define $\sigma^{2}_{i}(f)$ to be the asymptotic variance of an MCMC trajectory average approximating $\pi_{i}[f]$ for $i=K,K+1$. We assume

[TABLE]

for some $D$ independent of $M$ and $f$.

In fact, since Assumption 4.2 implies that the overdamped Langevin dynamics is ergodic for $\pi(x)=\exp(-V(x))$ on the unbounded domain $\Omega=\mathbb{R}$, we fully expect (but do not prove here) that under Assumption 4.2, Assumption 4.2 holds for overdamped Langevin constrained (by reflection as in (19)) to remain in the support of $\pi_{K}$ or $\pi_{K+1}.$ Alternatively the reader may simply assume that we draw i.i.d. samples from the biased distributions. All our results hold in that case.

We show in Theorem 4.3 that the relative asymptotic variance of the EMUS estimate of $p_{M}$ grows only polynomially with $M$ for a broad class of target distributions $\pi$. Therefore, EMUS may be dramatically more efficient than direct MCMC sampling when the goal is to compute tail probabilities. We observe that while the hypotheses of the theorem are somewhat restrictive, similar results hold more generally; for example, see Section 5 where we compute tails of a marginal density.

Theorem 4.3**.**

Let Assumptions 3.2.1, 4.2, and 4.2 hold. Set

[TABLE]

Define a family of $K+2$ bias functions $\psi_{i}$ by (30). Take $\kappa_{i}=1/(K+2)$. Let $\sigma^{2}_{M,\text{US}}$ denote the asymptotic variance of the EMUS estimate of $p_{M}$. We have

[TABLE]

for some constant $C>0$ depending on $V$ but not on $M$.

For example, suppose that

[TABLE]

where $\tilde{V}$ has bounded support and $r\geq 1$. Then $\lvert V^{\prime}(x)\rvert\leq C(1+M^{r-1})$, and so

[TABLE]

Proof 4.4**.**

*The proof is similar to the low temperature limit, Theorem 4.1, but with complications arising because not all strata are bounded and because here we consider the relative variance instead of the variance; see Appendix F. In particular, we require Assumption 4.2 to show that one can in fact choose $h$ so that all nonzero entries of $F$ are bounded above zero uniformly as $M$ increases; cf. Lemma F.3. This is the only part of the proof relying on Assumption 4.2. *

5 EMUS for tails: An example from Bayesian inference

We demonstrate the use of EMUS for efficiently exploring and visualizing distributions. In particular, we show how EMUS may be used to efficiently compute both marginal densities and also tail probabilities of the form $\mathrm{\mathbf{P}}[\eta(Z)\geq\varepsilon^{-1}]$ where $\eta(Z)$ is a real valued function of a high-dimensional random variable $Z$. For both tails and marginals, there is a natural and easy to implement choice of strata, which we describe in Section 5.1.

In Section 5.3, we calculate two different one-dimensional marginals of the posterior distribution of the hierarchical Bayesian mixture model described in Section 5.2. For one marginal, the natural stratification suffices. For the other, it does not, but a preliminary computation made with the natural stratification suggests a better choice of strata. We use this example to explain how to diagnose and correct problems related to poorly chosen strata: Our results will serve to guide the practice of stratified MCMC.

5.1 The natural stratification for tails and marginals

Here, we briefly explain how EMUS can be used to estimate tail probabilities and low-dimensional marginals of high-dimensional distributions. Let $\Omega\subset\mathbb{R}^{d}$; let $\pi$ be a probability distribution on $\Omega$; and let $\eta:\Omega\rightarrow\mathbb{R}$. Suppose that one wishes to estimate the very small tail probability $\mathrm{\mathbf{P}}[\eta(Z)\geq\varepsilon^{-1}]$. In this case, it is natural to stratify in $\eta$ only. That is, one may choose a partition of unity $\{\phi_{i}\}_{i=1}^{L}$ on $\mathbb{R}$ and define bias functions

[TABLE]

depending only on $\eta$. For a partition of unity, one might choose the regular grid of piecewise constant functions defined in Section 4.2. We refer to (31) as the natural stratification. To compute the tail probability, one uses EMUS to estimate $\pi(\mathrm{\mathbf{1}}_{[\varepsilon^{-1},\infty)}\circ\eta)$.

Computing marginal densities is similar; in fact, computing tails may be understood as a special case of computing a marginal density. Suppose now that $\eta:\Omega\rightarrow\mathbb{R}^{\ell}$. To estimate the marginal $\pi_{\eta}$ of $\pi$ in $\eta$, one chooses a partition of unity $\{\phi_{i}\}_{i=1}^{L}$ on $\mathbb{R}^{\ell}$, again defining bias functions by (31). One then uses EMUS to compute averages of histogram bins, which are functions of the form

[TABLE]

We have

[TABLE]

so for small $h$ the averages of the histogram bins approximate $\pi_{\eta}$.

By the argument in Section 4.2, EMUS with the natural stratification will be dramatically more efficient than direct sampling as long the biased distributions are no harder to sample than the target distribution $\pi$. Essentially, this is because, with the natural stratification, very small averages like $\mathrm{\mathbf{P}}[\eta(Z)\geq\varepsilon^{-1}]$ over the target distribution $\pi$ are expressed as functions of much larger averages over the biased distributions $\pi_{i}$. Unfortunately, however, for general functions $\eta$, the biased distributions of the natural stratification need not be easy to sample. In Section 5.3, we give one example where the natural stratification works and one where it does not. In the case where it does not, we explain how to make a better choice of strata.

5.2 A hierarchical Bayesian mixture model

Here, we review the hierarchical Bayesian mixture model proposed in [39], and we discuss the difficulties which complicate inference under this model. As a tutorial in the use of EMUS, we present a numerical investigation of these difficulties in Section 5.

In the hierarchical mixture model, the data vector $\mathbf{y}=(y_{1},\dots,y_{n})\in\mathbb{R}^{n}$ consists of independent identically distributed samples drawn from a mixture distribution of the form

[TABLE]

where $K$ is the number of mixture components, $q_{k}$ is the weight of the $k$’th mixture component, $\nu(\cdot;\mu_{k},\lambda_{k}^{-1})$ is the normal density with mean $\mu_{k}$ and variance $\lambda_{k}^{-1}$, and $\phi$ is the vector of parameters

[TABLE]

(Since $p(y_{i}|\phi)$ is a probability distribution, $q_{1}+\dots+q_{K}=1$, and ${q_{1},\dots,q_{K-1}}$ determine $q_{K}$.) The following prior distribution is imposed on $\phi$:

[TABLE]

As in [23, 10], we choose

[TABLE]

where $R$ and $M$ are the range and the mean of the observed data, respectively. The posterior density is

[TABLE]

where $\theta=(\phi,\beta)$ denotes the vector of all parameters to be inferred, including the hyperparameter $\beta$.

Several factors complicate inference based on this model. First, the mixture components are not identifiable; that is, the posterior distribution is invariant under permutation of the labels of the mixture components. Consequences of non-identifiability are discussed at length in [23, 10]. In our computations in Section 5.3, we impose the constraint

[TABLE]

to ensure that the components are identifiable. Second, in Lemma 5.1, we show that the posterior density may be unbounded, introducing spurious modes with infinite density. Finally, even with identifiability constraints, the posterior distribution may have multiple modes of finite posterior density. For example, see the modes reported in [10]. In Section 5.3, we use EMUS to efficiently visualize the posterior, assessing the effects of multimodality and unboundedness.

We suspect that the unboundedness of the posterior for this model is well known. However, we are unable to find a reference, so we now explain. It is certainly well known that the likelihood of a Gaussian mixture model is unbounded: Roughly speaking, the likelihood is infinite when any mixture component is collapsed on a single data point [1]. Nonetheless, one might expect the posterior density $p(\theta|\mathbf{y})$ to be bounded, since the prior penalizes large values of the precisions $\lambda_{i}$. This is not always the case when the data vector contains repeated entries:

Lemma 5.1**.**

If any datum $y_{i}$ has frequency $N_{i}$ greater than

[TABLE]

*then the posterior density $p(\theta|\mathbf{y})$ is unbounded. *

Proof 5.2**.**

*Take the limit of $p(\theta|\mathbf{y})$ as $\lambda_{1}\rightarrow\infty$ with $\mu_{1}=y_{i}$, $\beta=\lambda_{1}^{-1}$, and all other variables held fixed. *

The reader will observe that under the model, the set of data vectors with repeated entries has probability zero. However, in practice, the data consist of measurements with finite precision, and therefore repeated entries occur commonly, cf. the Hidalgo stamp data used in Section 5.3.

5.3 Numerical experiments: Choosing strata, computing tails, diagnosis of problems

In this section, we explain how to recognize and correct problems related to poor choices of strata, and we demonstrate the use of EMUS to investigate the multimodality and unboundedness of the posterior in the mixture model. We first compute two one-dimensional marginals of the high-dimensional posterior density $p(\theta|\mathbf{y})$ using the natural stratification (31). The natural stratification works in one case but not the other. In the case where the natural stratification does not work, preliminary calculations based on the natural stratification suggest a better choice of strata.

Here, we let $\mathbf{y}$ be the Hidalgo stamp data set first studied in [22], consisting of the thicknesses of $485$ stamps, ranging between $60\ \mu\rm{m}$ and $130\ \mu\rm{m}$. We let there be three mixture components ($K=3$), following previous computational studies [10, 23]. In our first calculation, we estimated the marginal in $\mu_{2}$ using the natural stratification with a grid of 201 bias functions covering the range $[7,11]$, with the support of the leftmost and rightmost bias functions reaching to $-\infty$ and $\infty$, respectively. For the middle strata, define $\phi_{1}:\mathbb{R}\rightarrow\mathbb{R}$ by

[TABLE]

We used the bias functions

[TABLE]

for $i=2,\dots,200$. Now, define $\phi_{2}:\mathbb{R}\rightarrow\mathbb{R}$ by

[TABLE]

The first and last bias functions were

[TABLE]

where $h=0.02$ as before.

We chose the total number of bias functions based on the sizes of the off-diagonal entries in the overlap matrix. For any bias functions of the form (34), the overlap matrix is tridiagonal. Thus, by Remark 3.12, if the superdiagonal and subdiagonal entries $F_{i,i+1}$ and $F_{i,i-1}$ are sufficiently large, then the EMUS estimator is not too sensitive to statistical errors in $\bar{F}$. For our choice of bias functions,

[TABLE]

We sampled the biased distributions using the affine invariant ensemble sampler with $100$ walkers, as implemented in the emcee package [15]. Due to computational restrictions on memory, only every tenth sample point was saved. As a check on the sampling, the average acceptance probability over all walkers in the ensemble sampler was calculated for each biased distribution. Averaging over biased distributions gave a total average acceptance probability of 0.31. The minimum acceptance probability over all distributions was 0.12.

To initialize sampling, we computed an unbiased test trajectory; that is, a trajectory having ergodic distribution $\pi$. We then started by sampling a single biased distribution $\pi_{k}$, initializing with points drawn randomly from the unbiased trajectory. We sampled the other biased distributions in sequence, initializing with points drawn randomly from samples of adjacent biased distributions. Thus, we sampled $\pi_{k}$ first, then $\pi_{k-1}$ and $\pi_{k+1}$, then $\pi_{k-2}$ and $\pi_{k+2}$, etc. We equilibrated the sampler in each $\pi_{i}$ for 3000 Monte Carlo steps, and collected data for an additional 100000 Monte Carlo steps. Each step of the ensemble sampler involves perturbing the positions of each of the 100 walkers.

We computed the marginal in $\mu_{2}$ using a grid of 200 histogram bins, covering the region $[7,11]$; this corresponds to taking ${h=0.01}$ in (32). The result is the curve labeled EMUS in Figure 3a. The marginal in $\mu_{2}$ has two modes, labeled $1$ and $2$ in Figure 3a. We plot the mixture distributions corresponding to these modes in Figure 4. (To be precise, the distributions in Figure 4 correspond to means over histogram bins centered at the labeled points.)

For comparison, we also estimated the marginal in $\mu_{2}$ from multiple long, unbiased trajectories. We computed $100$ unbiased trajectories of the affine invariant ensemble sampler in parallel. For each trajectory, the ensembles were first equilibrated for 10000 Monte Carlo steps, and then data were collected for 100000 steps. These trajectories were combined and binned to produce the density labeled Unbiased in Figure 4. We estimated the relative asymptotic variance of the marginal density for the unbiased calculation using ACOR [14], and we estimated the relative asymptotic variance for the EMUS calculation using the method outlined in Appendix G. We present the results in Figure 3a. Note that near the mode, unbiased MCMC performs slightly better than EMUS, but in the tails, EMUS performs dramatically better.

After computing the marginal in $\mu_{2}$, we tried computing the marginal in $\log_{10}\lambda_{1}$. We used the natural stratification with a grid of 50 bias functions with maxima equally spaced between –1 and 3.2 constructed as

[TABLE]

where

[TABLE]

We used the same initialization scheme as for the marginal in $\mu_{2}$, beginning with a single biased distribution initialized from an unbiased test trajectory. We call this the center sample. The result of this calculation was the density labeled “1D Center” in Figure 5a. When we tried to compute the asymptotic variance of this density estimate, we noticed very slow convergence of the sampler for some biased distributions. To investigate, we performed another EMUS calculation using a similar initialization procedure, but starting from $\pi_{1}$, the biased distribution at the extreme left, covering the lowest values of $\lambda_{1}$. We call this the left sample. The result of this second calculation was the density labeled “1D Left” in Figure 5a. For both the center and left samples, the strata were equilibrated for $3000$ steps and sampled for another $200000$. We observe that the two densities differ significantly in the region $-1\leq\log_{10}\lambda_{1}\leq 0.5$. They should be the same up to sampling errors; for example, we observe that different initializations have no effect on the calculation of the marginal in $\mu_{2}$, cf. Figure 3a.

Figure 6 explains the problem and suggests a solution: In the region $0.2\leq\log_{10}\lambda_{1}\leq 0.7$, the center and left samples cover entirely different ranges of $\log_{10}\lambda_{2}$. This suggests that the biased distributions corresponding to the range $0.2\leq\log_{10}\lambda_{1}\leq 0.7$ are multimodal, with barriers in $\lambda_{2}$ impeding sampling.

To confirm the hypothesis that barriers in $\lambda_{2}$ were responsible for the poor convergence observed in the center and left samples, we performed a third calculation, stratifying in both $\log_{10}\lambda_{1}$ and $\log_{10}\lambda_{2}$. We used a $50\times 50$ grid of bilinear bias functions, with maxima equally spaced between $-1$ and $3.2$. To be precise, for $i,j=1,\dots,50$, we defined the bias functions

[TABLE]

with $h$ as before. Let $\eta_{ij}$ denote the biased distribution corresponding to $\psi_{ij}$.

We performed the two-dimensional EMUS calculation twice, initializing from the center and left samples drawn from the natural stratification in $\log_{10}\lambda_{1}$. For each $i=1,\dots,L$, to sample the row $\{\eta_{ij}:j=1,\dots,50\}$ of biased distributions, we began by initializing sampling of a single biased distribution $\eta_{ik}$ with points from the either the center or left sample of $\pi_{i}$. We then sampled the other distributions $\eta_{ij}$ for $j\neq k$ in sequence, again initializing with points from samples of adjacent distributions, either $\eta_{i,j+1}$ or $\eta_{i,j-1}$ in this case. If no samples were found inside the support of a biased distribution, that distribution was ignored. For each biased distribution, sampling was burned in for 4500 steps, and samples were collected for an additional 2500 steps. Ultimately, 1397 of the 2500 biased distributions were sampled; the unsampled distributions correspond to the white space in Figure 7a.

We computed the marginal in $\log_{10}\lambda_{1}$ and $\log_{10}\lambda_{2}$ using a $200\times 200$ grid of histogram bins, covering the region ${-1\leq\log_{10}\lambda_{1}\leq 3.2}$ and ${-1\leq\log_{10}\lambda_{2}\leq 3.2}$; this corresponds to taking ${h=(3.2-(-1))/200}$ in (32); the result from the center calculation appears in Figure 7a. In Figure 8, we show the mixture distributions corresponding to the modes of the two-dimensinoal marginal in Figure 7a. The two-dimensional marginals were essentially the same for the center and left initializations; see Figure 9. We also estimated the one-dimensional marginal in $\log_{10}\lambda_{1}$ using the two-dimensional stratification; see the results labeled “2D Center” and “2D Left” in Figure 5a. Finally, we estimated the relative asymptotic variance of the marginal in $\log_{10}\lambda_{1}$ computed by two-dimensional stratification. Again, we observe that EMUS performs much better than unbiased sampling in the tails, cf. Figure 5b.

The marginal in $\log_{10}\lambda_{1}$ and $\log_{10}\lambda_{2}$ confirms that barriers in $\lambda_{2}$ caused the problems observed in calculating the marginal in $\log_{10}\lambda_{1}$ using the natural stratification. In fact, we see that computing the marginal in either $\lambda_{1}$ or $\lambda_{2}$ requires stratifying both variables, as stratifying only one leads to barriers that impede sampling in the other. In particular, there are barriers in $\lambda_{2}$ along the line $\log_{10}\lambda_{1}=0.45$ and a barrier in $\lambda_{1}$ along $\log_{10}\lambda_{2}=0.6$: In Figure 10, we plot an estimate of the conditional distribution of $\log_{10}\lambda_{2}$ with $\log_{10}\lambda_{1}=0.45$ fixed. This distribution is multimodal with a region of very low probability separating the modes, which explains the poor sampling depicted in Figure 6.

To conclude, we have confirmed that EMUS can be extremely efficient for computing tails. However, one must exercise care in the choice of strata. The natural stratification often suffices, but in some cases, like computing the marginal in $\log_{10}\lambda_{1}$, the biased distributions of the natural stratification may be very difficult to sample. We propose the use of different initializations, like the center and left samples, as a method of identifying problems related to poorly chosen strata. Careful inspection of simulations performed with these different initializations can identify problems and suggest better strata.

6 Conclusions

We have analyzed the Eigenvector Method for Umbrella Sampling (EMUS), an especially simple and effective stratified MCMC method sharing many features with the popular WHAM [27] and MBAR [44] methods of computational chemistry. We have demonstrated the advantages of EMUS for sampling from multimodal distributions and computing tail probabilities, and we have explained how to identify and resolve the problems which may occur if the method is implemented poorly. We have also given a tutorial intended to explain how to diagnose and correct problems related to poorly chosen strata.

Our purpose was to explain the benefits of stratified MCMC analytically, with the ultimate goal of introducing stratified MCMC to a diverse audience of statisticians, engineers, and scientists. Since stratified MCMC had previously been applied only to a particular class of statistical mechanics calculations without any general justification, we began by developing a general theory. We hope that our theory will serve as the basis for further developments. For example, it may now be possible to undertake a comparison of EMUS and other so-called reaction coordinate methods such as Wang–Landau sampling [52] or Metadynamics [28]. Despite some similarities with EMUS, these methods work by a substantially different mechanism and understanding the relative advantages of the two approaches is non-trivial. We also note that there is potential to apply stratification to problems that lie outside the scope of the present work. For example, we present a stratification method capable of computing dynamical quantities such as mean first passage times [11].

Appendix A Derivation of (1) and (3)

To see that (1) holds, observe that

[TABLE]

Therefore, we have

[TABLE]

To prove (3), observe that

[TABLE]

Appendix B Proof of Lemma 2.1

Proof B.1**.**

We prove only the second statement; proof of the first is similar. By definition, a non-negative matrix $M\in\mathbb{R}^{L\times L}$ is irreducible if and only if for every subset $A\subset\{1,2,\dots,L\}$ of the indices, there exist indices $i\in A$ and $j\notin A$ so that $M_{ji}>0$. Now assume that for every $A\subset\{1,2,\dots,L\}$, there exist $j\notin A$ and $t\geq 0$ so that

[TABLE]

*Then for some $i\in A$, $\psi_{i}(X^{j}_{t})>0$, so $\bar{F}_{ji}>0$, hence $\bar{F}$ is irreducible. *

Appendix C Sparse Grid of Strata

One may define a uniform grid of strata so that (15) increases only as $d^{2}$ with dimension, not exponentially: For any $\mathbf{i}\in\mathbb{Z}^{d}$, let $V_{\mathbf{i}}^{\prime}:=h\mathbf{i}+h\left[-\frac{1}{2},\frac{1}{2}\right]^{d}$. For $\mathbf{i}\neq\mathbf{j}$, define

[TABLE]

to be the $d$-dimensional pyramid consisting of all points in $V^{\prime}_{\mathbf{j}}$ closer to $V^{\prime}_{\mathbf{i}}$ than to any other cube $V^{\prime}_{\mathbf{k}}$. Now let $e_{n}$ denote the $n$’th standard basis vector in $\mathbb{R}^{d}$, and define

[TABLE]

to be the cube $V_{\mathbf{i}}^{\prime}$ enlarged by all the neighboring pyramids $W_{\mathbf{i}\mathbf{j}}$. The strata $V_{\mathbf{i}}$ are convex, and the corresponding bias functions $\psi_{\mathbf{i}}=\frac{1}{2}\mathrm{\mathbf{1}}_{V_{\mathbf{i}}}$ are a partition of unity. Each stratum $V_{\mathbf{i}}$ intersects only the $2d$ neighboring strata $V_{\mathbf{i}\pm e_{n}}$ for $n=1,\dots,d$. Moreover, each intersection between neighboring strata $V_{\mathbf{i}}$ and $V_{\mathbf{j}}$ consists of the pair of pyramids $W_{\mathbf{i}\mathbf{j}}$ and $W_{\mathbf{j}\mathbf{i}}$, and it has volume $1/d$. Therefore, by Lemma 3.10, for this choice of bias functions, the nonzero entries of $F$ decrease as $1/d$. It follows that (15) increases as $d^{2}$.

Appendix D Proofs of Theorem 3.1 and Theorem 3.4

Our proof of Theorem 3.1 (the CLT for EMUS) is based on the delta method. To apply the delta method, we require the following result ensuring the differentiability of $w(G)$:

Lemma D.1**.**

*The function $w(G)$ admits an extension $\tilde{w}:\mathbb{R}^{L\times L}\rightarrow\mathbb{R}^{L}$ which is differentiable on the set of irreducible stochastic matrices. *

Proof D.2**.**

*By [47, Lemma 3.1], $w(G)$ admits a continuously differentiable extension to an open set $U\subset\mathbb{R}^{L\times L}$. We further extend the domain of $w(G)$ to $\mathbb{R}^{L\times L}$ by arbitrarily defining $w(G)=0$ whenever $G\in\mathbb{R}^{L\times L}\setminus U$. *

The extension in Lemma D.1 resolves two technicalities: First, the set of stochastic matrices is not a vector space but a compact, convex subset of $\mathbb{R}^{L\times L}$ with empty interior. Therefore, the derivative of $w$ is undefined. Second, $\bar{F}$ may be reducible for some values of $N$ and some realizations of the processes sampling the biased distributions. In that case, the invariant distribution of $\bar{F}$ is not unique, so $w(\bar{F})$ is undefined. Throughout the remainder of this work, $w(G)$ will denote the extension guaranteed by the lemma.

We now prove the CLT for EMUS.

Proof D.3** (Proof of Theorem 3.1).**

The proof is based on the delta method [6, Proposition 6.2] and a formula for $w^{\prime}(\bar{F})$ given in [19].

By Lemma D.1, $w(\bar{F})$ is differentiable at $F$, so the function

[TABLE]

is differentiable at $\left(F,\{\pi_{i}[g^{\ast}]\}_{i=1}^{L},\{\pi_{i}[\mathrm{\mathbf{1}}^{\ast}]\}_{i=1}^{L}\right)$. Let $\partial_{i}B\in\mathbb{R}^{L+2}$ be the derivative of $B$ with respect to those quantities computed from $X^{i}_{t}$: That is,

[TABLE]

where $\frac{\partial B}{\partial\bar{F}_{i:}}\in\mathbb{R}^{L}$ denotes the partial derivative of $B$ with respect to the $i$’th row of $\bar{F}$ and

[TABLE]

To simplify notation, we will assume throughout the remainder of this argument that all derivatives are evaluated at $\left(F,\{\pi_{i}[g^{\ast}]\}_{i=1}^{L},\{\pi_{i}[\mathrm{\mathbf{1}}^{\ast}]\}_{i=1}^{L}\right)$. In formulas involving matrix multiplication, we will treat $\partial_{i}B$, $\frac{\partial B}{\partial\bar{F}_{i:}}$, and $\frac{\partial B}{\partial\bar{G}_{i:}}$ as row vectors.

Since we assume that the processes $X^{i}_{t}$ sampling the different measures $\pi_{i}$ are independent, [5, Chapter 1, Theorem 2.8] implies that

[TABLE]

where $\Sigma$ is the covariance matrix of the product of the distributions $N\left(0,\kappa_{i}^{-1}\Sigma_{i}\right)$. (That is, ${\Sigma\in\mathbb{R}^{L(L+2)\times L(L+2)}}$ is the block diagonal matrix with the matrices $\kappa_{i}^{-1}\Sigma_{i}$ along the diagonal.) Therefore, by the delta method,

[TABLE]

where

[TABLE]

Now we observe that for any column vector $v\in\mathbb{R}^{L}$ having mean zero,

[TABLE]

by [19, Theorem 3.1]. (In the formula above, $e_{i}\in\mathbb{R}^{L}$ denotes the $i$’th standard basis vector.) Therefore, we have

[TABLE]

where

[TABLE]

(Equality (42) above follows from (1) and the definition (LABEL:eq:_definition_of_Psi) of $\Psi$.) Also,

[TABLE]

Thus,

[TABLE]

*and the result follows by (41). *

We now prove of Theorem 3.4. To begin, we give some upper bounds on the partial derivatives of the weight vector $w(F)$ with respect to the entries of the overlap matrix $F$.

Definition D.4**.**

Let $e_{i}\in\mathbb{R}^{L}$ denote the $i$’th standard basis vector. For $i,j\in\{1,2,\dots,L\}$ with $i\neq j$, define the logarithmic partial derivatives

[TABLE]

*(These partial derivatives must be understood as derivatives of the extension guaranteed by Lemma D.1; otherwise, they are defined only when $F_{ij}>0$ and $F_{ii}>0$.) *

Our definition of logarithmic partial derivatives in (D.4) is not standard. However, we observe that a version of the standard formula relating the total and partial derivatives of $\log w$ holds: For all matrices $H$ whose rows sum to zero,

[TABLE]

We need only consider matrices whose rows sum to zero, since these are the only perturbations for which $F+\varepsilon H$ can be stochastic.

The following result appears in [47, Theorem 3.6]. It is crucial in our proof of Theorem 3.4.

Lemma D.5**.**

Recall $\mathrm{\mathbf{P}}_{i}[t_{j}<t_{i}]$ and $\frac{\partial\log w_{k}}{\partial{F}_{ij}}$ from Definitions 3.3 and D.4. For all stochastic and irreducible matrices $F$, we have

[TABLE]

We also require the following lemma in the proof of Theorem 3.4.

Lemma D.6**.**

The asymptotic covariance matrix $\sigma^{i}$ has the properties:

The rows and colums of $\sigma^{i}$ sum to zero. That is, for $e\in\mathbb{R}^{L}$ the vector of all ones,

[TABLE] 2. 2.

For all $j=1,\dots,L$,

[TABLE]

Proof D.7**.**

Since the rows of $\bar{F}$ sum to one with probability one, we have

[TABLE]

for any fixed number of samples $N_{i}$. Therefore, the asymptotic variance $\sigma^{i}$ has $e^{\mathrm{t}}\sigma^{i}e=0$, and it follows that $e^{\mathrm{t}}\sigma^{i}=\sigma^{i}e=0$ since $\sigma^{i}$ is symmetric and positive semidefinite.

Let $k$ be such that $F_{ik}=0$. Since $\bar{F}_{ik}=0$ with probability one, we have

[TABLE]

*for any $j=1,\dots,L$, and therefore $\sigma^{i}_{jk}=0$. *

We now prove Theorem 3.4.

Proof D.8** (Proof of Theorem 3.4).**

We begin with formula (41):

[TABLE]

Since the asymptotic covariance matrix $\Sigma_{i}$ is symmetric and positive semidefinite, the Cauchy inequality holds:

[TABLE]

for all $a,b\in\mathbb{R}^{L+1}$. Therefore,

[TABLE]

(Here, $\mathbf{0}$ denotes the zero vector in $\mathbb{R}^{L}$, interpreted as a column vector.)

We now estimate the term $A_{0}$ defined above. By (42), we have

[TABLE]

where

[TABLE]

(The third equality above follows from formula (45) relating the total and partial derivatives of $\log w$, since the rows and colums of $\sigma^{i}$ sum to zero by Lemma D.6.)

We claim that

[TABLE]

To prove this, we observe that $R^{i}$ is symmetric and positive semidefinite since $\sigma^{i}$ is symmetric and positive semidefinite. Therefore, $R^{i}$ has the spectral decomposition

[TABLE]

with eigenvalues $\lambda_{i,j}>0$ and corresponding eigenvectors $v^{i,j}$ such that $\lVert v^{i,j}\rVert=1$. Thus, for any $a\in\mathbb{R}^{L}$,

[TABLE]

Inequality (49) follows from (50) by setting

[TABLE]

Finally, combining (48), (49), and Lemma D.5 yields

[TABLE]

Moreover, we have

[TABLE]

by (1), and therefore

[TABLE]

We now observe that by (43)

[TABLE]

$\operatorname{C}(\bar{h}_{i})$* denotes the asymptotic covariance of the trajectory average $\bar{h}_{i}$ of $h$ over the biased process $X^{i}_{t}$. Therefore, combining (46) and (51), we find*

[TABLE]

*as desired. *

When the bias functions are a partition of unity, both the EMUS method and the statements of Theorems 3.1 and 3.4 simplify considerably. (The bias functions are a partition of unity if and only if $\sum_{i=1}^{L}\psi_{i}(x)=1$ for all $x$.) In this case, $f^{\ast}=f$ for all functions $f$, and the EMUS method reduces to

[TABLE]

where

[TABLE]

Corollary D.9**.**

*Suppose that the bias functions are a partition of unity. In that case, Theorem 3.4 holds with either $\operatorname{var}_{\pi}(g)$ or $\pi[\lvert g\rvert]^{2}$ in place of $\pi[\lvert h\rvert]^{2}$. In addition, $\Psi=1$, and one can replace $\operatorname{C}(\bar{h}_{i})$ with the asymptotic variance $\operatorname{C}(\bar{g}_{i})$ of $\bar{g}_{i}$. *

Proof D.10**.**

When the bias functions are a partition of unity,

[TABLE]

*and so we may replace $\pi[\lvert h\rvert]^{2}$ with $\operatorname{var}_{\pi}(g)$. In addition, equation (42) holds with $\mathfrak{g}_{k}=\pi_{k}[g]$. Thus, following the argument above, one may verify that the result also holds with $\pi[\lvert g\rvert]^{2}$ in place of $\pi[\lvert h\rvert]^{2}$. *

Appendix E Proof of Theorem 3.6

In the arguments below, for any probability measure $\nu$ on a set $\Omega$, we let

[TABLE]

and we define the $L^{2}(\nu)$ inner product

[TABLE]

with the corresponding norm

[TABLE]

Given a set $U\subset\mathbb{R}^{d}$, we define $L^{2}(U)$, $\lVert\cdot\rVert_{L^{2}(U)}$, $\langle,\rangle_{U}$ to be the analogous function space, norm, and inner product for Lebesgue measure on $U$.

Our proof of Theorem 3.6 requires a Poincaré inequality, Lemma E.1. We refer to [31, Section 3] for an introduction to Poincaré inequalities and their role in the theory of diffusion processes.

Lemma E.1**.**

Assume that the Poincaré inequality holds for $U$ with constant $\Lambda$; that is, assume that for all weakly differentiable $f:U\rightarrow\mathbb{R}$ so that $\nabla f\in L^{2}(U)$,

[TABLE]

We have a similar Poincaré inequality for $\pi_{h}$:

[TABLE]

Proof E.2**.**

By a standard scaling argument, the Poincaré inequality holds for $U_{h}$ with constant $h\Lambda$. To see this, let $A_{h}:U\rightarrow U_{h}$ be the affine transformation

[TABLE]

For any $f:U_{h}\rightarrow\mathbb{R}$ with $\nabla f\in L^{2}(U_{h})$, using the change of variable formula and the chain rule, we have

[TABLE]

Now observe that for any $f\in L^{2}(\pi_{h})$,

[TABLE]

since $\pi_{h}[f]$ is the $L^{2}(\pi_{h})$ orthogonal projection of $f$ onto the space of constant functions. Therefore, we have

[TABLE]

*and the result follows. *

Remark E.3**.**

*The Poincaré inequality for the Lebesgue measure on a set $U$ holds under very weak conditions on $U$. For example, when $U$ is convex, the Poincaré inequality holds with constant $\Lambda(U)=D/\pi$, where $D$ is the diameter of the domain [38]. *

We now prove Theorem 3.6:

Proof E.4** (Proof of Theorem 3.6).**

We begin by stating a simple consequence of the functional central limit theorem for reversible, continuous time Markov processses: Let $Y_{t}$ be a reversible, stationary Markov process with ergodic distribution $\pi$ and generator $L$. Let $g\in L^{2}(\pi)$, and define

[TABLE]

By [25, Corollary 1.9],

[TABLE]

where

[TABLE]

Here, $L^{-1}(g-\pi[g])$ denotes any function in the domain of $L$ with

[TABLE]

and $\pi[L^{-1}(g-\pi[g])]=0$. Such a function must exist when $g\in L^{2}(\pi)$ and $X_{t}$ is reversible [25].

We now show that the process $X^{h}_{t}$ meets the conditions above for the central limit theorem. First, we recall that the generator of $X^{h}_{t}$ is the operator

[TABLE]

with domain

[TABLE]

see [2, Proposition 3.2] for the case of a convex polyhedron or [13, Chapter 8] for a domain with $C^{3}$ boundary.

By [24, Theorem 4.3.3], a process $Y_{t}$ with invariant distribution $\pi$ is reversible if its generator is symmetric and it has the strong continuity property

[TABLE]

where $T_{t}f(x):=\mathrm{\mathbf{E}}_{x}[f(Y_{t})]$ denotes the backwards semigroup associated with $Y_{t}$. The generator $L_{h}$ of $X^{h}_{t}$ is symmetric, since for all $f,g\in D(L_{h})$, using integration by parts, we have

[TABLE]

(Here, $z_{h}^{-1}:=\int_{U_{h}}\exp(-\beta V)\,dx$ is the normalizing constant for $\pi_{h}$.) Since $\langle\nabla f,\nabla g\rangle_{\pi}$ is invariant under exchanging $f$ and $g$, $\langle f,L_{h}g\rangle_{\pi}=\langle L_{h}f,g\rangle_{\pi}$ and $L_{h}$ is symmetric. We postpone discussion of the strong continuity of $X^{h}_{t}$ to the end of the proof.

We now use the Poincaré inequality (Lemma E.1) and (54) to prove that $X^{h}_{t}$ is ergodic and to estimate the term $L_{h}^{-1}(g-\pi_{h}[g])$ appearing in the formula for $\sigma^{2}_{h}(g)$; in essence, we adapt the approach outlined in [31, Section 3] to the family of reflected processes $X^{h}_{t}$. We prove ergodicity first. By [4, Proposition 2.2], a process is ergodic if and only if [math] is a simple eigenvalue of its generator. By the Poincaré inequality (Lemma E.1) and (54), for all $u\in D(L_{h})$,

[TABLE]

where

[TABLE]

Now if $u$ is not constant, $\lVert u-\pi_{h}[u]\rVert_{L^{2}(\pi_{h})}^{2}>0$, so $\lVert L_{h}u\rVert_{L^{2}(\pi_{h})}>0$ and $u$ is not an eigenvector with eigenvalue [math]. Hence, [math] is a simple eigenvalue of $L_{h}$, and $X^{h}_{t}$ is ergodic.

Finally, we estimate $\sigma^{2}_{h}(g)$. We have

[TABLE]

Taking $u=L_{h}^{-1}(g-\pi_{h}[g])$ in the above yields

[TABLE]

which implies

[TABLE]

using the Cauchy–Schwarz inequality.

It remains to show that the process $X^{h}_{t}$ has the strong continuity property (53). We only sketch an argument, since the basic ideas are standard. First, one can use the Lipschitz continuity of strong solutions of the reflected process [2, Lemma 4.1] to show that $X^{h}_{t}$ has the Feller property. (That is, one can show that $T_{t}u$ is continuous whenever $u$ is continuous.) In addition, since the process $X^{h}_{t}$ has an infinitesimal generator, we have the pointwise continuity property

[TABLE]

for all $x\in U_{h}$ and all $u\in D(L_{h})$. Now we have $\lVert T_{t}\rVert_{\infty}\leq 1$ for all $t\geq 0$, where $\lVert T_{t}\rVert_{\infty}$ is the operator norm of $T_{t}$ on the space of continuous functions with the sup-norm, and therefore by a density argument the limit (56) holds for all continuous $u$. Hence, by [8, Lemma 1.4], we have

[TABLE]

*for all continuous $u$. The strong continuity property (53) then follows by another density argument, using that $\lVert T_{t}\rVert_{L^{2}(\pi_{h})}\leq 1$ for all $t\geq 0$. *

Appendix F Proof of Theorems 4.1 and 4.2

Proof F.1** (Proof of Theorem 4.1).**

By Corollary D.9, since the bias functions are a partition of unity, we have

[TABLE]

To prove the desired upper bound, we substitute estimates of $C(\bar{g}_{\mathbf{i}})$, $R^{\mathbf{i}}$, and $F_{\mathbf{i}\mathbf{j}}$ into the inequality above.

First, we consider the asymptotic covariances $R^{\mathbf{i}}$ and $C(\bar{g}_{\mathbf{i}})$. Let

[TABLE]

The diameter of $U_{\mathbf{i}}$ is $2\sqrt{d}h$, so by Assumption 3.2.1

[TABLE]

(The second inequality follows since $h\beta\leq 1$ by definition.) Similarly,

[TABLE]

Second, by Lemma 3.10, the nonzero entries of the overlap matrix $F$ are bounded below as $\beta$ tends to infinity:

[TABLE]

for all $\mathbf{i},\mathbf{j}$ so that $F_{\mathbf{i},\mathbf{j}}>0$. We also observe that each row of $F$ has $3^{d}$ nonzero entries, since $F_{\mathbf{i},\mathbf{i}+\mathbf{k}}>0$ only when all entries of $\mathbf{k}$ belong to $\{-1,0,1\}$.

We now estimate the term involving $\operatorname{C}(\bar{g}_{\mathbf{i}})$ in (57). By (59), we have

[TABLE]

Now we have

[TABLE]

Therefore,

[TABLE]

(The inequality follows since $0\leq z_{\mathbf{i}}\leq 1$ for all $\mathbf{i}$; the second to last equality follows using (4) and that $\{\psi_{\mathbf{i}}\}_{\mathbf{i}\in\mathbb{Z}^{d}/K\mathbb{Z}^{d}}$ is a partition of unity.) Thus,

[TABLE]

It remains to address the term involving $R^{\mathbf{i}}$ in (57): Using (58), (60), and that each row of $F$ has $3^{d}$ nonzero entries, we have

[TABLE]

for every $\mathbf{i}\in\mathbb{Z}^{d}/K\mathbb{Z}^{d}$. Finally, using (63), (64), and $\kappa_{\mathbf{i}}^{-1}=K^{d}=\lceil\beta\rceil^{d}$, we conclude

[TABLE]

*where the constants $D$ and $E$ depend on $d$ and $V$, but not on $g$ or $\beta$. *

We note that if one uses the bias functions proposed in Appendix C, then the constants $D$ and $E$ in the proof of Theorem 4.1 grow only polynomially with the dimension $d$, not exponentially. However, we do not claim that those bias functions perform better than the uniform grid (17) or the bias functions of Section 5.3 in practice.

We now prove Theorem 4.3:

Proof F.2** (Proof of Theorem 4.3).**

Take $g:=\mathrm{\mathbf{1}}_{x\geq M}$. Since the bias functions are a partition of unity, by Corollary D.9, we have

[TABLE]

First, we estimate

[TABLE]

for all $i=1,\dots,K-1$ by Assumption 3.2.1. By Assumption 4.2,

[TABLE]

and $R^{K}_{jj}=0$ for $j\neq K-1,K$, since $\psi_{K+1}$ is constant over the support of $\pi_{K}$. In addition,

[TABLE]

since all bias functions $\psi_{i}$ take a constant value over the support of $\pi_{K+1}$. Likewise,

[TABLE]

and $\operatorname{C}(\bar{g}_{i})=0$ for all $i\neq K$.

We now show that the nonzero entries of the overlap matrix are bounded below independent of $M$. First, we estimate the entries which are averages over the biased distributions with bounded support. By Lemma 3.10, we have

[TABLE]

for all ${i=0,\dots,K-1}$ and $j$ so that $F_{ij}>0$. whenever $F_{ij}>0$. It remains to address those entries related to biased distributions with unbounded support, so with $i=K,K+1$. By Lemma F.3, $F_{K,K+1}$ and $F_{K,K-1}$ are bounded below by some $\theta>0$ independent of $M$, for this choice of bias functions. (Lemma F.3 and its proof appear in Appendix F. Lemma F.3 is the only part of the proof which relies on Assumption 4.2.) In addition, for any $i=0,\dots,K+1$, we have $F_{ii}=\frac{1}{2}$, which implies $F_{K+1,K}=1-F_{K+1,K+1}=\frac{1}{2}$ since $F$ is stochastic when the bias functions are a partition of unity.

Finally, we substitute the above estimates of the overlap matrix and the variances into (65). Let $c=\min\{\theta,1/2\exp(2)\}$. Observe that $h$ decreases with $M$, so $Ceh^{a}\leq E$ for some constant $E$, uniformly in $M$. Let $F=\max\{D,E\}$. We have

[TABLE]

We now observe that

[TABLE]

Therefore,

[TABLE]

*which proves the result. *

We now prove Lemma F.3, which is used in the proof of Theorem 4.3.

Lemma F.3**.**

Under the hypotheses of Theorem 4.3, there exist constants ${M_{1},\theta_{+},\theta_{-}>0}$ depending on $V$ but not on $M$ so that

[TABLE]

*whenever $M\geq M_{1}$. *

Proof F.4**.**

We consider $F_{K,K-1}$ first. We have

[TABLE]

By the integral mean value theorem,

[TABLE]

for some $\xi_{M-h,M}\in[M-h,M]$. Moreover, by (28), we have

[TABLE]

Therefore, when $M-h\geq M_{0}$,

[TABLE]

It follows that

[TABLE]

using the definition $h=M/K$.

To estimate the quotient in expression (66), we distinguish two cases: By (28), $V^{\prime}$ is nondecreasing on $[M_{0},\infty)$, so either ${\lim_{x\rightarrow\infty}V^{\prime}(x)=C_{2}<\infty}$ or ${\lim_{x\rightarrow\infty}V^{\prime}(x)=\infty}$. In the first case, $V^{\prime}$ is bounded, and we have

[TABLE]

whenever $M-h\geq M_{0}$. In the second case, for $M$ sufficiently large,

[TABLE]

Therefore, applying in succession the mean value theorem, the monotonicity of $V^{\prime}$, assumption (29), and the hypothesis $\lim_{x\rightarrow\infty}V^{\prime}(x)=\infty$, we have that for all $M$ sufficiently large,

[TABLE]

(In the second and third lines above, $\eta_{M-h,M}\in[M-h,M]$ denotes the point guaranteed by the mean value theorem so that ${V^{\prime}(M)-V^{\prime}(M-h)=hV^{\prime\prime}(\eta_{M-h,M})}$.) It follows from (66), (67), and (68) that there exist $M_{-},\theta_{-}>0$ so that

[TABLE]

whenever $M\geq M_{-}$.

Now we prove that $F_{K,K+1}$ is bounded below. We have

[TABLE]

As above, to bound the quantity appearing in the exponent in (70), we distinguish the two cases $\lim_{x\rightarrow\infty}V^{\prime}(x)=C_{1}<\infty$ and ${\lim_{x\rightarrow\infty}V^{\prime}(x)=\infty}$. In the first case, for $M$ sufficiently large that ${2C_{1}\geq\lvert V^{\prime}(x)\rvert\geq C_{1}/2}$ whenever ${x\geq M-h}$, we have

[TABLE]

In the second case, for $M$ sufficiently large,

[TABLE]

By (29), we have the differential inequality

[TABLE]

This implies

[TABLE]

for

[TABLE]

the solution of the initial value problem

[TABLE]

Therefore,

[TABLE]

so by (72),

[TABLE]

It follows from (70), (71), and (73) that there exist $M_{+},\theta_{+}>0$ so that

[TABLE]

*whenever $M\geq M_{+}$. *

Appendix G Improved Method of Computing Error Bars

In [48, Section VII.B.1], we proposed a practical method of estimating the asymptotic standard deviations (error bars) of averages computed by EMUS. Using the notation established in Appendix D, our method proceeds as follows:

Compute $\bar{F}$, $\{\bar{g}^{\ast}_{i}\}_{i=1}^{L}$, and $\{\bar{1}^{\ast}_{i}\}_{i=1}^{L}$. 2. 2.

Compute $w(\bar{F})$ and the group inverse $(I-\bar{F})^{\#}$. 3. 3.

Evaluate $\partial_{i}B$ at $\bar{F}$, $\{\bar{g}^{\ast}_{i}\}_{i=1}^{L}$, and $\{\bar{\mathrm{\mathbf{1}}}^{\ast}_{i}\}_{i=1}^{L}$. 4. 4.

Compute the time series

[TABLE] 5. 5.

Compute an estimate $\bar{\chi}_{i}^{2}$ of the integrated autocovariance of $\bar{\zeta}^{i}_{t}$ using an algorithm such as ACOR [14]. 6. 6.

Compute as an estimate $\sigma^{2}$ the quantity

[TABLE]

We originally proposed computing the group inverse $(I-\bar{F})^{\#}$ using the method of [19] based on the QR factorization. We have since discovered that this method does not always yield sufficiently accurate results. For example, when computing error bars for the marginal in $\mu_{2}$ in Section 5.3, we observed a highly oscillatory numerical error affecting some entries of $(I-\bar{F})^{\#}$. That the sign pattern in Figure 11a fails to be symmetric is evidence of this numerical error. We note that since the exact overlap matrix $F$ is in detailed balance with $w(F)$, we have $\operatorname*{diag}(w(F))F\operatorname*{diag}(w(F))^{-1}=F^{\mathrm{t}}$. (Here, $\operatorname*{diag}(w(F))$ denotes the diagonal matrix with $w(F)$ along the diagonal.) Therefore,

[TABLE]

which implies that the sign pattern of $(I-F)^{\#}$ is symmetric since $w(F)$ is positive. As a result of these numerical errors, we were unable to accurately compute error bars for the EMUS estimate of the marginal density.

We therefore propose computing the group inverse by a new method combining QR factorization with power iteration. We first compute an estimate $G^{0}$ of $(I-\bar{F})^{\#}$ by the method of [19]. We then iterate

[TABLE]

where $e\in\mathbb{R}^{L}$ denotes the column vector of all ones and $\tilde{F}:=(I-ew(F)^{t})\bar{F}$. We observe that $(I-\bar{F})^{\#}$ is a fixed point of this iteration, since

[TABLE]

Above, we use well known properties of the group inverse, including that the spectral projector $I-ew(\bar{F})^{t}$ commutes with $\bar{F}$, that $(I-ew(\bar{F})^{t})(I-\bar{F})^{\#}=(I-F)^{\#}$, and that $(I-\bar{F})(I-\bar{F})^{\#}=I-ew(\bar{F})^{t}$.

Moreover, when $\bar{F}$ is irreducible, $\mathscr{I}^{K}$ is a contraction for $K$ sufficiently large. By the Perron-Frobenius theorem, the spectral radius of $\tilde{F}$ is smaller than $1-\varepsilon$, for some $\varepsilon>0$. Therefore, by Gelfand’s formula, for any matrix norm $\left\lVert\cdot\right\rVert$, we have ${\lim_{k\rightarrow\infty}\left\lVert\tilde{F}^{k}\right\rVert^{1/k}<1-\varepsilon/2}$, and so for some $K$,

[TABLE]

Now

[TABLE]

Thus, assuming that the norm $\lVert\cdot\rVert$ is submultiplicative,

[TABLE]

Therefore, the power iteration converges and its limit is the group inverse $(I-\bar{F})^{\#}$.

Using this new method, we computed $(I-\bar{F})^{\#}$ for $\bar{F}$ the overlap matrix involved in estimating the marginal in $\mu_{2}$ in Section 5.3. We performed $10^{6}$ power method iterates. Observe that the sign pattern of the group inverse computed with power iteration is symmetric; see Figure 11b.

The power iteration (76) converges slowly when the spectral gap of $\bar{F}$ is small. We have shown in [47] that the spectral gap may be very small: It decreases exponentially with a temperature parameter in a limit similar to the one analyzed in Section 4.1 above. However, even when the spectral gap is small, we conjecture that a modest number of power iterations will significantly reduce the numerical error in the group inverse, since the error in the initial calculation seems to be highly oscillatory and the power iteration has a smoothing effect.