Advanced Multilevel Monte Carlo Methods

TL;DR

This paper reviews advanced Multilevel Monte Carlo techniques, focusing on their application with Markov chain and sequential Monte Carlo methods when exact sampling is not feasible, highlighting strategies to reduce computational costs.

Contribution

It introduces strategies for applying MLMC to MCMC and SMC methods without requiring exact sampling, expanding MLMC's applicability.

Findings

Strategies for MLMC in MCMC and SMC without exact sampling

Reduction in computational cost for biased expectation estimation

Enhanced applicability of MLMC in complex probabilistic models

Abstract

This article reviews the application of advanced Monte Carlo techniques in the context of Multilevel Monte Carlo (MLMC). MLMC is a strategy employed to compute expectations which can be biased in some sense, for instance, by using the discretization of a associated probability law. The MLMC approach works with a hierarchy of biased approximations which become progressively more accurate and more expensive. Using a telescoping representation of the most accurate approximation, the method is able to reduce the computational cost for a given level of error versus i.i.d. sampling from this latter approximation. All of these ideas originated for cases where exact sampling from couples in the hierarchy is possible. This article considers the case where such exact sampling is not currently possible. We consider Markov chain Monte Carlo and sequential Monte Carlo methods which have been introduced in the literature and we describe different strategies which facilitate the application of MLMC within these methods.

Figures

Tables

Table 1. Table 1: The key rates of convergence required for MLSMC samplers.

Rate parameter	Relevant quantity
$\alpha$	$\left({𝔼}_{{\pi }_L}-{𝔼}_{\pi }\right)\left(\phi \right)$
$\beta$	${sup}_{u\in 𝖤}\left|\frac{{\kappa }_l\left(u\right)Z_{l-1}}{{\kappa }_{l-1}\left(u\right)Z_l}-1\right|$

Table 2. Table 2: The key rates of convergence required for PMMH using MLMC.

Rate parameter	Relevant quantity
$\alpha$	$\left({𝔼}_{{\pi }_L^k}-{𝔼}_{{\pi }^k}\right)\left(\phi \right)$
$\beta$	${\int }_{\Theta \times {ℝ}^{2kd}}{\left|\phi (\theta ,{\overline{u}}_{1:k})-\phi (\theta ,{\underset{¯}{u}}_{1:k})\right|}^2{\check{\pi }}_{l-1:l}^k(\theta ,{\check{u}}_{1:k})𝑑\theta 𝑑{\check{u}}_{1:k}$

Table 3. Table 3: The key rates of convergence required for MLPF.

Rate parameter	Relevant quantity
$\alpha$	$\left({𝔼}_{{\pi }_L^k}-{𝔼}_{{\pi }^k}\right)\left(\phi \right)$
$\beta$	${\left({\int }_{{ℝ}^{2d}}{\left|\phi \left({\overline{u}}_k\right)-\phi \left({\underset{¯}{u}}_k\right)\right|}^2{\check{\pi }}_{l-1:l}^k\left({\check{u}}_k\right)𝑑{\check{u}}_k\right)}^2$

Table 4. Table 4: The key rates of convergence required for MLEnKF, for all p ≥ 1 𝑝 1 p\geq 1 , where π ^ ˇ l − 1 : l k superscript subscript ˇ ^ 𝜋 : 𝑙 1 𝑙 𝑘 \check{\widehat{\pi}}_{l-1:l}^{k} denotes the coupled measure resulting from the algorithm above.

Rate parameter	Relevant quantity
$\alpha$	$\left({𝔼}_{{\hat{\pi }}_L^k}-{𝔼}_{{\pi }^k}\right)\left(\phi \right)$
$\beta$	${\left({\int }_{{ℝ}^{2d}}{\left|\phi \left({\overline{u}}_k\right)-\phi \left({\underset{¯}{u}}_k\right)\right|}^p{\check{\hat{\pi }}}_{l-1:l}^k({\overline{u}}_k,{\underset{¯}{u}}_k)𝑑{\overline{u}}_k𝑑{\underset{¯}{u}}_k\right)}^{2/p}$

Equations

$${\mathbb{E}}_{\pi }[\phi (U)]={\int }_{\mathsf{E}}\phi (u)\pi (du)$$

$$\pi (du)=\frac{\kappa (u)}{Z}du,$$

$$\mathbb{E}\Big{[}\Big{(}\frac{1}{N}\sum_{i=1}^{N}\varphi(U^{i})-\mathbb{E}_{\pi}[\varphi(U)]\Big{)}^{2}\Big{]},$$

$$\frac{1}{N}\mathbb{V}{\text{ar}}_{{\pi }_L}[\phi (U)]+({\mathbb{E}}_{{\pi }_L}[\phi (U)]-{\mathbb{E}}_{\pi }[\phi (U)]{)}^2,$$

$${\mathbb{E}}_{{\pi }_L}[\phi (U)]=\sum _{l=1}^L\{{\mathbb{E}}_{{\pi }_l}-{\mathbb{E}}_{{\pi }_{l-1}}\}[\phi (U)],$$

$$\mathbb{V}{\text{ar}}_{({\pi }_l,{\pi }_{l-1})}[\phi (U_l(l))-\phi (U_{l-1}(l))]=\mathcal{O}(h_l^{\beta }).$$

$$\frac{1}{N_l}\sum _{i=1}^{N_l}\{\phi (u_l^i(l))-\phi (u_{l-1}^i(l))\}\approx \{{\mathbb{E}}_{{\pi }_l}-{\mathbb{E}}_{{\pi }_{l-1}}\}[\phi (U)].$$

$$\frac{1}{N_{1}}\sum_{i=1}^{N_{1}}\varphi(u_{1}^{i}(1))+\sum_{l=2}^{L}\Big{(}\frac{1}{N_{l}}\sum_{i=1}^{N_{l}}\{\varphi(u_{l}^{i}(l))-\varphi(u_{l-1}^{i}(l))\}\Big{)}.$$

$$\frac{1}{N_1}\mathbb{V}{\text{ar}}_{{\pi }_1}[\phi (U)]+\sum _{l=2}^L\frac{1}{N_l}\mathbb{V}{\text{ar}}_{({\pi }_l,{\pi }_{l-1})}[\phi (U_l(l))-\phi (U_{l-1}(l))]+({\mathbb{E}}_{{\pi }_L}[\phi (U)]-{\mathbb{E}}_{\pi }[\phi (U)]{)}^2,$$

$$L\propto -\frac{\log (ϵ)}{\alpha \log (2)}.$$

$$\sum _{l=1}^L\frac{h_l^{\beta }}{N_l}\propto {ϵ}^2.$$

$$N_l={ϵ}^{-2}{h}_l^{(\beta +\zeta )/2}{K}_L,$$

$$-\nabla \cdot (\hat{u}\nabla p)$$

$$p$$

$$\hat{u}(x)=\bar{u}(x)+\sum _{k=1}^K{u}_k{\sigma }_k{\varphi }_k(x).$$

$$\mathcal{G}(p)=[g_1(p),\cdots ,g_M(p){]}^{\top },$$

$$Y=\mathcal{G}(p)+\xi ,\xi \sim \mathcal{N}(0,\Gamma ),\xi \perp u,$$

$$\kappa (u)=e^{-\Phi [\mathcal{G}(p(\cdot ;u))]},\Phi (\mathcal{G})=\frac{1}{2}|\mathcal{G}-y{|}_{\Gamma }^2.$$

$$\pi_{l}(u)=\frac{\kappa_{l}(u)}{Z_{l}}=\frac{e^{-\Phi[\mathcal{G}(p^{l}(\cdot;u))]}}{\int_{E}e^{-\Phi[\mathcal{G}(p^{l}(\cdot;u))]}du}\ ,\quad l=1,\ldots,L.\$$

$$d{U}_t$$

$$\prod _{i=1}^{n}G(u_i,y_i)Q^{\infty }(u_{(i-1)},u_i),$$

$${\pi }_{\infty }(u_k|y_{1:k})={\pi }_{\infty }^k(u_k)=\frac{{\int }_{{\mathbb{R}}^{(k-1)d}}{\prod }_{i=1}^{k}G(u_i,y_i)Q^{\infty }(u_{(i-1)},u_i)d{u}_{1:k-1}}{{\int }_{{\mathbb{R}}^{kd}}{\prod }_{i=1}^{k}G(u_i,y_i)Q^{\infty }(u_{(i-1)},u_i)d{u}_{1:k}},$$

$$Z_{\infty }(y_{1:k})=Z_{\infty }^k{\int }_{{\mathbb{R}}^{kd}}\prod _{i=1}^{k}G(u_i,y_i)Q^{\infty }(u_{(i-1)},u_i)d{u}_{1:k}.$$

