An ergodic theorem for the weighted ensemble method

TL;DR

This paper proves an ergodic theorem for the weighted ensemble method, ensuring its reliability for long-term average computations in Markov chain sampling, and provides variance formulas for optimization.

Contribution

It establishes the first ergodic theorem for weighted ensemble, confirming its consistency for long-time averages and deriving variance formulas for potential improvements.

Findings

Proves ergodic theorem for weighted ensemble

Derives explicit variance formulas

Confirms method's consistency for long-term sampling

Abstract

We study weighted ensemble, an interacting particle method for sampling distributions of Markov chains that has been used in computational chemistry since the 1990s. Many important applications of weighted ensemble require the computation of long time averages. We establish the consistency of weighted ensemble in this setting by proving an ergodic theorem for time averages. As part of the proof, we derive explicit variance formulas that could be useful for optimizing the method.

Equations

$$\displaystyle\begin{split}\begin{array}[]{ccc}\text{parents }\{\xi_{t}^{i}\}^{i=1,\ldots,N}&\text{parents' weights }\{\omega_{t}^{i}\}^{i=1,\ldots,N}&\textup{user-chosen}\\ \downarrow\text{selection}&\downarrow\text{selection}&\textup{parameters}\\ \text{children }\{\hat{\xi}_{t}^{i}\}^{i=1,\ldots,N}&\text{children's weights }\{\hat{\omega}_{t}^{i}\}^{i=1,\ldots,N}&\leftarrow\\ \downarrow\text{mutation}&\downarrow\text{mutation}\\ \text{new parents }\{\xi_{t+1}^{i}\}^{i=1,\ldots,N}&\text{new parents' weights }\{\omega_{t+1}^{i}\}^{i=1,\ldots,N}\end{array}\end{split}$$

$${\omega }_t(u)=\sum _{i:{\xi }_t^i\in u}{\omega }_t^i.$$

$$\{N_t^i:{\xi }_t^i\in u\}\sim \text{Multinomial}\left(N_t(u),\left\{\frac{{\omega }_t^i}{{\omega }_t(u)}:{\xi }_t^i\in u\right\}\right).$$

$${\hat{\omega }}_t^j=\frac{{\omega }_t(u)}{N_t(u)},\text{if }{\hat{\xi }}_t^j\in u.$$

$$({\xi }_{t+1}^1,\dots ,{\xi }_{t+1}^N)\sim K({\hat{\xi }}_t^1,\cdot )\times \dots \times K({\hat{\xi }}_t^N,\cdot ).$$

$${\omega }_{t+1}^j={\hat{\omega }}_t^j,j=1,\dots ,N.$$

$$Pr(\text{sample }{\xi }_t^j\text{ in bin }u)=\frac{{\omega }_t^j}{{\omega }_t(u)},$$

$${\hat{\omega }}_t^j=\frac{{\omega }_t(u)}{N_t(u)},\text{if }{\hat{\xi }}_t^j\in u.$$

$$\mathbb{E}\left[\sum _{i=1}^N{\omega }_t^{i}g({\xi }_t^i)\right]=\int K^{t}gd\nu ,$$

$$\Vert K^t(\xi ,\cdot )-\mu (\cdot ){\Vert }_{TV}\le c{\lambda }^t,\text{for all }\xi \text{ and all }t\ge 0.$$

$$\underset{T\to \infty }{\lim }\frac{1}{T}\sum _{t=0}^{T-1}\sum _{i=1}^N{\omega }_t^{i}f({\xi }_t^i)=\int fd\mu .$$

$${\eta }_t^u=\sum _{i:{\xi }_t^i\in u}\frac{{\omega }_t^i}{{\omega }_t(u)}{\delta }_{{\xi }_t^i},$$

$$h_{t,T}=\sum _{s=0}^{T-t-1}{K}^{s}f.$$

$$\eta (g)=\int gd\eta ,{\text{Var}}_{\eta }(g)=\eta (g^2)-\eta (g{)}^2,$$

$$\text{Var}\left(\frac{1}{T}\sum _{t=0}^{T-1}\sum _{i=1}^N{\omega }_t^{i}f({\xi }_t^i)\right)$$

$$=\frac{1}{T^2}\text{Var}\left(\sum _{i=1}^N{\omega }_0^i{h}_{0,T}({\xi }_0^i)\right)\text{(initial condition variance)}$$

$$+\frac{1}{T^2}\sum _{t=0}^{T-2}\mathbb{E}\left[\sum _u\frac{{\omega }_t(u{)}^2}{N_t(u)}{\text{Var}}_{{\eta }_t^u}(K{h}_{t+1,T})\right]\text{(selection variance)}$$

$$+\frac{1}{T^2}\sum _{t=0}^{T-2}\mathbb{E}\left[\sum _u\frac{{\omega }_t(u{)}^2}{N_t(u)}{\eta }_t^u({\text{Var}}_K{h}_{t+1,T})\right]\text{(mutation variance)}.$$

$$\frac{1}{T^2}\sum _{t=0}^{T-2}\mathbb{E}\left[\sum _u\frac{{\omega }_t(u{)}^2}{N_t(u)}{\eta }_t^u({\text{Var}}_K{h}_{t+1,T})\right].$$

$$N_t(u)\approx \frac{N{\omega }_t(u)\sqrt{{\eta }_t^u({\text{Var}}_K{h}_{t+1,T})}}{{\sum }_u{\omega }_t(u)\sqrt{{\eta }_t^u({\text{Var}}_K{h}_{t+1,T})}}.$$

$$\approx \frac{1}{N{T}^2}\sum _{t=0}^{T-2}\mathbb{E}\left[{\left(\sum _u{\omega }_t(u)\sqrt{{\eta }_t^u({\text{Var}}_K{h}_{t+1,T})}\right)}^2\right].$$

$$\frac{1}{N{T}^2}\sum _{t=0}^{T-2}\mathbb{E}\left[\sum _u{\omega }_t(u){\eta }_t^u({\text{Var}}_K{h}_{t+1,T})\right],$$

$$K(1,2)=K(2,3)=\delta ,K(1,1)=K(2,1)=1-\delta ,K(3,1)=1,$$

$$\frac{1}{T^2}\sum _{t=0}^{T-2}\mathbb{E}\left[\sum _{i=1}^3\frac{{\omega }_t(i{)}^2}{N_t(i)}{\text{Var}}_K{h}_{t+1,T}(i)\right]\approx \frac{1}{T}\sum _{i=1}^3\frac{\mu (i{)}^2}{N/3}{\text{Var}}_K{h}_{t+1,T}(i)\approx \frac{6{\delta }^3}{NT},$$

$$\frac{1}{N{T}^2}\sum _{t=0}^{T-2}\nu (K^t({\text{Var}}_K{h}_{t+1,T}))\approx \frac{1}{NT}\mu ({\text{Var}}_K{h}_{t+1,T})\approx \frac{{\delta }^2-{\delta }^3}{NT}$$

$$\text{par}({\hat{\xi }}_t^j)={\xi }_t^i⟺{\xi }_t^i\text{ is the parent of }{\hat{\xi }}_t^j.$$

$$\mathbb{E}\left[{\omega }_{t+1}^{i}g({\xi }_{t+1}^i)|{\hat{\mathcal{F}}}_t\right]={\hat{\omega }}_t^{i}Kg({\hat{\xi }}_t^i).$$

$$\mathbb{E}\left[\sum _{i:{\hat{\xi }}_t^i\in u}{\hat{\omega }}_t^{i}g({\hat{\xi }}_t^i)|{\mathcal{F}}_t\right]=\sum _{i:{\xi }_t^i\in u}{\omega }_t^{i}g({\xi }_t^i).$$

$$\mathbb{E}\left[\sum _{i:{\hat{\xi }}_t^i\in u}{\hat{\omega }}_t^{i}g({\hat{\xi }}_t^i)|{\mathcal{F}}_t\right]$$

$$=\sum _{i:{\xi }_t^i\in u}\frac{{\omega }_t(u)}{N_t(u)}g({\xi }_t^i)\mathbb{E}[N_t^i|{\mathcal{F}}_t]=\sum _{i:{\xi }_t^i\in u}{\omega }_t^{i}g({\xi }_t^i).$$

$$\mathbb{E}\left[\sum _{i=1}^N{\omega }_{t+1}^{i}g({\xi }_{t+1}^i)|{\hat{\mathcal{F}}}_t\right]=\sum _{i=1}^N{\hat{\omega }}_t^{i}Kg({\hat{\xi }}_t^i),$$

$$\mathbb{E}\left[\sum _{i=1}^N{\hat{\omega }}_t^{i}g({\hat{\xi }}_t^i)|{\mathcal{F}}_t\right]=\sum _{i=1}^N{\omega }_t^{i}g({\xi }_t^i).$$

$$\mathbb{E}\left[\sum _{i=1}^N{\omega }_t^{i}g({\xi }_t^i)|{\mathcal{F}}_0\right]=\sum _{i=1}^N{\omega }_0^i{K}^{t}g({\xi }_0^i).$$

Full text

\authornames

David Aristoff

\authorone

[Colorado State University]David Aristoff

\addressone

841 Oval Drive, Fort Collins, CO 80523 \emailonearistoff@math.colostate.edu

\ams

65C0565C35;65C40

1 Introduction

Weighted ensemble [6, 11, 16, 17, 23, 25, 35, 39, 40, 52, 53, 56] is an importance sampling method, based on interacting particles, for distributions associated with a Markov chain. In this article, we focus on sampling the average of a function with respect to the steady state distribution of a generic Markov chain. By generic, we mean that the only thing we might know about the Markov chain is how to sample it; in particular, we may not know its stationary distribution up to a normalization factor.

Weighted ensemble consists of a collection of evolving particles with associated weights. In this sense, weighted ensemble can be understood as a kind of sequential Monte Carlo method [29, 18, 19, 21, 33, 4, 28, 49, 50]. In weighted ensemble, the particles evolve between selection steps according to the law of the underlying Markov chain. In each selection step, some of the particles are copied while others are killed; the resulting particles are given new weights so that weighted ensemble is statistically unbiased [53].

The selection step is based on dividing the particles into bins, where the particles in each bin are resampled according to their relative weights. In practice, the binning, and the number of copies in each bin, should be chosen so that important particles survive and irrelevant particles are killed. The definition of the bins, and how many copies to maintain in each, requires some care. With appropriate choices, weighted ensemble can have drastically smaller variance than direct Monte Carlo, or independent particles; see the references above, or [55] for a more complete list.

Weighted ensemble was developed for applications in computational chemistry [35] ranging from state space exploration [23] to protein association [35] and protein folding [57]. One important application we have in mind is the computation of the mean time for a protein to unfold [6]. This time can be reformulated as the inverse of the steady state flux into the unfolded state of the underlying Markovian molecular dyamics, with an added sink in the unfolded state and source in the folded state [34]. This dynamics can approach its steady state on time scales significantly smaller than the mean unfolding time [54]. As the flux into the unfolded state is usually very small, importance sampling is needed to estimate it with substantial precision [6, 44].

Other unbiased methods, differing from weighted ensemble in that they usually sample finite-time quantities rather than ergodic averages, include Adaptive Multilevel Splitting [7, 8, 10], Forward Flux Sampling [1], and some sequential Monte Carlo methods [21, 13, 50]. This unbiased property allows for a relatively straightforward study of variance using martingale techniques [2, 8, 18, 21]. In this article, we extend these techniques to study the long-time stability of weighted ensemble.

Our main contribution here is a proof of the consistency of weighted ensemble via an ergodic theorem. We believe that this is the first ergodic theorem for an interacting particle system in which the interactions come from resampling.

A secondary contribution comes from explicit formulas for the variance of weighted ensemble at finite particle number. The proof of the ergodic theorem is a straighforward consequence of these formulas. On the theoretical side, our variance formulas are handy for understanding the rate of weighted ensemble convergence, and on the practical side, they could be used for optimizing the method. We mostly leave this discussion to other works, including our companion paper [3]; see also [2] and the references above.

This article is organized as follows. In Section 2, we describe weighted ensemble in detail. In Section 3, we state our main results, including the unbiased property (Theorem 3.1), the ergodic theorem (Theorem 3.2), and the variance formulas (Theorem 3.3). In Section 4, we compare weighted ensemble to direct Monte Carlo, and give a simple example illustrating the potential gain. All of our proofs are in Section 5.

2 Description of the method

Weighted ensemble consists of a fixed number, $N$, of particles belonging to a common state space, each carrying a positive scalar weight, and undergoing repeated selection and mutation steps. In the selection step, some of the particles are copied, and others are killed, according to a stratification or binning scheme. In the mutation step, the particles evolve according to an underlying Markov kernel $K$.

At time $t$ before selection, the particles, called parents, are $\xi_{t}^{1},\ldots,\xi_{t}^{N}$. At time $t$ after selection, the particles, called children, are $\hat{\xi}_{t}^{1},\ldots,\hat{\xi}_{t}^{N}$. The weights of the parents and children are $\omega_{t}^{1},\ldots,\omega_{t}^{N}$ and $\hat{\omega}_{t}^{1},\ldots,\hat{\omega}_{t}^{N}$, respectively. The following diagram illustrates weighted ensemble evolution:

[TABLE]

The initial particles $\xi_{0}^{1},\ldots,\xi_{0}^{N}$ can be arbitrary. The initial weights must be strictly positive and sum to one: $\omega_{0}^{i}>0$ for all $i$, and $\omega_{0}^{1}+\ldots+\omega_{0}^{N}=1$. The children are initially just copies of their parents, but they evolve forward in time conditionally independently. When we say conditionally independent, we mean conditional on the $\sigma$-algebra representing the information from (1) up to the current time.

Weighted ensemble requires the user to choose, before the selection step at time $t$, a collection of nonempty bins that partition the set of parents, as well as a particle allocation that defines the number of children in each bin. We write $u$ for the bins, and $N_{t}(u)\geq 1$ for the number of children in bin $u$ at time $t$. We require that $\sum_{u}N_{t}(u)=N$. The bins can change in time, but for simpler notation we leave this implicit.

In the selection step, the children in each bin are obtained by sampling with replacement from the parents in the bin, according to their weight distribution, as many times as the particle allocation specifies. The children’s weights in each bin are all the same after selection, and the total weight in the bin is preserved [16].

In more detail, define the total weight in bin $u$ at time $t$ as

[TABLE]

The numbers, $N_{t}^{i}$, of children of the parents $\xi_{t}^{i}\in u$ are conditionally multinomial:

[TABLE]

Children are assigned to the same bins as their parents, with weights

[TABLE]

Selections in distinct bins are conditionally independent.

In the mutation step, the children evolve conditionally independently via $K$:

[TABLE]

The weights do not change during the mutation step. Thus

[TABLE]

We summarize weighted ensemble in the following algorithm.

**Algorithm **\thetheorem

Choose initial weights $\omega_{0}^{1},\ldots,\omega_{0}^{N}>0$ summing to $1$ and initial particles $\xi_{0}^{1},\ldots,\xi_{0}^{N}$. Then iterate over $t\geq 0$:

• Selection step**

Partition the parents $\xi_{t}^{1},\ldots,\xi_{t}^{N}$ into a collection of bins.

• 2.

Assign a number $N_{t}(u)\geq 1$ of children to the parents in each bin $u$.

• 3.

Sample $N_{t}(u)$ children from the parents in bin $u$, with replacement, using

[TABLE]

where $\omega_{t}(u):=\sum_{i\,:\,\xi_{t}^{i}\in u}\omega_{t}^{i}$ is the total weight in bin $u$.

• 4.

Give all the children in bin $u$ the same weight

[TABLE]

• Mutation step**

Evolve the children $\hat{\xi}_{t}^{1},\ldots,\hat{\xi}_{t}^{N}$ conditionally independently using $K$ to get the next parents $\xi_{t+1}^{1},\ldots,\xi_{t+1}^{N}$. Keep the same weights $\omega_{t+1}^{j}=\hat{\omega}_{t}^{j}$, $j=1,\ldots,N$.

2.1 Algorithm details and remarks

A few remarks are in order to clarify Algorithm 2.

• •

We abuse notation by writing $\xi_{t}^{i}\in u$ or $\hat{\xi}_{t}^{j}\in u$ to indicate that $\xi_{t}^{i}$ or $\hat{\xi}_{t}^{j}$ is in bin $u$, even though the bins form a partition of the particles, and not (necessarily) a partition of the state space.

• •

A child is simply a copy of its parent: if $\hat{\xi}_{t}^{j}$ is a child of $\xi_{t}^{i}$, then $\hat{\xi}_{t}^{j}=\xi_{t}^{i}$. The indices of the children are not important, so specifying the number of children of each parent is enough to define them.

• •

Since the weights don’t change in the mutation step and the the selection step preserves the total weight, the total weight is constant in time: $\omega_{t}^{1}+\ldots+\omega_{t}^{N}=1$ for all $t\geq 0$. We discuss the importance of this in Remark 5.11.

• •

We assume that the bins and particle allocation at time $t$ are included in ${\mathcal{F}}_{t}$, the $\sigma$-algebra generated by the information from Algorithm 2 just before the $t$-th selection step. We also write $\hat{\mathcal{F}}_{t}$ for the $\sigma$-algebra generated by the information from Algorithm 2 just after the $t$-th selection step. See Section 5.1 for details.

• •

We assume multinomial sampling in the bins because it leads to simple explicit variance expressions in terms of intrabin variances. In Remark 5.13, we comment on residual sampling, which performs much better than multinomial resampling and still admits nice variance formulas.

For our ergodic theorem, the bins and particle allocation can be arbitrary. To actually do better than direct Monte Carlo, they must be judiciously chosen. The most common strategy is to define bins based on a carefully constructed partition of state space – particles occupy the same bin when they belong to the same element of the partition – and then allocate children approximately uniformly among these bins. Some knowledge about the underlying problem is needed to choose the bins, but this strategy has had considerable success, as the references in the introduction attest (see [55] for a mostly current list of application papers). We propose a different strategy in our companion paper [3] that uses our variance analysis below. We summarize that strategy in Section 4 below.

3 Main results

3.1 Unbiased property

We begin with the unbiased property of weighted ensemble. This property was previously noted in [53], and proved in a slightly different setting in [2].

Theorem 3.1** **(Unbiased property)

For each $t\geq 0$ and all bounded measurable $g$,

[TABLE]

where $\nu$ is the weighted ensemble initial distribution, $\int g\,d\nu:={\mathbb{E}}[\sum_{i=1}^{N}\omega_{0}^{i}g(\xi_{0}^{i})]$.

We could interpret Theorem 3.1 as follows. If $(X_{t})$ is a Markov chain with kernel $K$ and initial distribution $\nu$, then ${\mathbb{E}}\left[\sum_{i=1}^{N}\omega_{t}^{i}g(\xi_{t}^{i})\right]={\mathbb{E}}[g(X_{t})]$. In this sense, weighted ensemble gives unbiased estimates of the law of the underlying Markov chain.

3.2 Ergodic theorem

To ensure that weighted ensemble is ergodic, the underlying Markov kernel $K$ must be ergodic in some sense. We assume that $K$ is uniformly ergodic [27]:

{assumption}

There is $c>0$, $\lambda\in[0,1)$ and a probability measure $\mu$ such that

[TABLE]

Here and below, $f$ is a fixed bounded measurable function.

Theorem 3.2** **(Ergodic theorem)

If Assumption 3.2 holds, then with probability $1$,

[TABLE]

Convergence of the mean of the time average in (7), at the same rate as direct Monte Carlo, follows from Assumption 3.2 and the unbiased property (Theorem 3.1). For the ergodic theorem to hold, and for weighted ensemble to beat direct Monte Carlo, the variance of the time average should be sufficiently small. Well-behaved variance is not automatic for unbiased methods; see Remark 5.11 below.

3.3 Variance formulas

Here, we give exact, finite $N$ formulas for the variance of weighted ensemble, based on a martingale decomposition. To get nice concise formulas, we need some notation. Define the intrabin distributions

[TABLE]

where $\delta_{\xi}$ is the Dirac delta distribution centered at $\xi$. Define also

[TABLE]

For a probability measure $\eta$ and bounded measurable function $g$, define

[TABLE]

and in particular, let $\textup{Var}_{K}g(\xi):=\textup{Var}_{K(\xi,\cdot)}(g)=Kg^{2}(\xi)-(Kg)^{2}(\xi)$.

Theorem 3.3** **(Variance formulas)

For each time $T>0$,

[TABLE]

The expression in (9) can be interpreted as the variance coming from the initial condition, while the expressions in (10) and (11) can be understood as the variances arising from each selection and mutation step, respectively.

Using Theorem 3.3, the proof of the ergodic theorem is straightforward. Under Assumption 3.2, we can show that variances of $h_{t,T}$ and $Kh_{t,T}$ are uniformly bounded in $t$ and $T$. This makes the weighted ensemble variance $\textup{Var}\left(\frac{1}{T}\sum_{t=0}^{T-1}\sum_{i=1}^{N}\omega_{t}^{i}f(\xi_{t}^{i})\right)=O(1/T)$. The ergodic theorem then follows from standard arguments. Note that the variance expressions (8)-(11) by themselves do not require Assumption 3.2.

Beyond the ergodic theorem, these variance formulas are interesting in their own right, since they could be used to design binning and particle allocation schemes that minimize the weighted ensemble variance. Indeed, that is what we have done in our companion paper [3] (see also [2]). We discuss this more in the next section.

4 Comparison to direct Monte Carlo

There are many existing works showing that weighted ensemble can provide significant gains over direct Monte Carlo: see for instance the references list in [55], and our companion paper [3]. The main goal of this article is to prove the consistency of weighted ensemble from explicit variance formulas, and not to reaffirm this point. We include, however, a brief discussion here.

Weighted ensemble works by reducing the mutation variance (11), compared to that of direct Monte Carlo (see (29) below), via the selection step. This comes however at the cost of a positive selection variance (10), compared to direct Monte Carlo which has selection variance equal to zero. Thus, weighted ensemble beats direct Monte Carlo if the reduction in mutation variance is greater than the selection variance cost. (The initial condition variances may be ignored, since they are $O(1/T^{2})$ for weighted ensemble and direct Monte Carlo, while the overall variances are $O(1/T)$.)

In more detail, the weighted ensemble mutation variance is

[TABLE]

A Lagrange multiplier calculation shows that the expression inside the expectation in (12) is minimized over $N_{t}(u)$, subject to the constraint $\sum_{u}N_{t}(u)=N$, when

[TABLE]

Plugging this into (12) makes the weighted ensemble mutation variance

[TABLE]

By Jensen’s inequality, (14) is less than or equal to

[TABLE]

which, by the unbiased property of weighted ensemble (Theorem 3.1), is exactly the direct Monte Carlo mutation variance (see Remark 5.12 below). This form of optimal mutation variance gain was originally observed in [2, Remark 4.1].

The selection variance (10) is small whenever, at each time $t$, bins $u$ are chosen inside which $Kh_{t+1,T}$ does not vary too much. With enough particles and bins, it is possible to keep the selection variance arbitrarily small, while also controlling the mutation variance by keeping the particle allocation close to the optimal (13). Even with modest numbers of particles and bins, weighted ensemble has proven useful; see the list of applications in [55], most of which use relatively small $N$.

In our companion paper [3], we propose an optimization strategy based on choosing the particle allocation to minimize mutation variance, and the bins to minimize selection variance. There, the particle allocation is a simplified version of (13), in which we use the limit $h:=\lim_{T\to\infty}(h_{t,T}-(T-t)\int f\,d\mu)$ in place of $h_{t+1,T}$, and the bins are chosen to make the intrabin variances of $Kh$ small. Of course, estimating $Kh$ and $\textup{Var}_{K}h$ is a difficult problem. In [3] we propose using a Markov state model to get “cheap” approximations of $Kh$ and $\textup{Var}_{K}h$; such models are already commonly used for preconditioning weighted ensemble simulations [6, 14, 15].

4.1 Example

Below is a simple example illustrating the variance reduction in weighted ensemble, compared to direct Monte Carlo, in the context of our variance analysis above. Consider a Markov chain on $3$ states, with transition matrix

[TABLE]

where $\delta>0$ is small, and we take $f(1)=f(2)=0$ and $f(3)=1$.

For weighted ensemble, we will assign particles to the same bin if and only if they occupy the same point in space. Following the usual method in the literature, we allocate an approximately equal number of children to each bin.

In this case, the variance of weighted ensemble can be estimated as

[TABLE]

where the first approximation in (16) holds for large enough $N$ and $T$, and the second approximation uses direct calculations, dropping terms of higher order than $\delta^{3}$. The variance of direct Monte Carlo (see Remark 5.12) can be estimated by

[TABLE]

for large $T$. Figure 1 shows numerical confirmation of these estimates. Note the smaller variance (higher order in $\delta$) for weighted ensemble, compared to direct Monte Carlo.

5 Derivations

5.1 Notation

Below, ${\mathcal{F}}_{t}$ is the $\sigma$-algebra generated by the parents and their weights from times $0\leq s\leq t$, the children and their weights from times $0\leq s\leq t-1$, and the bins and particle allocations from times $0\leq s\leq t$. Meanwhile, $\hat{\mathcal{F}}_{t}$ is the $\sigma$-algebra generated by ${\mathcal{F}}_{t}$ together with the children and their weights at time $t$. Throughout, $g$ denotes a bounded measurable function, and $c$ a positive constant whose value can change between different equations. We will use the notation

[TABLE]

5.2 One-step means

Lemma 5.1

For each $i=1,\ldots,N$ and $t\geq 0$,

[TABLE]

At each time $t\geq 0$, for each bin $u$,

[TABLE]

Proof 5.2

By (5), ${\mathbb{E}}[g(\xi_{t+1}^{i})|\hat{\mathcal{F}}_{t}]=Kg(\hat{\xi}_{t}^{i})$. Now (19) follows from this and (6). Meanwhile, by (3), ${\mathbb{E}}[N_{t}^{i}|{\mathcal{F}}_{t}]=N_{t}(u)\omega_{t}^{i}/\omega_{t}(u)$ if $\xi_{t}^{i}\in u$. Thus by (4) and (18),

[TABLE]

Lemma 5.3** **(One-step means)

For each time $t\geq 0$,

[TABLE]

Proof 5.4

This follows immediately from Lemma 5.1, by summing over the particles in (19) to get (21), and summing over the bins in (20) to get (22).

5.3 Proof of the unbiased property

Proof 5.5** **(Proof of Theorem 3.1)

By repeated application of Lemma 5.3 with the tower property,

[TABLE]

Taking expectations in (23) gives the result.

5.4 Doob martingale and variance decomposition

Below, define the Doob martingale

[TABLE]

Proposition 5.6

For $0\leq t\leq T-1$, we have

[TABLE]

Proof 5.7

This comes from Lemma 5.3 by repeated application of the tower property.

Proposition 5.8** **(Martingale variance decomposition)

For each $T>0$,

[TABLE]

Proof 5.9

It is straightforward to check that all the martingale differences $D_{t+1}-\hat{D}_{t}$ and $\hat{D}_{t}-D_{t}$ are uncorrelated with each other and with $D_{0}$. The proof is finished by writing $D_{T-1}=\sum_{t=0}^{T-1}\sum_{i=1}^{N}\omega_{t}^{i}f(\xi_{t}^{i})$ as a telescoping sum of the martingale differences, computing ${\mathbb{E}}[D_{T-1}^{2}]$ in terms of the martingale differences, and subtracting ${\mathbb{E}}[D_{T-1}]^{2}={\mathbb{E}}[D_{0}]^{2}$ from the resulting expression.

In the proof of Theorem 3.3 below, we will see that the initial condition variance, selection variance, and mutation variance in (25) are the same as those in (9)-(11).

5.5 Proof of the variance formulas

Proof 5.10** **(Proof of Theorem 3.3)

Using the formula (24) for the Doob martingale, together with the one-step mean formula (21), the weight update formula (6), and the conditional independence of particle evolution in (5),

[TABLE]

By (26), the weight update formula (4), and the bin mean formula (20),

[TABLE]

In light of Proposition 5.8, this gives (11). Now we turn to (10). Using the formula (24) for the Doob martingale, the one-step mean formula (22), and the fact that selections in distinct bins are conditionally independent,

[TABLE]

Using (27), the weight update formula (4), and the conditional independence property of multinomial resampling (see e.g. [26], equation (6)),

[TABLE]

Taking expectations in this expression and appealing to Proposition 5.8 gives (10).

5.6 Remarks on the variance formulas

Remark 5.11

We briefly comment on the variance for other methods of the selection and mutation type (1). Consider a method with the same mutation step, but a different selection step that is still unbiased in the sense of (22).

In this case, the same variance decomposition (25) applies, with selection variance

[TABLE]

For a method like weighted ensemble in which the total weight is always $1$, $Kh_{t+1,T}$ can be replaced with $Kh_{t+1,T}-(T-t-1)\int f\,d\mu$ in (28) without otherwise changing the equation. Assumption 3.2 shows that $Kh_{t+1,T}-(T-t-1)\int f\,d\mu$ is uniformly bounded in $t$ and $T$. As a result, the selection variance (28) is $O(1/T)$, and the ergodic theorem remains valid.

We observed numerically that if the total weight varies at each time, then the weights tend to approach zero, and the variance (8) is of order $T$, as $T\to\infty$. Indeed, $Kh_{t+1,T}$ is typically of order $T$ as $T\to\infty$, which suggests that in this case, the selection variance (28) is on the order of $\sum_{t=0}^{T-2}{\mathbb{E}}[\textup{Var}(\sum_{i=1}^{N}\hat{w}_{t}^{i}|{\mathcal{F}}_{t})]$ as $T\to\infty$. Of course, the ergodic theorem fails if the variance (8) goes to infinity as $T\to\infty$.

Remark 5.12

Note that direct Monte Carlo – which we define as independent, equally-weighted particles – is a special case of weighted ensemble in which each particle always has weight $1/N$, every parent always gets its own bin $u$, and $N_{t}(u)$ is always $1$. In this case the selection variance is zero, while the mutation variance is

[TABLE]

Remark 5.13

If the selections in the bins use residual multinomial resampling [26] instead of multinomial resampling, then the selection variance (10) becomes

[TABLE]

where $r_{t}(u)$ and $\gamma_{t}^{u}$ are defined from the residuals r_{t}^{i}=\frac{N_{t}(u)\omega_{t}^{i}}{\omega_{t}(u)}-\big{\lfloor}\frac{N_{t}(u)\omega_{t}^{i}}{\omega_{t}(u)}\big{\rfloor} by

[TABLE]

We omit proof, but include this formula in case it is useful for optimizations, since residual multinomial resampling performs much better than multinomial resampling, yet still admits simple explicit variance expressions.

5.7 Proof of the ergodic theorem

Lemma 5.14

If Assumption 3.2 holds, then as $T\to\infty$,

[TABLE]

Proof 5.15

By Assumption 3.2, we have

[TABLE]

where now $c>0$ is a different constant. Thus

[TABLE]

This shows that $\textup{Var}_{\eta}h_{t+1,T}\leq C^{2}$, and similarly $\textup{Var}_{\eta}Kh_{t+1,T}\leq C^{2}$, for any probability distribution $\eta$. As a result, the selection and mutation variances (10)-(11) are both $O(1/T)$ as $T\to\infty$. By similar arguments, $\textup{Var}(\sum_{i=1}^{N}\omega_{0}^{i}h_{0,T}(\xi_{0}^{i}))\leq C^{2}$, which makes the initialization variance (9) $O(1/T^{2})$. Thus, the variance in (8) is $O(1/T)$.

Proof 5.16** **(Proof of Theorem 3.2)

Define $\theta_{T}=\frac{1}{T}\sum_{t=0}^{T-1}\sum_{i=1}^{N}\omega_{t}^{i}f(\xi_{t}^{i})$. Then for $0<S\leq T$,

[TABLE]

By Lemma 5.14, $\textup{Var}(\theta_{T})=O(1/T)$. So by Chebyshev’s inequality, there is $c>0$ so that

[TABLE]

for large enough $T$. With $T_{n}=n^{2}$, by the Borel-Cantelli lemma, there is $n_{0}$ such that

[TABLE]

By the unbiased property and Assumption 3.2,

[TABLE]

Now given $S>0$, we can choose $T_{n}$ so that $T_{n}\leq S\leq T_{n+1}$ and write

[TABLE]

By (31)-(33), with probability $1$, the right hand side of (34) vanishes as $S\to\infty$.

\acks

The author would like to acknowledge Frédéric Cérou, Peter Christman, Josselin Garnier, Gideon Simpson, Gabriel Stoltz and Brian Van Koten for helpful comments, and especially Jeremy Copperman, Matthias Rousset, Robert J. Webber, and Dan Zuckerman for interesting discussions and insights. The author thanks Robert J. Webber for pointing out errors and making many helpful suggestions concerning a previous version of the manuscript.

\fund

The author also gratefully acknowledges support from the National Science Foundation via the awards NSF-DMS-1818726 and NSF-DMS-1522398.