$$U_k^l(m+1)$$

$${\xi }_k(m)$$

$${\pi }_L(u_k|y_{1:k})=\frac{{\int }_{{\mathbb{R}}^{(k-1)d}}{\prod }_{i=1}^{k}G(u_i,y_i)Q^L(u_{(i-1)},u_i)d{u}_{1:k-1}}{{\int }_{{\mathbb{R}}^{kd}}{\prod }_{i=1}^{k}G(u_i,y_i)Q^L(u_{(i-1)},u_i)d{u}_{1:k}},$$

$$Z_L(y_{1:k})={\int }_{{\mathbb{R}}^{kd}}\prod _{i=1}^{k}G(u_i,y_i)Q^L(u_{(i-1)},u_i)d{u}_{1:k}.$$

$$d{U}_t=a_{\theta }(U_t)dt+b_{\theta }(U_t)d{W}_t,$$

$${\pi }_{\infty }(d\theta |y_{1:k})={\pi }_{\theta }(d\theta )\frac{{\int }_{{\mathbb{R}}^{kd}}{\prod }_{i=1}^k{G}_{\theta }(u_i,y_i)Q_{\theta }^{\infty }(u_{(i-1)},u_i)d{u}_{1:k}}{{\int }_{{\mathbb{R}}^{kd}\times \Theta }{\prod }_{i=1}^k{G}_{\theta }(u_i,y_i)Q_{\theta }^{\infty }(u_{(i-1)},u_i){\pi }_{\theta }(d\theta )d{u}_{1:k}}$$

$$\frac{1}{N}\sum _{i=1}^N\phi (u^i)$$

$$U^{\prime }=u^i+Z,$$

Full text

Key words: Multilevel Monte Carlo, Markov chain Monte Carlo, Sequential Monte Carlo, Ensemble Kalman filter, Coupling.

Advanced Multilevel Monte Carlo Methods

BY AJAY JASRA1, KODY LAW2, & CARINA SUCIU3

1Department of Statistics & Applied Probability, National University of Singapore, Singapore, 117546, SG. E-Mail: staja@nus.edu.sg

2Computer Science and Mathematics Division, Oak Ridge National Laboratory, Oak Ridge, 37934 TN, USA. E-Mail: kodylaw@gmail.com

3Center for Uncertainty Quantification in Computational Science & Engineering, King Abdullah University of Science and Technology, Thuwal, 23955-6900, KSA.

E-Mail: oana.suciu@kaust.edu.sa

1 Introduction

Let $(\mathsf{E},\mathcal{E})$ be a measurable space, $\pi$ be a probability measure on $(\mathsf{E},\mathcal{E})$ and $\varphi:\mathsf{E}\rightarrow\mathbb{R}$ be a measurable and $\pi-$integrable function. In this article we are concerned with the computation of

[TABLE]

for many different $\pi-$integrable functions $\varphi$. In addition, if $\pi$ admits a density w.r.t. a dominating $\sigma-$finite measure $du$ and if one can write for $\kappa:\mathsf{E}\rightarrow\mathbb{R}_{+}$

[TABLE]

where $Z$ is not known, but one can obtain $\kappa$ up-to a non-negative unbiased estimator, then one is also interested in the computation of $Z$. These problems occur in a wide variety of real applications, often concerning Bayesian statistical inference. For instance, the computation of (1) can be associated to posterior expectations, and the value of $Z$ can be used for Bayesian model selection; see  [76]. Many of these problems are found in many real applications, such as meteorology, finance and engineering; see [65, 76]. Later in this article, we will expand upon the basic problem here.

We focus on the case when $\pi$ is associated to some complex continuum problem, for instance, a continuous-time stochastic process or the solution of a partial differential equation (PDE), although the methodology described in this article is not constrained to such examples. Also, we will assume that:

One must resort to numerical methods to approximate (1) or $Z$. 2. 2.

One can, at best, hope to approximate expectations w.r.t. some biased version of $\pi$, call it $\pi_{L}$. It is explicitly assumed that this bias is associated with a scalar parameter $h_{L}\in\mathbb{R}_{+}$ and that the bias disappears as $h_{L}\rightarrow 0$. 3. 3.

When using Monte Carlo methods, exact sampling from $\pi_{L}$ is not possible; that is, one cannot sample i.i.d. from $\pi_{L}$.

Examples of models satisfying 1. & 2. include laws of stochastic differential equations (SDE) for which one cannot sample exactly (e.g.,  [64]), and one resorts to Euler or Milstein discretization. In addition, the law of a quantity of interest (QOI) resulting from the solution of a PDE associated to random input parameters, which cannot be solved exactly and needs to be numerically approximated. Examples satisfying 1.-3. include for example Bayesian instances of the above, where one updates the prior probability distribution based on noisy data to obtain the posterior conditional on the observed data (e.g.,  [48]) or general models where approximate Bayesian computation (e.g.,  [67]) must be used.

1.1 Monte Carlo and Multilevel Monte Carlo

For now, let us assume that only 1. & 2. apply. In this context, one could, for instance, sample $U^{1},\dots,U^{N}$ i.i.d. from $\pi_{L}$ and one could use the standard Monte Carlo estimator $\frac{1}{N}\sum_{i=1}^{N}\varphi(u^{i})$, which approximates $\mathbb{E}_{\pi_{L}}[\varphi(U)]=\int_{\mathsf{E}}\varphi(u)\pi_{L}(du)$. To explain what will follow, we will suppose the following.

Assumption 1.1** (Cost and discretization error).**

There are $\alpha,\zeta>0$ such that

• •

The cost of simulating one sample is $\mathcal{O}(h_{L}^{-\zeta})$.

• •

The bias is of order $\mathcal{O}(h_{L}^{\alpha})$.

This scenario would occur under an Euler discretization of a suitably regular SDE and for appropriate $\varphi$, with $\alpha=\zeta=1$. For simplicity suppose, $h_{L}=2^{-L}$. Consider the mean square error (MSE) associated to the Monte Carlo estimator

[TABLE]

where $\mathbb{E}[\cdot]$ is the expectation operator w.r.t. the distribution of the samples $(U^{1},\dots,U^{N})$. Note that $\mathbb{E}[\varphi(U^{i})]=\mathbb{E}_{\pi_{L}}[\varphi(U)]$. Adding and subtracting $\mathbb{E}_{\pi_{L}}[\varphi(U)]$, and assuming $\varphi$ has a second moment w.r.t. $\pi_{L}$, one has that the MSE is equal to

[TABLE]

which is the standard variance plus bias squared. Now let $1>\epsilon>0$ be given, and suppose one wants to control the MSE, e.g., so that it is $\mathcal{O}(\epsilon^{2})$. One begin by controlling the bias, by setting $L$. The constraint that $2^{-2\alpha L}=\mathcal{O}(\epsilon^{2})$ can be satisfied by choosing $L\propto-\frac{\log(\epsilon)}{\alpha\log(2)}$. Then the constraint that the variance is $\mathcal{O}(\epsilon^{2})$ can be satisfied by choosing $N\propto\epsilon^{-2}$. The cost can then be controlled by $2^{\zeta L}\epsilon^{-2}=\mathcal{O}(\epsilon^{-2-\frac{\zeta}{\alpha}})$. In the case of an Euler discretization of a sufficiently regular SDE, one can asymptotically obtain an MSE of $\mathcal{O}(\epsilon^{2})$ for a cost of $\mathcal{O}(\epsilon^{-3})$.

The multilevel Monte Carlo (MLMC) method [39, 40, 46] is designed to improve over the cost of Monte Carlo. As above, suppose $h_{l}=2^{-l}$. The idea is to consider a hierarchy, $\infty>h_{1}>\cdots>h_{L}>0$ and consider the respresentation

[TABLE]

where for $1\leq l\leq L$, $\mathbb{E}_{\pi_{l}}$ is the expectation w.r.t. $\pi_{l}$ (i.e., the biased approximation with parameter $h_{l}$) and for $l=1$, $\mathbb{E}_{\pi_{l-1}}[\varphi(U)]:=0$. This will be referred to in what follows as the ML identity. Here, it is assumed that for each probability $\pi_{l}$ one is only interested in a marginal on $\mathsf{E}$, even if the entire space must be enlarged to facilitate the biased approximation. So, for instance, $\pi_{l}$ may be defined on a larger space than $\mathsf{E}$, but it admits a marginal on $\mathsf{E}$ which approaches $\pi$ as $l$ grows. Furthermore, it is explicitly assumed that the cost of sampling or evaluating $\pi_{l}$ grows with $l$; again this would occur in most applications mentioned above. To approximate the first term in the summation of (3), one samples $U^{1}(1),\dots,U^{N_{1}}(1)$ i.i.d. from $\pi_{1}$ and one uses the standard Monte Carlo estimator $\frac{1}{N}\sum_{i=1}^{N}\varphi(u^{i}(1))$. For the remainder of the terms $2\leq l\leq L$, we suppose that it is possible to sample a (dependent) coupling of $(\pi_{l},\pi_{l-1})$ with samples $(U_{l}(l),U_{l-1}(l))$ such that the following holds.

Assumption 1.2** (Variance).**

There is a $\beta>0$ such that the variance w.r.t. the coupling of $(\pi_{l},\pi_{l-1})$,

[TABLE]

Note that by coupling, we mean that $U_{l}(l)\sim\pi_{l}$ and $U_{l-1}(l)\sim\pi_{l-1}$ (the random variables are generally dependent). Couplings which satisfy the above bullet exist, for instance, in the context of SDE. If the SDE is suitably regular and is approximated by Euler method, then $\beta=1$, or $\beta=2$ if the diffusion coefficient is constant. In order to approximate the summands in (3) for $2\leq l\leq L$, draw $N_{l}$ i.i.d. samples $(U_{l}^{1}(l),U_{l-1}^{1}(l)),\dots,(U_{l}^{N_{l}}(l),U_{l-1}^{N_{l}}(l))$ from the coupling $(\pi_{l},\pi_{l-1})$, and use the unbiased estimator

[TABLE]

The multilevel estimator is thus

[TABLE]

One can analyze the MSE as above. It is equal to

[TABLE]

and the associated cost is $\sum_{l=1}N_{l}h_{l}^{-\zeta}$, where we assume that the cost of sampling the coupling $(\pi_{l},\pi_{l-1})$ is at most the cost of sampling $\pi_{l}$. Since we have assumed $h_{1}=\mathcal{O}(1)$, then $\frac{1}{N_{1}}\mathbb{V}\textrm{ar}_{\pi_{1}}[\varphi(U)]\leq\frac{C}{N_{1}}$, for $\infty>C>0$ a constant independent of $l$. Now let $1>\epsilon>0$ be given, and suppose one wants to control the MSE, e.g.,  so that it is $\mathcal{O}(\epsilon^{2})$. One controls the bias as above by letting

[TABLE]

Then one seeks to minimize the cost $\sum_{l=1}N_{l}h_{l}^{-\zeta}$ in terms of $N_{1},\dots,N_{L}$, subject to the constraint

[TABLE]

This constrained optimization problem is solved in [39] and has the solution $N_{l}\propto h_{l}^{(\beta+\zeta)/2}$ to obtain a MSE of $\mathcal{O}(\epsilon^{2})$. Solving for the Lagrange multiplier, with equality above, one has that

[TABLE]

where $K_{L}=\sum_{l=1}^{L}h_{l}^{(\beta-\zeta)/2}$. Note that $K_{L}$ may depend upon $L$, depending upon the values of $\beta,\zeta$. In the Euler case $\beta=\zeta$. So, one is able to obtain an MSE of $\mathcal{O}(\epsilon^{2})$ for the cost $\mathcal{O}(\epsilon^{-2}\log(\epsilon)^{2})$. In the special case in which the diffusion coefficient is constant, one obtains the Milstein method with $\beta>\zeta$, so the cost can be controlled by $\mathcal{O}(\epsilon^{-2})$.

The MLMC framework discussed above is considered in various different guises in this paper. The cost will always scale as in Assumption 1.1, and some analogue of the bias from Assumption 1.1 will determine $L$ as defined in (4). We will then require rates on different quantities analogous to Assumption 1.2, in order to ensure the choice of $N_{l}$ in (5) is optimal.

1.2 Methodology Reviewed

In the above section, we have supposed only the points 1. and 2. However, in this article we are considering all three points. In other words, it is not possible to exactly sample from any of $\pi_{L}$ or $\pi_{1},\dots,\pi_{L-1}$. One of the critical ingredients of the MLMC method is sampling dependent couples of the pairs $(\pi_{l},\pi_{l-1})$, which one might argue is even more challenging than sampling from $\pi_{L}$ for a single given $L$. In the context of interest, one might use Markov chain Monte Carlo (MCMC – see e.g., [76, 78]) or Sequential Monte Carlo (SMC – see e.g., [23, 24, 32, 31]) to overcome the challenges of not being able to sample $\pi_{L}$. However, a simple procedure of trying to approximate the ML identify (3) by sampling independent MCMC chains targeting $\pi_{1},\dots,\pi_{L}$ would seldom lead to improvement over just sampling from $\pi_{L}$. So, in such contexts where also using the MLMC approach makes sense, the main issue is how can one utilize such methodology so that one reduces the cost relative to exact sampling from $\pi_{L}$, for a given MSE. There have been many works on this topic and the objective of this article is to review these ideas as well as to identify important areas which could be investigated in the future.

The challenge lies not only in the design and application of the method, but in the subsequent analysis of the method, i.e., verifying that indeed it yields an improvement in cost for a given level of MSE. For instance, the analysis of MCMC and SMC rely upon techniques in Markov chains (e.g., [68, 78]) and Feynman-Kac formulae (e.g., [23, 24]). We highlight these techniques during our review.

1.3 Structure

This article is structured as follows. In Section 2, we give a collection of motivating examples from applied mathematics and statistics, for which the application of multilevel methods would make sense, and are of interest from a practical perspective. These examples are sufficiently complex that standard independent sampling is not currently possible, but advanced simulation methods such as those described in the previous subsection can be used. In Section 3, a short review of some of the computational methods for which this review is focussed on is given. In Section 4, we review several methods which have been adopted in the literature to date, mentioning the benefits and drawbacks of each approach. In Section 5, some discussion of the potential for future work is provided.

We end this introduction by mentioning that this review is not intended to be comprehensive. For instance, we do not discuss quasi-Monte Carlo methods or debiasing methods (e.g., [75]). An effort is of course made to discuss as much work as possible that exists under the umbrella of advanced MLMC methods.

2 Motivating Examples

2.1 Bayesian Inverse Problems

We consider the following example as it is described in [8] (see also [48] and the references therein).

We introduce the nested spaces $V:=H^{1}(\Omega)\subset L^{2}(\Omega)\subset H^{-1}(\Omega)=:V^{*}$, where the domain $\Omega$ will be defined later. Furthermore, denote by $\langle\cdot,\cdot\rangle,\|\cdot\|$ the inner product and norm on $L^{2}$, and by $\langle\cdot,\cdot\rangle,|\cdot|$ the finite dimensional Euclidean inner product and norms. Denote weighted norms by adding a subscript as $\langle\cdot,\cdot\rangle_{A}:=\langle A^{-\frac{1}{2}}\cdot,A^{-\frac{1}{2}}\cdot\rangle$, with corresponding norms $|\cdot|_{A}$ or $\|\cdot\|_{A}$ for Euclidean and $L^{2}$ spaces, respectively (for symmetric, positive definite $A$ with $A^{\frac{1}{2}}$ being the unique symmetric square root).

Let $\Omega\subset\mathbb{R}^{D}$ with $\partial\Omega\in C^{1}$ convex. For $f\in V^{*}$, consider the following PDE on $\Omega$:

[TABLE]

where

[TABLE]

Define $u=\{u_{k}\}_{k=1}^{K}$, with $u_{k}\sim\mathcal{U}[-1,1]$ i.i.d. (the uniform distribution on $[-1,1]$). This determines the prior distribution for $u$. Assume that $\bar{u},\phi_{k}\in C^{\infty}$ for all $k$ and that $\|\phi_{k}\|_{\infty}=1$. In particular, assume $\{\sigma_{k}\}_{k=1}^{K}$ decay with $k$. The state space is $\mathsf{E}=\prod_{k=1}^{K}[-1,1]$. Assume the following property holds: $\inf_{x}\hat{u}(x)\geq\inf_{x}\bar{u}(x)-\sum_{k=1}^{K}\sigma_{k}\geq u_{*}>0$ so that the operator on the left-hand side of (6) is uniformly elliptic. Let $p(\cdot;u)$ denote the weak solution of (6) for parameter value $u$. Define the following vector-valued function

[TABLE]

where $g_{m}$ are elements of the dual space $V^{*}$ for $m=1,\ldots,M$. It is assumed that the data take the form

[TABLE]

where $\mathcal{N}(0,\Gamma)$ denotes the Gaussian random variable with mean [math] and covariance $\Gamma$, and $\perp$ denotes independence. The unnormalized density for $u\in\mathsf{E}$ is then is given by:

[TABLE]

2.1.1 Approximation

Consider the triangulated domains (with sufficiently regular triangles) $\{\Omega^{l}\}_{l=1}^{\infty}$ approximating $\Omega$, where $l$ indexes the number of nodes $d_{l}\propto h_{l}^{-D}$, for triangulation diameter $h_{l}$, so that we have $\Omega^{1}\subset\cdots\subset\Omega^{l}\subset\Omega^{\infty}:=\Omega$. Furthermore, consider a finite element discretization on $\Omega^{l}$ consisting of $H^{1}$ functions $\{\psi_{\ell}\}_{\ell=1}^{d_{l}}$. Denote the corresponding space of functions of the form $\varphi=\sum_{\ell=1}^{d_{l}}v_{\ell}\psi^{l}_{\ell}$ by $V^{l}$, and notice that $V^{1}\subset V^{2}\subset\cdots\subset V^{l}\subset V$. By making the further Assumption 7 of [48] that the weak solution $p(\cdot;u)$ of (6)-(7) for parameter value $u$ is in the space $W=H^{2}\cap H^{1}_{0}\subset V$, one obtains a well-defined finite element approximation $p^{l}(\cdot;u)$ of $p(\cdot;u)$, with a rate of convergence in $V$ or $L^{2}$, independently of $u$. Thus, the sequence of distributions of interest in this context is:

[TABLE]

One is also interested in computing $Z_{L}$, for instance, to perform model selection or averaging.

Exact sampling of this sequence of posterior distributions is not possible in general, and one must resort to an advanced method such as MCMC. But it is not obvious how one can leverage the MLMC approach for this application. Several strategies are suggested later on in the article.

2.2 Partially Observed Diffusions

The following model is considered, as described in [54, 59]. Consider the partially-observed diffusion process:

[TABLE]

with $U_{t}\in\mathbb{R}^{d}$, $t\geq 0$, $U_{0}$ given $a:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$ (denote the $j^{th}-$element as $a^{j}(U_{t})$), $b:\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times d}$ (denote the $j^{th},k^{th}-$element as $b^{j,k}(U_{t})$) and $\{W_{t}\}_{t\in[0,T]}$ a Brownian motion of $d-$dimensions. Some assumptions are made in [54, 59] to ensure that the diffusion has an appropriate solution; see [54, 59] for details.

It will be assumed that the data are regularly spaced (i.e., in discrete time) observations $y_{1},\dots,y_{n}$, with $y_{k}\in\mathbb{R}^{m}$. It is assumed that conditional on $U_{k\delta}=u_{k\delta}$, for $1\geq\delta>0$, $Y_{k}$ is independent of all other random variables and has density $G(u_{k\delta},y_{k})$. For simplicity of notation, let $\delta=1$ (which can always be done by rescaling time), so $U_{k}=U_{k\delta}$. The joint probability density of the observations and the unobserved diffusion at the observation times is then

[TABLE]

where $Q^{\infty}(u_{(i-1)},u)$ is the transition density of the diffusion process as a function of $u$, i.e., the density of the solution $U_{1}$ of Eq. (10) at time $1$ given initial condition $U_{0}=u_{(i-1)}$.

In this problem, one wants to sequentially approximate a probability on a fixed space. For $k\in\{1,\dots,n\}$, the objective is to approximate the filter

[TABLE]

with $u_{1:k}=(u_{1},\dots,u_{k})$ and $y_{1:k}=(y_{1},\dots,y_{k})$. The shorthand notation $\pi^{k}(\cdot)=\pi(\cdot|y_{1:k})$ is used above and in what follows. Note that we will use $\pi_{{\infty}}$ as the notation for measure and density, with the use clear from the context. It is also of interest, to estimate the normalizing constant, or marginal likelihood

[TABLE]

Note that the filtering problem has many applications in engineering, statistics, finance, and physics (e.g., [12, 21, 31] and the references therein)

2.2.1 Approximation

There are several issues associated to the approximation of the filter and marginal likelihood, sequentially in time. Even if one knows $Q^{\infty}$ pointwise, up-to a non-negative unbiased estimator, and/or can sample exactly from the associated law, advanced computational methods, such as particle filters (e.g., [32, 37]) – an exchangeable term for SMC when used in a filtering context – are often adopted in order to estimate the filter. In the setting considered in this paper, it is assumed that one cannot

• •

evaluate $Q^{\infty}$ pointwise, up-to a non-negative unbiased estimator ;

• •

sample from the associated distribution of $Q^{\infty}$.

$Q^{\infty}$ and its distribution must be approximated by some discrete time-stepping method [64] (for time-step $h_{l}=2^{-l}$).

For simplicity and illustration, Euler’s method [64] will be considered. One has

[TABLE]

for $m=0,\dots,k_{l}-1$, where $k_{l}=2^{l}$ and $\mathcal{N}_{d}(0,I_{d})$ is the $d-$dimensional normal distribution with mean zero and covariance the identity (when $d=1$ we omit the subscript). Here $U_{k}^{l}(k_{l})=U^{l}_{k}$, $U^{l}_{k}(0)=U^{l}_{k-1}=U_{k-1}^{l}(k_{l})$. The numerical scheme gives rise to its own transition density between observation times $Q^{l}(u^{l}_{k-1},u^{l}_{k})$.

Therefore, one wants to approximate for $k\in\{1,\dots,n\}$ the filter

[TABLE]

and marginal likelihood

[TABLE]

First, we consider how this task can be performed using SMC and how that in turn can be extended to the MLMC context.

2.2.2 Parameter Estimation

Suppose that there is a static parameter $\theta\in\Theta\subseteq\mathbb{R}^{d_{\theta}}$ in the model, so

[TABLE]

and $G_{\theta}$ is the likelihood function above. If one assumes a prior $\pi_{\theta}$ on $\theta$, then one might be interested in, for $k$ fixed:

[TABLE]

and the associated discretization. We consider how the latter task is possible using MCMC for a given level and how that is extended for the MLMC case.

3 Some Computational Methods

The following section gives a basic introduction to some computational methods that can be used for the examples of the previous section. The review is quite basic, in the sense that there are numerous extensions in the literature, but, we try to provide the basic ideas, with pointers to the literature, where relevant. We first consider basic MCMC in Section 3.1. Then, in Section 3.2, SMC is discussed in context of particle filtering. Here, a certain SMC algorithm, SMC samplers (e.g., [25]) is approached, which uses MCMC algorithms within it. In Section 3.3, we discuss particle MCMC [1], which combines SMC within MCMC proposals. Finally, in Section 3.4, we discuss ensemble Kalman filter approaches.

3.1 Markov chain Monte Carlo

We consider a target probability $\pi_{L}$ on measurable space $(\mathsf{E},\mathcal{E})=(\mathbb{R}^{d},\mathcal{B}(\mathbb{R}^{d}))$, with $\mathcal{B}(\mathbb{R}^{d})$ the Borel sets on $\mathbb{R}^{d}$; extensions to other spaces is simple, but it is omitted for brevity. We write the target density of $\pi_{L}$ w.r.t. Lebesgue measure as $\pi_{L}$ also, following standard practice.

The idea of MCMC is to build an ergodic Markov kernel of invariant measure $\pi_{L}$ or, at least, that $\pi_{L}$ is a marginal of the invariant measure of the Markov chain – we concentrate on the former case. That is, samples of the Markov chain $U^{1},\dots,U^{N}$ have the property that for $\varphi:\mathsf{E}\rightarrow\mathbb{R}$, $\varphi$ $\pi-$integrable, one has that estimates of the form

[TABLE]

will converge almost surely as $N\rightarrow\infty$ to $\mathbb{E}_{\pi_{L}}[\varphi(U)]$.

There are many ways to produce the Markov chain. Here we will only describe the standard random walk Metropolis-Hastings algorithm. At the end of the description, references for more recent work are given. Suppose at time $i$ of the algorithm, $u^{i}$ is the current state of the Markov chain. A new candidate state $U^{\prime}|u^{i}$ is proposed according to

[TABLE]

where $Z\sim\mathcal{N}_{d}(0,\Sigma)$ independently of all other random variables. The proposed value is accepted ($u^{i+1}=u^{i}$) with probability:

[TABLE]

otherwise it is rejected and $u^{i+1}=u^{i}$. The scaling of the proposal, i.e., the proposal covariance $\Sigma$, is often chosen so that the acceptance rate is about 0.234 [79], although there are adaptive methods for doing this; see [2, 44].

The algorithm mentioned here is the most simple approach. The ideas can be extended to alternative proposals, Langevin [77], Hamiltonian Monte Carlo ([34], see also [70]), pre-conditioned Crank Nicholson (e.g., [20] and the references therein). The algorithm can also be used in infinite dimensions (e.g., Hilbert spaces [20]) and each dimension need not be updated simultaneously - for example, one can use Gibbs and Metropolis-within-Gibbs approaches; see [76, 78] for some coverage. There are also population-based methods (e.g., [53] and the references therein) and non-reversible approaches (e.g., [72, 73, 11] and the references therein). Even this list is not complete: the literature is so vast, that one would require a book length introduction, which is well out of the scope of this work - the reader is directed to the previous papers and the references therein for more information. There is also a well established convergence theory; see [68, 78] for information. We remark also that the cost of such MCMC algorithms can be quite reasonable if $d$ is large, often of polynomial order (e.g., [79]) and can be dimension free when there is a well-defined limit as $d$ grows [9, 20]. Note that finally, it is not simple to use ‘standard’ MCMC to estimate normalizing constants.

3.2 Sequential Monte Carlo

Here we will consider the standard SMC algorithm for filtering, also called the particle filter in the literature. To assist our discussion, we focus on the filtering density from Section 2.2.1

[TABLE]

Here there are no static parameters to be estimated. To facilitate the discussion, we will suppose that $Q^{L}$ can be evaluated pointwise, although of course, it is generally not the case in our application. Moreover, it will be assumed that $G(u_{i},y_{i})Q^{L}(u_{(i-1)},u_{i})>0$ for each $i\geq 1$, $u_{i-1},u_{i}$.

We suppose we have access to a collection of proposals $q_{1}(u_{1}),q_{2}(u_{1},u_{2}),q_{3}(u_{2},u_{3}),\dots$, where $q_{j}(u_{j-1},u_{j})$ is a positive probability density in $u_{j}$, for each value of $u_{j-1}$. The particle filter is then as follows:

• •

Initialize. Set $k=1$, for $i\in\{1,\dots,N\}$ sample $u_{1}^{i}$ from $q_{1}$ and evaluate the weight

[TABLE]

• •

Iterate: Set $k=k+1$,

• –

Resample $(\hat{u}_{k-1}^{1},\dots,\hat{u}_{k-1}^{N})$ according to the weights $(w_{k-1}^{1},\dots,w_{k-1}^{N})$.

• –

Sample $u_{k}^{i}|\hat{u}_{k-1}^{i}$ from $q_{k}$, for $i\in\{1,\dots,N\}$, and evaluate the weight

[TABLE]

The resampling step can be performed using a variety of schemes, such as systematic, multinomial residual etc; the reader is referred to [32] for more details. For $\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$, and $\varphi$ $\pi^{k}_{{L}}-$integrable, one has the consistent estimate:

[TABLE]

In addition, the marginal likelihood is unbiasedly [23] estimated by

[TABLE]

with the abuse of notation that $\hat{u}_{0}^{i}=u_{0}$. In principle, for $\varphi:\mathbb{R}^{kd}\rightarrow\mathbb{R}$, and $\varphi$ $\pi^{k}_{{L}}-$integrable, one could also try to estimate $\mathbb{E}_{\pi^{k}_{{L}}}[\varphi(U_{1:k})]$ but this does not work well in practice due to the well-known path degeneracy problem; see [32, 62].

The algorithm given here is one of the most basic and many modifications can enhance the performance of this algorithm; see [23, 24, 32, 62] for some ideas. The theoretical validity of the method has been established in many works; see e.g., [13, 17, 23, 24, 30]. The algorithm performs very well w.r.t. the time parameter $k$. Indeed $L_{p}-$errors for estimates such as (13) are $\mathcal{O}(N^{-1/2})$ where the constant is independent of time and the relative variance of (14) is $\mathcal{O}(k/N)$ (if $N>Ck$ for $C$ some constant independent of $k$); see [24] and the references therein. One of the main issues with particle filters/SMC methods in this context is that they do not perform well in high dimensions (i.e., for large $d$) often having an exponential cost in $d$ [86]. Note however, that if there is a well-defined limit as $d$ grows, SMC methods can be designed to perform quite well in practice on a finite time horizon (see e.g., [61] and the next section). There have been some methods developed for high-dimensional filtering (e.g., [7, 69, 74]), however, they are only useful for a small class of models.

3.2.1 Sequential Monte Carlo Samplers

Consider a sequence of distributions $\pi_{1},\dots,\pi_{L}$ on a common measurable space. In addition to this suppose we have Markov kernels $M_{2},\dots,M_{L}$ of invariant measures $\pi_{2},\dots,\pi_{L}$. This is possible if the densities are known up-to a constant (and potentially also a non-negative unbiased estimator - although this is not considered at the moment), simply by using using MCMC. The SMC sampler algorithm (e.g., [25]) can be used to approximate expectations w.r.t. $\pi_{1},\dots,\pi_{L}$, as well as to estimate ratios of normalizing constants. The un-normalized densities (assumed to exist w.r.t. a common dominating measure) of $\pi_{1},\dots,\pi_{L}$ are written $\kappa_{1},\dots,\kappa_{L}$. To ease the notational burden, we suppose one can sample from $\pi_{1}$, but this is not necessary.

The algorithm is as follows:

• •

Initialize. Set $l=1$, for $i\in\{1,\dots,N\}$ sample $u_{1}^{i}$ from $\pi_{1}$.

• •

Iterate: Set $l=l+1$. If $l=L+1$ stop.

• –

Resample $(\hat{u}_{l-1}^{1},\dots,\hat{u}_{l-1}^{N})$ using the weights $(w_{l}^{1},\dots,w_{l}^{N})$ where, for $i\in\{1,\dots,N\}$,

[TABLE]

• –

Sample $u_{l}^{i}|\hat{u}_{l-1}^{i}$ from $M_{l}$ for $i\in\{1,\dots,N\}$.

One can estimate expectations w.r.t. $\pi_{l}$, for $\varphi:\mathsf{E}\rightarrow\mathbb{R}$, $\varphi$ $\pi_{l}-$integrable. The consistent estimator

[TABLE]

converges almost surely as $N\rightarrow\infty$. In addition, for any $l\geq 2$, we have the unbiased estimator

[TABLE]

which converges almost surely as $N\rightarrow\infty$.

The basic algorithm goes back to at least [51]. Several versions are found in [16, 71], with a unifying framework in [25] and a rediscovery in [15]. Subsequently, several refined and improved versions of the algorithm have appeared [19, 26, 47, 58, 81], including those which allow algorithmic parameters to be set adaptively, that is, without user specification.

Contrary to particle filters, when $\mathsf{E}=\mathbb{R}^{d}$, this method indeed performs quite well w.r.t. the dimension $d$ with only polynomial cost in $d$; see [3, 4]. Whilst the underlying theory for this algorithm is very similar to particle filters and it is covered in [23, 24], there are some additional results in [5, 82, 90]. In particular, [5] establish that when one updates parameters adaptively, such as in [58, 81], then the algorithm is still theoretically correct. The method is very useful in the following scenaria: (i) if one wishes to compute ratios of normalizing constants, (ii) the available MCMC kernels do not mix particularly well, and/or (iii) the target is multimodal and the modes are separated by regions of very low probability.

3.3 Particle Markov chain Monte Carlo

We now consider the scenario of Section 2.2.2. In this context, the standard approach is to consider the extended target with density

[TABLE]

Sampling this distribution is notoriously challenging. One recent method and it’s extensions can be considered the gold standard as we describe now. Note that the SMC algorithm (in Section 3.2) used below uses the Euler discretized dynamics as the proposal.

The particle marginal Metropolis-Hastings (PMMH) algorithm of [1] proceeds as follows.

• •

Initialize. Set $i=0$ and sample $\theta^{0}$ from the prior. Given $\theta^{0}$ run the SMC algorithm in Section 3.2 and record the estimate of $\widehat{Z}_{L,\theta^{0}}^{k}$ from eq. (14).

• •

Iterate:

• –

Set $i=i+1$ and propose $\theta^{\prime}$ given $\theta^{i-1}$ from a proposal $r(\theta^{i-1},\cdot)$.

• –

Given $\theta^{\prime}$ run the SMC algorithm in Section 3.2 and record the estimate $\widehat{Z}_{L,\theta^{\prime}}^{k}$.

• –

Set $\theta^{i}=\theta^{\prime}$ with probability

[TABLE]

otherwise $\theta^{i}=\theta^{i-1}$.

The samples of this algorithm can be used to estimate expectations such as $\int_{\Theta}\varphi(\theta)\pi_{{L}}(\theta|y_{1:k})d\theta$ with

[TABLE]

Note that this is consistent as $N$ grows, in the sense that it recovers the true expectation with probability 1, under minimal conditions. The algorithm can also be extended to allow estimation of the hidden states $u_{1:k}$ as well. There are many parameters of the algorithm, such as the number of samples of the SMC algorithm, and tuning them has been discussed in [1, 33], for example.

The PMMH algorithm is the most basic in [1]. Several enhancements are in [1] and numerous algorithms that improve upon this method can be found in [29, 85].

3.4 Ensemble Kalman filter

The idea of the ensemble Kalman filter (EnKF) is to approximate the filtering distribution (12) using an ensemble of particles and their sample covariance [35]. As such they are sometimes also referred to as sequential Monte Carlo methods, but we believe it is important to distinguish them from the methods described in subsection 3.2. The observations are incorporated as though the process were linear and Gaussian, hence requiring only the covariance approximation. Hence, the method is consistent only in the case of a linear Gaussian model [65]. However, it is robust even in high dimensions [36] and can be tuned to perform reasonably well in tracking and forecasting. It has therefore become very popular among practitioners.

It will be assumed here for simplicity that the observation selection function is given by:

[TABLE]

where $H$ is a linear operator. The linearity assumption is without any loss of generality since if $h$ is nonlinear, one can always extend the system to $(u,v)^{\top}$, where $v=h(u)$.

The EnKF is executed in a variety of ways and only one will be considered here, the perturbed observation EnKF:

[TABLE]

[TABLE]

Here $\xi_{j}^{(n)}$ are i.i.d. draws from $N(0,\Gamma)$. Perturbed observation refers to the fact that each particle sees an observation perturbed by an independent draw from $N(0,\Gamma)$. This procedure ensures the Kalman Filter is obtained in the limit of infinite ensemble in the linear Gaussian case [65]. Notice that the ensemble is not prescribed to be Gaussian, even though it is updated as though it were, so the limiting target is some non-Gaussian $\widehat{\pi}_{L}^{k}$, which is in general not equal to the density defined by (12) (see e.g.,  [66]).

4 Approaches for MLMC Estimation

We now consider various ways in which the MLMC method can be used in these challenging situations, where it is non-trivial to construct couplings of the targets.

4.1 Importance Sampling

In this case, we investigate the ML identity where the sequence of targets $\pi_{1},\dots,\pi_{L}$ are defined on a common measurable space and are known up-to a normalizing constant; i.e., $\pi_{l}(u)=\kappa_{l}(u)/Z_{l}$ as in Section 2.1 and 3.2.1.

In this scenario [8] (see also [6, 27, 28, 52]) investigate the simple modification

[TABLE]

The idea here is simple. If one does not know how to construct a coupling of the targets, then one replaces coupling by importance sampling. The key point is that as the targets $\pi_{l}$ and $\pi_{l-1}$ are very closely related by construction, and therefore the change of measure formula above should facilitate an importance sampling procedure that performs well. Just as for ‘standard’ MLMC (for instance as described in Section 1.1) where the coupling has to be ‘good enough’, the change of measure needs to be chosen appropriately to ensure that this approach can work well. Recall from Section 3.2.1 that SMC samplers can be designed to sequentially approximate $\pi_{1},\dots,\pi_{L}$, and the ratios $Z_{l}/Z_{l-1}$. Therefore the change of measure in (16) is very natural here.

The approach in [8] is to simply run the algorithm of Section 3.2.1, except at step $l$ one resamples $N_{l+1}<N_{l}$ particles, where the schedule of numbers $N_{0:L-1}$ is chosen using a similar principle as for standard MLMC. The identity (16) can be approximated via:

[TABLE]

Note that the algorithm need only be run up-to level $L-1$. [8] not only show that this is consistent, but also give a general MLMC theorem using the theory in [23] with some additional work and assumptions (which are relaxed in [28]). In the context of the example of Section 2.1, the authors show that the work to compute expectations relative to standard SMC samplers (as in Section 3.2.1) is reduced to achieve a given MSE, under the following assumptions.

Assumption 4.1** (MLSMC samplers).**

There is a $\varepsilon>0$ such that for all $l=1,\dots,L$, $u,v\in\mathsf{E}$, and $A\in\mathcal{E}$

• •

$\varepsilon<\kappa_{l}(u)<\varepsilon^{-1}$*   ;*

• •

$\varepsilon M_{l}(v,A)<M_{l}(u,A)<\varepsilon^{-1}M_{l}(v,A)$, where we recall from subsection 3.2.1 that $M_{l}$ is the Markov kernel with invariant measure proportional to $\kappa_{l}$, used to mutate the updated population of samples at step $l$.

The main reason why this approach can work well, can be explained by terms that look like

[TABLE]

In the context of the problem in Section 2.1, this term will tend to zero at a rate $h_{l}^{\beta}$, under suitable assumptions and in an appropriate norm, just as the variance terms in the coupling of Section 1.1. In more standard importance sampling language, the weight tends to one as the sequence of target distributions gets more precise. This particular approach exploits this property, and the success of the method is dependent upon it. In particular, the key quantities for which one needs to obtain rates $\alpha$ and $\beta$ for are summarized in table 1, and the convergence rate is illustrated in Figure 1.

The method of [8] has been extended to the computation of normalizing constants [28] and has been applied to other examples, such as non-local equations [52] and approximate Bayesian computation [60]. The method has also been extended to the case that the accuracy of the approximation improves as the dimension of the target grows; see [6].

The importance sampling idea has also been considered in other articles such as [48, 83]. In [48], the change of measure appears in the context of MCMC and is not too dissimilar to the one presented in (16), although it is with respect to the higher level in the couple, and multiple discretization indices are considered as well (see also [45]). In [83], the numerator and denominator of (1), arising from the form (2), are approximated independently in terms of expectation w.r.t. the prior, in a Bayesian set up, and they find reasonably decent performance. This corresponds to a change of measure for Monte Carlo and Quasi Monte Carlo, although for MLMC it is slightly different, since MLMC is used separately for the numerator and the denominator. In general such a procedure is not advisable unless the prior is very informative; one expects that the estimators for each of the numerator and the denominator will have very high variance.

4.2 Approximate Coupling

The case of Section 2.2.2 will be considered here in order to illustrate this idea. In particular, there is a stochastic process that is partially observed. It is assumed that the dynamics of the associated discretized processes can be easily coupled. Given this, an approximate coupling is devised, which: (i) can be sampled from (using MCMC/SMC), and (ii) has marginals which are similar to, but not exactly equal to, the pair $(\pi_{l},\pi_{l-1})$. Then the difference $\{\mathbb{E}_{\pi_{l}}-\mathbb{E}_{\pi_{l-1}}\}[\varphi(U)]$ is replaced with an importance sampling formula (change of measure w.r.t. the approximate coupling) and is approximated by sampling from the approximate coupling. The motivation for the idea will become clear as it is described in more detail. The approach is considered in [55] (see also [56]).

We begin by considering:

[TABLE]

for two levels $l,l-1\geq 1$. We know that one can construct a good coupling of the two discretized kernels $Q_{\theta}^{l},Q_{\theta}^{l-1}$ for any fixed $\theta$ by sampling the finer Gaussian increments and concantenating them for the coarser discretization (e.g., [39] or as written in [54, 59]). More precisly, given $\check{u}_{i-1}=(\underline{u}_{i-1},\overline{u}_{i-1})\in\mathbb{R}^{2d}$ and $\theta\in\Theta$, there is a Markov kernel $\check{Q}_{\theta}^{l,l-1}$ such that for any $A\in\mathcal{B}(\mathbb{R}^{d})$

[TABLE]

and

[TABLE]

Note that under the coupling considered, the discretized processes are not independent.

We consider the joint probability on $\Theta\times\mathbb{R}^{2kd}$:

[TABLE]

for any non-negative function $\check{G}_{\theta}(\check{u}_{i},y_{i})$. Whilst this function can be ‘arbitrary’, up-to some constraints, we set it as

[TABLE]

This will be explained below. Let $\varphi:\Theta\times\mathbb{R}^{2kd}\rightarrow\mathbb{R}$ be $\pi_{l}^{k}$ and $\pi_{{l-1}}^{k}-$integrable. Then we have (supressing the conditioning on $y_{1:k}$ in the expectations)

[TABLE]

[TABLE]

where

[TABLE]

The difference can then be approximating by sampling from $\check{\pi}_{{l-1:l}}^{k}$, e.g., by using the PMMH 111More precisely, one samples from a suitably extended measure, as described in [1]. from Section 3.3. This is done independently for each summand in the ML identity, with the first summand (the coarsest discretization) sampled by PMMH.

We now explain the idea in more detail. The basic idea is that one knows how to construct an exact coupling of the discretizations of the prior, i.e., the stochastic forward dynamics here, and it is natural to leverage this. However, as noted previously, exact couplings of the posterior are not trivial to sample. Instead, one aims to construct a joint probability that should have marginals which are close to, but not exactly equal to, the correct ones. As the coupling is not exact, one must correct for this fact and, as in Section 4.1, use importance sampling. Just as argued in that section, the associated weights of the importance sampling, that is, the terms $(\overline{H}_{\theta},\underline{H}_{\theta})$ should be well behaved in some sense. This can be ensured by choosing the function $\check{G}_{\theta}$ so that the variance of the weights w.r.t. any probability measure will remain bounded uniformly in time. Hence the reason for its selection as 18. [55] are able to prove, under suitable assumptions on the model and PMMH kernel, that the computational effort to estimate a class of expectations is reduced versus a single PMMH algorithm on the finest level, for a given MSE sufficiently small. The reduction in cost is a direct consequence of the prior coupling and well-behaved importance weights, in connection with the ML identity. Note that the results of [55] do not consider the dependence on the time parameter $k$, and that is something that should be addressed. In particular, the required assumptions in this context are given below.

Assumption 4.2** ((PMMH using MLMC)).**

There is a $\varepsilon>0$ and probability $\nu$ over $\Theta\times\mathbb{R}^{2kd}$, such that for all $l=1,\dots,L$, $u\in\mathbb{R}^{d}$, $y\in\mathbb{R}^{m}$, $\theta\in\Theta$, $A\in\sigma(\Theta)$, and any $\varphi:\Theta\times\mathbb{R}^{2kd}$ bounded and Lipschitz, and $w\in\mathsf{W}$ (the space of all auxiliary variables involved in the higher-dimensional chain), the following hold

• •

$\varepsilon<G_{\theta}(u,y)<\varepsilon^{-1}$*   ;*

• •

$Q^{l}(u,v)>\varepsilon$*   ;*

• •

The final Metropolis kernel $K$ on the extended space satisfies $\int_{\mathsf{W}}\varphi(w^{\prime})K(w,dw^{\prime})\geq\varepsilon\int_{\Theta\times\mathbb{R}^{2kd}}\varphi(v)\nu(dv)$,

The quantities for which rates $\alpha$ and $\beta$ need to be obtained in this context are given in table 2, and the results are illustrated in Figure 2 for 2 example SDE of the type introduced in section 2.2.2, each with 2 unknown parameters $\sigma$ and $\theta$.

We end this section by mentioning that the strategy described in this section is not the only one that could be adopted. For instance, the importance sampling approach of the previous section might also be considered. Indeed that may be considered as an extreme example of approximate coupling in which $\check{G}_{\theta}(\check{u}_{i},y_{i})=G_{\theta}(\underline{u}_{i},y_{i})\delta_{\overline{u}_{i},\underline{u}_{i}}.$ There are some reasons why the approximate coupling method described in this section might be preferred, for the example considered here. Firstly, the terms $\kappa_{l}(u)/\kappa_{l-1}(u)$ are not available pointwise for this example; this could possibly be dealt with by random weight ideas (e.g., [38, 80]), but it is still an issue. Secondly, there is a well-designed MCMC algorithm for the target (i.e., a PMMH kernel) and hence one would like to use this, as it is one of the gold standards for such models. If one elects to use an SMC sampler to approximate (19), then PMMH kernels can be used as described in [19]. The algorithm of [19] can also be used for dynamic (sequential) inference on the parameter, in an MLMC context. The same principle as described above can be generalized to MCMC, and this has been done in the work [57].

4.3 Coupling Algorithms

We now consider the case where one seeks to approximate the differences in the ML identity exactly. This is achieved by somehow trying to correlate or couple stochastic algorithms, rather than by constructing any joint coupling, either exact or approximate. We refer to this approach as coupling algorithms and it is explained further below.

4.3.1 Coupling MCMC

We begin by considering the method in [63], which is a Markov chain approach. Consider two probability measures $(\pi_{l},\pi_{l-1})$, where the support of $\pi_{l-1}$ is $\mathsf{E}$ and the support of $\pi_{l}$ is $\mathsf{E}\times\mathsf{U}$. We focus on computing $\mathbb{E}_{\pi_{l}}[\varphi_{l}(U)]-\mathbb{E}_{\pi_{l-1}}[\varphi_{l-1}(U)]$ where $\varphi_{l}:\mathsf{E}\times\mathsf{U}\rightarrow\mathbb{R}$ and $\varphi_{l-1}:\mathsf{E}\rightarrow\mathbb{R}$ are $\pi_{l}$ and $\pi_{l-1}$ integrable.

Suppose we have a current state $(u_{l},u_{l-1})\in\mathsf{E}\times\mathsf{U}\times\mathsf{E}$. The approach consists first by sampling from $\pi_{l-1}$ exactly. Given this sample and the current $u_{l}\in\mathsf{E}$ a new state is proposed from a proposal $r$ and is accepted or rejected according to the standard Metropolis-Hastings method. It is clear that the samples are coupled and seemingly that they have the correct invariant distribution. [63] show that this approach indeed can obtain the advantages of MLMC estimation for some examples. As noted in that article, supposing exact sampling from $\pi_{l-1}$ is feasible, is not realistic and the authors propose a subsampling method to assist with this. See [63] for more details.

Before concluding this section, we mention the recent related works [50, 84, 88]. The work [88] (see also [41]) considers coupled Stochastic Gradient Langevin algorithms for some i.i.d. models in Bayesian statistics. Note that there the levels are from an Euler discretization associated to the algorithm, not the model per-se as described here.

4.3.2 Coupling Particle Filters

We consider the context of Section 2.2.1 of filtering a discretely and partially observed diffusion processes. We describe an approach in [54] (see also [42, 50, 59, 84]). For $l\geq 2$, $\varphi:\mathbb{R}^{d}\rightarrow\mathbb{R}$, $\pi^{k}_{{l}}$ and $\pi^{k}_{{l-1}}-$integrable we consider approximating the difference $\mathbb{E}_{\pi^{k}_{{l}}}[\varphi(U_{k})]-\mathbb{E}_{\pi^{k}_{{l-1}}}[\varphi(U_{k})]$ sequentially in time. Some of the notations of Section 4.2 are also used. The parameter $\theta$ is also dropped from the notations, as it is assumed to be fixed here.

The multilevel particle filter (MLPF) is described as follows. First, for $l=1$, run a particle filter for the coarsest discretization. Now, run the following procedure independently for each $l\geq 2$.

For $i=1,\dots,N_{l}$, draw $(\overline{U}^{l,i}_{1},\underline{U}_{1}^{l,i})\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\check{Q}^{l,l-1}((u_{0},u_{0}),\cdot)$.

Initialize $k=1$. Do

• (i)

For $i=1,\dots,N_{l}$, draw $(\overline{I}_{k}^{l,i},\underline{I}_{k}^{l,i})$ according to the coupled resampling procedure below. Set $k=k+1$.

• (ii)

For $i=1,\dots,N_{l}$, independently draw $(\overline{U}^{l,i}_{k},\underline{U}^{l,i}_{k})|(\overline{u}_{k-1}^{l,\overline{I}_{k}^{l,i}},\underline{u}_{k-1}^{l,\underline{I}_{k}^{l,i}})\sim\check{Q}^{l,l-1}((\overline{u}_{k-1}^{l,\overline{I}_{k}^{l,i}},\underline{u}_{k-1}^{l,\underline{I}_{k}^{l,i}}),\leavevmode\nobreak\ \cdot\leavevmode\nobreak\ )$.

The coupled resampling procedure for the indices $(\overline{I}^{l,i}_{k},\underline{I}^{l,i}_{k})$ is described below. First let

[TABLE]

Now

• a.

with probability $\alpha_{k}^{l}=\sum_{i=1}^{N_{l}}\overline{w}_{k}^{l,i}\wedge\underline{w}_{k}^{l,i}$, draw $\overline{I}^{l,i}_{k}$ according to

[TABLE]

and let $\underline{I}^{l,i}_{k}=\overline{I}^{l,i}_{k}$.

• b.

otherwise, draw $(\overline{I}^{l,i}_{k},\underline{I}^{l,i}_{k})$ independently according to the probabilities

[TABLE]

Note that by using the coupled kernel $\check{Q}^{l,l-1}$, one is sampling from the exact coupling of the discretized process, $(\overline{U}^{l,i}_{k},\underline{U}^{l,i}_{k})$. Now one wants to maintain as much dependence as possible in the resampling, since resampling is necessary in particle filters. The coupled resampling described above maximizes the probability (conditional on the history) that the pair of samples remain coupled (see also [18]).

In the work [54], it is shown that

[TABLE]

consistently approximates $\mathbb{E}_{\pi^{k}_{{l}}}[\varphi(U_{k})]-\mathbb{E}_{\pi^{k}_{{l-1}}}[\varphi(U_{k})]$. The MLPF estimator of $\mathbb{E}_{\pi^{k}_{L}}[\varphi(U_{k})]$ is therefore given by

[TABLE]

In the case of Euler-Maruyama discretization, [54] it is shown that under suitable assumptions and for finite time the standard choice of $L$ and $N_{1:L}$ as in (4) and (5) provides an MSE of $\mathcal{O}(\epsilon^{2})$ for a cost of $\mathcal{O}(\epsilon^{-2.5})$. For a particle filter the cost required is $\mathcal{O}(\epsilon^{-3})$. The theory is not limited to Euler discretizations, but the ultimate bound on the cost will depend on the convergence rate of the numerical method.

Sufficient assumptions in this case are given by

Assumption 4.3** (MLPF).**

There is a $\varepsilon>0$ such that for all $l=1,\dots,L$, $u,v\in\mathbb{R}^{d}$, and $y\in\mathbb{R}^{m}$, the following hold

• •

$\varepsilon<G_{\theta}(u,y)<\varepsilon^{-1}$*   ;*

• •

$Q^{l}(u,v)>\varepsilon$.

The quantities for which rates $\alpha$ and $\beta$ need to be obtained in this context are given in table 2, and the complexity results are illustrated in Figure 3 for some example SDEs.

In [59] the approach and results are extended to the case of marginal likelihood estimation. We note that one drawback of the mathematical results in [54, 59] is that they do not consider the time-parameter. Note also that the method does not the preserve the standard squared strong convergence of the forward numerical method which defines $\check{Q}^{l,l-1}$. A power of $1/2$ is lost in the rate of convergence as a result of the coupled resampling (see (table 3). In [50, 84] the coupled resampling method is improved by using optimal transportation techniques [89]. Also [42] (see also [43]) obtain empirical results which indicate that more favorable convergence rates may be preserved in certain cases by replacing the resampling step with a deterministic linear transformation of the current population, derived from the optimal transportation. However, in [42, 43, 50, 84] there are no mathematical results which support the encouraging empirical results; one expects that this is mainly a technical challenge.

4.3.3 Coupling the EnKF

As discussed in subsection 3.4, the EnKF targets a different distribution than the filtering distribution in general, which has been denoted $\widehat{\pi}^{k}_{L}$. In between updates, this algorithm proceeds similarly to the MLPF of the previous section, propagating pairs of ensembles for each $l\geq 2$. The fundamental difference in the update results in approximations of increments $\mathbb{E}_{\widehat{\pi}^{k}_{{l}}}[\varphi(U_{k})]-\mathbb{E}_{\widehat{\pi}^{k}_{{l-1}}}[\varphi(U_{k})]$. Therefore, the MSE will ultimately depend upon the difference $\mathbb{E}_{\widehat{\pi}^{k}_{{L}}}[\varphi(U_{k})]-\mathbb{E}_{\pi^{k}_{{L}}}[\varphi(U_{k})]$, which includes a Gaussian bias in addition to the discretization bias. In order to preserve the coupling of this algorithm after the update, the sample covariance is approximated using the entire multilevel ensemble in [49], as follows. Recall the functions $G(u_{i},y_{i})$ are assumed to take the form given in (15).

For $l=2,\dots,L$, and $i=1,\dots,N_{l}$, draw $(\overline{U}^{l,i}_{1},\underline{U}_{1}^{l,i})\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\check{Q}^{l,l-1}((u_{0},u_{0}),\cdot)$. And draw $U^{1,i}_{1}\sim Q^{1}(u_{0},\cdot)$.

Initialize $k=1$. Do

• (i)

Compute the MLMC covariance estimator [10] :

[TABLE]

• (ii)

Compute $K_{k}^{\rm ML}={C}_{k}^{\rm ML}H^{T}(H{C}_{+,k}^{\rm ML}H^{T}+\Gamma)^{-1}$, where $C_{+,k}^{\rm ML}$ the positive semi-definite modification of ${C}_{k}^{\rm ML}$.

• (iii)

For $l=2,\dots,L$, and $i=1,\dots,N_{l}$, independently draw $Y_{k}^{l,i}\sim N(y_{k},\Gamma)$, and compute

[TABLE]

and similarly for $\widehat{\underline{u}}_{k}^{l,i}$ and $\widehat{u}_{k}^{1,i}$. Set $k=k+1$.

• (iv)

For $l=2,\dots,L$, and $i=1,\dots,N_{l}$, independently draw $(\overline{U}^{l,i}_{k},\underline{U}_{k}^{l,i})\stackrel{{\scriptstyle\textrm{i.i.d.}}}{{\sim}}\check{Q}^{l,l-1}((\widehat{\overline{u}}_{k-1}^{i},\widehat{\underline{u}}_{k-1}^{i}),\cdot)$. And draw $U^{1,i}_{k}\sim Q^{1}(\widehat{u}_{k-1}^{1,i},\,\cdot\,)$.

A sufficient assumption in this case is given by

Assumption 4.4** (MLEnKF).**

The coefficients of (10) are globally Lipschitz and the initial condition is in $L^{p}$ for all $p\geq 2$.

The quantities for which rates $\alpha$ and $\beta$ need to be obtained are given in table 4, and the complexity results are illustrated in Figure 4 for an example linear SDE of the form in (10). The work [49] established that slightly modified choices of $L$ and $N_{1:L}$ provide MSE at step $k$ of $\mathcal{O}(|\log\epsilon|^{2n}\epsilon^{2})$ for a cost of $\mathcal{O}(\epsilon^{-2}\tilde{K}_{L}^{3/2})$, where $\tilde{K}_{L}^{1/2}=\sum_{l=1}^{L}h_{l}^{(\beta-\zeta)/3}$. However, the numerical results indicate not only a time-independent rate of convergence without logarithmic penalty, but in fact also a time-uniform constant – see Figure 4. Presumably, the penalty on the MSE is mostly a technical hurdle. The recent work [14] has extended this method to spatial processes, for example given by stochastic partial differential equations. This is the context where the EnKF is typically applied, for example in numerical weather prediction.

4.3.4 Discussion

Some examples of coupling algorithms have been reviewed. The main issues with such techniques are (i) coupling the algorithms correctly, so that the coupling is ‘good enough’, and (ii) mathematical analysis of such couplings to prove that they indeed provide a benefit. These challenges are crucial and must be further studied. It was already mentioned that the works [50, 84] consider coupling the pair of particle filters arising in a PMCMC algorithm, and empirical results are promising. However, establishing that indeed (i)-(ii) would occur in practice is not so easy and at least does not appear to have been done in publicly available research. In the MLMC context, the theoretical and numerical results of [55] indicate a loss of a power of $1/2$ in the rate of strong convergence following from the coupled resampling of section 4.3.2, hindering the ultimate cost of the algorithm (see also table 3). However, the rates may be improved by the resampling based on optimal transportation from [50, 84]. Indeed the works [42, 43] numerically observe preservation of the strong rate of convergence using a deterministic transformation based on the optimal transportation coupling, in lieu of resampling.

5 Future Work and Summary

Here we examined some computational approaches to facilitate the application of the MLMC method in challenging examples, where standard (independent) sampling is not currently possible. Some review of the computational methods was provided, although as we have noted it is a large literature that one cannot hope to include a complete summary of all the methodology. We then detailed various approaches one can use to leverage MLMC within these methods.

There are many areas for possible exploration in future work. One strand consists of considering multi-dimension discretizations, such as in [45]. There are a small number of papers on this topic, such as [22, 57], but there seem to be many possible avenues for future work. Another direction consists of a general method for sampling (e.g., by MCMC/SMC) exact (dependent) couplings of the targets in the ML identity. As we have commented, it does not appear to be trivial, but it may be far from impossible. Such a method would be very beneficial, as one could then appeal to existing literature in order to prove complexity results about MLMC and MIMC versions. One final very interesting avenue for future research is exact coupling, using optimal transport and i.i.d. sampling. For some model structures, the ideas of [87] could be very useful.

Acknowledgements

AJ and CS were supported under the KAUST Competitive Research Grants Program-Round 4 (CRG4) project, “Advanced Multi-Level sampling techniques for Bayesian Inverse Problems with applications to subsurface.” KJHL was supported by ORNL LDRD Seed funding grant number 32102582.