Ensemble Kalman filter for multiscale inverse problems

TL;DR

This paper introduces a new ensemble Kalman filter-based algorithm for efficiently solving inverse problems involving multiscale elliptic PDEs, leveraging homogenization and finite element methods.

Contribution

It presents a novel multiscale inverse problem solver that combines numerical homogenization with ensemble Kalman filtering, including convergence analysis and Bayesian interpretation.

Findings

The method accurately recovers oscillatory tensors from measurements.

Convergence of the approximate solution and posterior distribution is established.

Numerical experiments demonstrate computational efficiency and robustness.

Abstract

We present a novel algorithm based on the ensemble Kalman filter to solve inverse problems involving multiscale elliptic partial differential equations. Our method is based on numerical homogenization and finite element discretization and allows to recover a highly oscillatory tensor from measurements of the multiscale solution in a computationally inexpensive manner. The properties of the approximate solution are analysed with respect to the multiscale and discretization parameters, and a convergence result is shown to hold. A reinterpretation of the solution from a Bayesian perspective is provided, and convergence of the approximate conditional posterior distribution is proved with respect to the Wasserstein distance. A numerical experiment validates our methodology, with a particular emphasis on modelling error and computational cost.

Equations

$$\{\begin{array}{rlrl}-\nabla \cdot (A_u^{\epsilon }\nabla p^{\epsilon })& =f,& & \text{ in }\Omega ,\\ p^{\epsilon }& =0,& & \text{ on }\partial \Omega ,\end{array}$$

$$\{\begin{array}{rlrl}-\nabla \cdot (A_u^0\nabla p^0)& =f,& & \text{ in }\Omega ,\\ p^0& =0,& & \text{ on }\partial \Omega ,\end{array}$$

$$\text{find }u\in X\text{ given observations }y=\mathcal{G}(u)+\eta \in Y,$$

$$\Xi (z)=\left[\begin{array}{c}u\\ \mathcal{G}(u)\end{array}\right],\text{ for }z=\left[\begin{array}{c}u\\ v\end{array}\right]\in Z,$$

$$z_{n+1}=\Xi (z_n),n=0,1,\dots$$

$$y_{n+1}=H{z}_{n+1}+{\eta }_{n+1},$$

$$z_0^{(j)}=\left[\begin{array}{c}{\psi }^{(j)}\\ \mathcal{G}({\psi }^{(j)})\end{array}\right],$$

$${\hat{z}}_{n+1}^{(j)}=\Xi (z_n^{(j)}).$$

$$z_{n+1}^{(j)}={\hat{z}}_{n+1}^{(j)}+K_{n+1}(y_{n+1}^{(j)}-H{\hat{z}}_{n+1}^{(j)}).$$

$$K_{n+1}=C_{n+1}{H}^{\ast }{R}_{n+1},R_{n+1}=(H{C}_{n+1}{H}^{\ast }+\Gamma {)}^{-1},$$

$$u_{\mathrm{EnKF}}=\frac{1}{J}\sum _{j=1}^J{H}^{\perp }{z}_N^{(j)}=\frac{1}{J}\sum _{j=1}^J{u}_N^{(j)},$$

$${\hat{\mu }}_n(du)=\frac{1}{J}\sum _{j=1}^J{\delta }_{u_n^{(j)}}(du),$$

$${\mu }_n(du)=\frac{1}{Z_n}{e}^{-n\Delta \Phi (u;y)}{\mu }_0(du),$$

$$\Phi (u;y)=\frac{1}{2}{\Vert {\Gamma }^{-1/2}(y-\mathcal{G}(u))\Vert }_2^2.$$

$$\mu (du)=\frac{1}{Z}{e}^{-\Phi (u;y)}{\mu }_0(du),$$

$$\text{find }u\in X\text{ given observations }y={\mathcal{G}}^{\epsilon }(u)+\eta \in Y,$$

$$\{\begin{array}{rlrl}-\nabla \cdot (A_u^{\epsilon }\nabla p^{\epsilon })& =f,& & \text{ in }\Omega ,\\ p^{\epsilon }& =0,& & \text{ on }\partial \Omega ,\end{array}$$

$$A_u^{\epsilon }(x)=A\left(u(x),\frac{x}{\epsilon }\right),$$

$$\{\begin{array}{rlrl}-\nabla \cdot (A_u^0\nabla p^0)& =f,& & \text{ in }\Omega ,\\ p^0& =0,& & \text{ on }\partial \Omega .\end{array}$$

$${\Vert A_{u_1}-A_{u_2}\Vert }_{L^{\infty }(\Omega ;{\mathbb{R}}^{d\times d})}\le M{\Vert u_1-u_2\Vert }_X,A_u\xi \cdot \xi \ge {\alpha }_0{\Vert \xi \Vert }_2^2,$$

$${\Vert \mathcal{O}(p_1)-\mathcal{O}(p_2)\Vert }_Y\le C_{\mathcal{O}}{\Vert p_1-p_2\Vert }_{L^2(\Omega )},$$

$$\Vert u\Vert :=\frac{1}{J}\sum _{j=1}^J{\Vert u^{(j)}\Vert }_2,$$

$$\mathbb{E}\left[\Vert u_N^{\epsilon }-u_{N,h}^0\Vert \right]\to 0\text{ as }\epsilon ,h\to 0.$$

$$\mathbb{E}\left[\Vert u_N^{\epsilon }-u_{N,h}^0\Vert \right]\le C(\epsilon +h^{s+1}),$$

$${\mu }^{\epsilon }=\frac{1}{J}\sum _{j=1}^J{\delta }_{u_N^{\epsilon ,(j)}}\text{ and }{\mu }_h^0=\frac{1}{J}\sum _{j=1}^J{\delta }_{u_{N,h}^{0,(j)}},$$

$${\mathbb{E}}_{\xi }\left[\left|{\int }_{E}fd{\mu }_n-{\int }_{E}fd\mu \right|\right]\to 0.$$

$$\{{\mu }^{\epsilon }-{\mu }_h^0{\}}_{\epsilon ,h}\stackrel{L^1}{\rightharpoonup }0\text{ as }\epsilon ,h\to 0.$$

$$e(\epsilon ,u)=\frac{1}{J}\sum _{j=1}^J{\Vert {\mathcal{G}}^{\epsilon }(u^{(j)})-{\mathcal{G}}^0(u^{(j)})\Vert }_2,$$

$$\tilde{e}(h,u)=\frac{1}{J}\sum _{j=1}^J{\Vert {\mathcal{G}}_h^0(u^{(j)})-{\mathcal{G}}^0(u^{(j)})\Vert }_2.$$

$$\{\begin{array}{rlrl}-\nabla \cdot (A_u\nabla p)& =f,& & \text{ in }\Omega ,\\ p& =0,& & \text{ on }\partial \Omega ,\end{array}$$

Full text

AMS subject classifications. 62G05, 65N21, 74Q05.

Key words. Inverse problems, Multiscale modelling, Homogenization, Ensemble Kalman filter, Bayesian inference, Modelling error.

1 Introduction

In this work we consider the application of techniques derived from the Kalman filter to inverse problems involving multiscale phenomena which can be modelled by means of partial differential equations (PDEs). Inverse problems arise in many fields, such as seismography, meteorology and tomography, all physical domains with a multiscale nature. Our reference mathematical model is given by multiscale elliptic PDEs of the form

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ is the physical domain, $A_{u}^{\varepsilon}$ is a tensor oscillating with an amplitude described by the parameter $\varepsilon$ and $u$ is a possibly infinite-dimensional unknown which parametrizes the tensor $A^{\varepsilon}_{u}$. We are then interested in the solution of inverse problems involving the retrieval of the parameter $u$ given noisy observations derived from the solution $p^{\varepsilon}$.

Multiscale inverse problems of this form have been recently introduced in [14] and analysed extensively in [2, 3]. In particular, in [2] Abdulle and Di Blasio build a coarse-graining approach to solve the inverse problem regularized with a Tikhonov technique. The main idea is replacing the computationally expensive solution of the highly-oscillating multiscale problem with an homogenized surrogate, which eliminates the fast variables and is therefore cheaper. In particular, the theory of homogenization guarantees under certain assumptions, which will be specified throughout this work, that there exists a PDE of the form

[TABLE]

such that the solution $p^{0}$ is the weak limit of the functions $p^{\varepsilon}$ in the vanishing limit for $\varepsilon$, and such that $A_{u}^{0}$ is independent of $\varepsilon$. In [2], the authors showed that employing this homogenized model to the multiscale inverse problem guarantees a good approximation to its solution if a Tikhonov regularization is employed. This framework has been successively enlarged by the same authors to the Bayesian case in [3], where the analysis involves posterior distributions arising from both the multiscale and the homogenized model. In the same work, a technique for estimating the modelling error which was developed in [6, 5] is successfully applied to multiscale inverse problems to account for the homogenization and discretization errors.

The ensemble Kalman filter (EnKF), first introduced in [10], is an algorithm which is widely employed in the engineering community for the estimation of the state of partially-observed dynamical systems whose dynamics are governed by a nonlinear agent. In particular, Kalman filters have long been used successfully in meteorology, oceanography and automation applications. In [11], Iglesias et al. propose the application of the EnKF method to obtain a point-wise solution to inverse problems involving PDEs, and an extension of their analysis giving a Bayesian interpretation of the filtering solution is presented in [17].

In this work, we present a combination of the well-established techniques of homogenization and filtering to build a novel scheme for solving multiscale inverse problems in an efficient and reliable manner. In the same spirit of [2, 3], we prove that it is possible to eliminate the fast scales from the PDE appearing in the inverse problem relying on the theory of homogenization, thus obtaining a solution which is accurate in the vanishing limit for the multiscale parameter $\varepsilon$. In our analysis, we both consider point-wise estimations as in [11] and Bayesian solutions as in [17], thus showing convergence results which are endowed with decay rates under special assumptions on the problem. Inspired by [6, 5, 3], we then consider offline and online techniques for estimating the modelling error and prove a novel result indicating the computational cost which is required for such an estimation for any given multiscale problem.

In general, the EnKF has two main advantages with respect to other approaches. First, a Bayesian interpretation of the solution to the inverse problem is obtained from the algorithm without any additional cost. The Bayesian paradigm, frequently adopted in the context of inverse problems involving PDEs, provides a full uncertainty quantification on the solution and is therefore preferable to a point-wise estimation. Secondly, the EnKF is easily parallelizable, thus allowing in practice to solve complex inverse problems faster than employing, e.g., Markov chain Monte Carlo methods.

The main contributions of this paper are:

• •

to introduce a new method based on filtering techniques and numerical homogenization, which is computationally efficient and easy parallelizable to solve multiscale inverse problems;

• •

to analyze theoretically the convergence properties of our method both from a point-wise and a Bayesian perspectives, proving the results of convergence of the EnKF scheme in the multiscale setting;

• •

to estimate the modelling error caused by homogenization and discretization, and prove a novel theoretical results which strengthens its value in practice.

The outline of the work is the following. In Section 2 we briefly summarize the technique of ensemble Kalman inversion, show how it can be applied to multiscale inverse problems and state our main theoretical results. In Section 3 we present the analysis of our theoretical results, and Section 4 is dedicated to the estimation of the modelling error. Finally, in Section 5 we present a series of numerical experiments which corroborate our analysis.

2 Ensemble Kalman inversion for multiscale problems

In this section, we present the ensemble Kalman inversion technique for multiscale inverse problems. First, we introduce a generic framework and illustrate how the EnKF is employed to solve an inverse problem. Then we particularize to a inverse problems involving multiscale elliptic PDEs, and we conclude this section by announcing our main theoretical results. For a more exhaustive treatment of the EnKF in a generic PDE context, we refer the reader to [11, 17].

2.1 Ensemble Kalman inversion

We first give a brief summary of the ensemble Kalman inversion for problems of the form

[TABLE]

where $X$ and $Y$ are Hilbert spaces, the operator $\mathcal{G}\colon X\to Y$ is a generic forward map and the noise $\eta$ follows the Gaussian distribution $\eta\sim\mathcal{N}(0,\Gamma)$ with a symmetric positive definite covariance $\Gamma$. Kalman filters are traditionally employed to estimate the state of a dynamical system given partial and noisy observations of its state. In order to approximate the solution of the otherwise static problem (3), it is therefore natural to introduce some artificial dynamics. Let us consider the space $Z=X\times Y$ and the map $\Xi\colon Z\to Z$ given by

[TABLE]

Given an initial value $z_{0}\in Z$, we define artificial discrete dynamics on $Z$ through the recursion

[TABLE]

The dynamics on $Z$ are completed consistently with the problem (3) by the observation equation

[TABLE]

where $H\colon Z\to Y$ is the projection operator defined by $H=\begin{bmatrix}0&I\end{bmatrix}$ and $\{\eta_{n}\}_{n\in\mathbb{N}}$ is an i.i.d. sequence of random variables distributed identically to the noise of the inverse problem (3), i.e., $\eta_{n}\sim\mathcal{N}(0,\Gamma)$. In fact, let us remark that combining (5) and (6) one gets $y_{n+1}=\mathcal{G}(u_{n})+\eta_{n+1}$, which is in law equivalent to the equality appearing in (3).

Kalman filters proceed recursively to estimate the state of dynamics of the form (5) when observations are provided by the model (6). At each time $n$, the estimation is performed in two steps. First, equation (5) is employed in the so-called prediction step, and then (6) is employed to correct the prediction in the update or analysis step. In case $\Xi$ is a linear map, both prediction and update steps admit a closed-form expression, often referred to in literature as the Kalman formulae. Conversely, in case $\Xi$ is nonlinear, there exist no explicit solution to the estimation problem and one has to recur to an approximation such as the EnKF method, which we briefly describe here.

Given a positive integer $J$, the EnKF method proceeds by propagating and updating an ensemble $\{z_{n}^{(j)}\}_{j=1}^{J}\subset Z$ of particles with discrete approximations of the Kalman formulae. Let $\mathcal{A}\subset X$ be such that $\dim(\mathcal{A})\leq J$, and let the initial ensemble $\{z_{0}^{(j)}\}_{j=1}^{J}$ to be given by

[TABLE]

where $\{\psi^{(j)}\}_{j=1}^{J}\subset\mathcal{A}$. At each time $n=0,1,\ldots,N-1$, and for each $j=1,\ldots,J$, the prediction step is simply given by

[TABLE]

In the analysis step, this partially-updated ensemble is updated given knowledge of the data $y$. For better exploring the space $Y$, the data is randomized and each particle $z_{n+1}^{(j)}$ is compared to i.i.d. versions of the data given by $y_{n+1}^{(j)}=y+\eta_{n+1}^{(j)}$, where $\eta_{n+1}^{(j)}\sim\mathcal{N}(0,\Gamma)$. The analysis step is then given by

[TABLE]

The operator $K_{n+1}\colon Y\to Z$, the Kalman gain, weighs the effects of dynamics and observations in this two-step procedure, and is defined as

[TABLE]

where $C_{n+1}\colon Z\to Z$ is the empirical covariance of the partially-updated ensemble $\{\hat{z}_{n+1}^{(j)}\}_{j=1}^{J}$, the operator $H^{*}\colon Y\to Z$ is the adjoint of $H$, which is given in (6), and we recall $\Gamma$ to be the covariance of the noise $y$, so that $R_{n+1}\colon Y\to Y$. Intuitively, one can notice that when the ensemble’s covariance $C_{n+1}$ is large with respect to the noise covariance $\Gamma$, i.e., the observation model is more precise than the dynamics, we will have $z_{n+1}^{(j)}\approx y_{n+1}^{(j)}$, while in the opposite case we will have $z_{n+1}^{(j)}\approx\hat{z}_{n+1}^{(j)}$. A more precise definition of the operators appearing above will be given in Section 3. At the final step $N$, we project the particles on the space $X$ and average the result to obtain the estimate

[TABLE]

where $H^{\perp}\colon Z\to X$ is defined by $H^{\perp}=\begin{bmatrix}I&0\end{bmatrix}$. The last detail missing to fully define the EnKF is its initialization, i.e., the choice of the space $\mathcal{A}$ defining the initial ensemble. We assume prior knowledge is available on the parameter $u\in X$ and that it is summarized by a probability measure $\mu_{0}$ on $X$. In this case, one can draw $J$ i.i.d. samples $\psi^{(j)}$ from $\mu_{0}$ and fix $\mathcal{A}=\mathrm{span}\{\psi^{(j)}\}_{j=1}^{J}$.

Remark 1*.*

The computational cost of the EnKF method is approximately equal to the number of evaluations of the forward operator, which in a PDE framework dominates with respect to the algebraic operations needed in the analysis step. Therefore, the complexity of the algorithm is $\mathcal{O}(JN)$. Nonetheless, let us remark that the prediction step (7) can be easily parallelized, since the forward operator is applied independently to each particle. Hence, for a reasonable number of particles (or a high number of computing units), we have that the overall cost is of order $\mathcal{O}(N)$.

As shown in [17], a slight modification of the EnKF algorithm allows to obtain with no additional cost a Bayesian solution to (3) from the evolving ensemble. Let $\mu_{0}$ be, as above, a prior probability measure on $X$ and let the initial ensemble $\{\psi^{(j)}\}_{j=1}^{J}$ consist of i.i.d. samples from $\mu_{0}$. Given a number of steps $N$, let $\Delta=1/N$ be a “stepsize”. Let us modify the algorithm above by taking instead of the covariance $\Gamma$ of the noise its scaled version $\Delta^{-1}\Gamma$ in formula (9). Moreover, let us define the empirical measure $\hat{\mu}_{n}$ on $X$ induced by the ensemble at the $n$-th step, i.e.

[TABLE]

where $\delta_{x}$ is the Dirac mass concentrated in $x\in U$. Then, it has been shown in [17] that $\hat{\mu}_{n}$ is a good approximation of the measure $\mu_{n}$ defined by

[TABLE]

where $Z_{n}$ is the normalization constant and $\Phi(u;y)$ is the least squares functional

[TABLE]

For $n=N$, we have by definition $N\Delta=1$ and the measure $\mu\coloneqq\mu_{N}$ given by

[TABLE]

where $Z$ is the normalization constant, is exactly the posterior measure of the parameter $u$ given the prior $\mu_{0}$ in the Bayesian sense (see, e.g., [18]). Summarizing, if one carefully modifies formula (9) for the Kalman gain, it is sufficient to run the EnKF method for $N$ steps and the empirical measure given by the particles is an approximation to the Bayesian posterior.

2.2 Multiscale ensemble Kalman inversion

In this work, we consider the application of ensemble Kalman inversion to a multiscale inverse problem of the form

[TABLE]

where $\varepsilon>0$ is the multiscale parameter, which often is $\varepsilon\ll 1$, the operator $\mathcal{G}^{\varepsilon}\colon X\to Y$ is the multiscale forward map and where, as above, $\eta\sim\mathcal{N}(0,\Gamma)$ for some symmetric positive definite covariance $\Gamma$ on $Y$. Let $\Omega\subset\mathbb{R}^{d}$ be an open bounded domain and let $H_{0}^{1}(\Omega)$ denote the space of functions $v\colon\Omega\to\mathbb{R}$ in $L^{2}(\Omega)$ with first order weak derivatives in $L^{2}(\Omega)$ and whose trace on $\partial\Omega$ vanishes. We consider the forward map $\mathcal{G}^{\varepsilon}$ to be the composition $\mathcal{G}^{\varepsilon}=\mathcal{O}\circ\mathcal{S}^{\varepsilon}$ of an observation operator $\mathcal{O}\colon H_{0}^{1}(\Omega)\to Y$ and a multiscale solution operator $\mathcal{S}^{\varepsilon}\colon X\to H_{0}^{1}(\Omega)$. In particular, for $u\in X$, the operator $\mathcal{S}^{\varepsilon}\colon u\mapsto p^{\varepsilon}\in H_{0}^{1}(\Omega)$ where $p^{\varepsilon}$ is the weak solution of the elliptic PDE

[TABLE]

for a right-hand side $f\in L^{2}(\Omega)$. We assume that the tensor $A^{\varepsilon}_{u}\colon\Omega\to\mathbb{R}^{d\times d}$ is a parametrized multiscale tensor admitting explicit scale separation between slow and fast spatial variables, i.e.,

[TABLE]

where the map $(t,x)\mapsto A(t,x/\varepsilon)$ is assumed to be known and where $A$ is periodic in its second argument. In other words, the unknown $u$ of the inverse problem (13) governs the slow-scale variations of the rapidly-oscillating tensor $A_{u}^{\varepsilon}$.

Let us consider now the application of ensemble Kalman inversion to the inverse problem (13). Since the PDE (14) does not in general admit a closed-form solution, one has to employ a numerical approximation to evaluate the forward map $\mathcal{G}^{\varepsilon}$. If $\varepsilon$ is small and we employ the finite element method (FEM), a fine discretization is needed to resolve the smallest scale and thus evaluate the forward operator $\mathcal{G}^{\varepsilon}$, which clearly leads to a high computational cost. Indeed, as for Remark 1, a run of the EnKF algorithm would lead to $\mathcal{O}(N)$ solutions of (14), which is indeed unfeasible.

In order to approach the multiscale problem more efficiently we recur to the theory of homogenization (see e.g. [8]), which ensures the existence of a non-oscillating homogenized tensor $A^{0}_{u}$, such that for $\varepsilon\to 0$ the solution $p^{\varepsilon}$ of (14) tends weakly in $H_{0}^{1}(\Omega)$ to the solution $p^{0}$ of the problem

[TABLE]

Hence, this homogenized problem is a good surrogate of (14) when $\varepsilon\ll 1$, and its non-oscillating nature allows us to discretize it with FEM on an arbitrarily coarse mesh, whose maximum diameter is denoted by $h$. Therefore, denoting by $\mathcal{G}^{0}_{h}\colon\mathcal{O}\circ\mathcal{S}^{0}_{h}$, where $\mathcal{S}^{0}_{h}\colon u\mapsto p_{h}^{0}$, the numerical solution of (15), we study in this paper the behavior of the EnKF when $\mathcal{G}^{\varepsilon}$ is replaced by its cheap approximation $\mathcal{G}^{0}_{h}$. Let us denote by $\{u_{n,h}^{0,(j)}\}_{j=1}^{J}$ the ensemble obtained after $n$ iterations of the EnKF algorithm with the forward operators $\mathcal{G}_{h}^{0}$ in the prediction step (7). With this notation, given an initial ensemble $\{u_{0,h}^{0,(j)}\}_{j=1}^{J}$, at each step $n=0,1,\ldots,N-1$, our algorithm proceeds as

1. (i)

for each $u_{n,h}^{0,(j)}$, compute the homogenized tensor $A_{u_{n}^{(j)}}^{0}$ and build the forward map $\mathcal{G}_{h}^{0}$, 2. (ii)

perform the prediction step (7) with $\mathcal{G}_{h}^{0}$ and the analyis step (8) to obtain the updated ensemble $\{u_{n+1,h}^{0,(j)}\}_{j=1}^{J}$.

The computation of the homogenized tensor relies as well on numerical procedures, here we use the finite element heterogeneous multiscale method (FE-HMM) [1, 4]. Let us finally remark that similar analyses have been carried on in [14, 2, 3] for different methodologies in the solution of (13).

2.3 Statement of main results

Let us first introduce some assumptions and notation which will be employed in the analysis. First, we introduce a regularity assumption on tensors which will be fulfilled by $A^{\varepsilon}_{u}$ and $A^{0}_{u}$.

Assumption 1*.*

The tensor $A_{u}\colon\Omega\to\mathbb{R}^{d\times d}$ satisfies for all $u,u_{1},u_{2}\in X$ and $\xi\in\mathbb{R}^{d}$

[TABLE]

where $M$ and $\alpha_{0}$ are positive constants.

We now introduce a regularity assumption on the observation operator.

Assumption 2*.*

The observation operator $\mathcal{O}\colon H^{1}_{0}(\Omega)\to Y$ satisfies for all $p_{1},p_{2}\in H^{1}_{0}(\Omega)$

[TABLE]

where $C_{\mathcal{O}}$ is a positive constant.

Note that since $\mathcal{O}$ is defined on $H^{1}_{0}(\Omega)\subset L^{2}(\Omega)$, Assumption 2 is stronger than Lipschitz continuity. Finally, we introduce an assumption on the algorithm which will be employed in the analysis.

Assumption 3*.*

All the particles in the ensemble lie at each iteration in a ball $B_{R}(u^{*})$ for some $R>0$ sufficiently big, where $u^{*}$ is the true value of the unknown.

For clarity, we present the analysis the finite-dimensional setting $X=\mathbb{R}^{M}$ and $Y=\mathbb{R}^{L}$ but claim that it can be readily generalized to the infinite-dimensional case. For an ensemble $u=\{u^{(j)}\}_{j=1}^{J}$ of particles in $\mathbb{R}^{M}$, we introduce the ensemble norm

[TABLE]

which is indeed a norm and where $\left\|\cdot\right\|_{2}$ is the Euclidean norm in $\mathbb{R}^{M}$. Moreover, given a scalar $\alpha$, we define the linear combination $w=u+\alpha v$ between two ensembles $u$ and $v$ with the same number of particles $J$ as $\{w^{(j)}=u^{(j)}+\alpha v^{(j)}\}_{j=1}^{J}$.

We can now present the first main result of this work, in which we show the convergence of the ensemble obtained by the EnKF employing $\mathcal{G}^{0}_{h}$ to the one obtained employing the exact operator $\mathcal{G}^{\varepsilon}$ linked to the PDE (14).

Theorem 1**.**

Let $u_{N,h}^{0}=\{u_{N,h}^{0,(j)}\}_{j=1}^{J}$, $u_{N}^{\varepsilon}=\{u_{N}^{\varepsilon,(j)}\}_{j=1}^{J}$ be the ensembles after $N$ iterations of the EnKF method with forward operators $\mathcal{G}^{0}_{h}$ and $\mathcal{G}^{\varepsilon}$ respectively. Then, if $A^{\varepsilon}_{u}$ and $A^{0}_{u}$ satisfy Assumption 1 and if Assumption 2 and Assumption 3 hold, we have

[TABLE]

In particular, if the exact solution $p^{0}$ of the homogenized problem (15) is in $H^{q+1}(\Omega)$ with $q\geq 1$ and we employ polynomials of degree $r$ for the finite element basis, then

[TABLE]

where $s=\min\{r,q\}$ and $C>0$ is a constant independent of $h$ and $\varepsilon$.

The proof of this result is the main focus of Section 3.1. The second main theoretical result concerns the Bayesian interpretation of the EnKF methodology for inverse problems in the multiscale setting. Let $\mu_{0}$ be a prior measure on $X$ and the ensembles $u_{N,h}^{0}=\{u_{N,h}^{0,(j)}\}_{j=1}^{J}$, $u_{N}^{\varepsilon}=\{u_{N}^{\varepsilon,(j)}\}_{j=1}^{J}$ resulting from the EnKF algorithms as in Theorem 1 both initialized with an i.i.d. sample from $\mu_{0}$. We consider the discrete probability measures

[TABLE]

i.e., the EnKF approximations of the posterior $\mu$ on $u$ defined in (12). Our goal is providing a measure on how far the two measures are from each other with respect to $\varepsilon$ and $h$. Let us remark that due to the randomization of the data at each step of the EnKF algorithm, both $\mu^{\varepsilon}$ and $\mu_{h}^{0}$ are random probability measures. We now introduce the metric we consider for comparing the two measures.

Definition 1**.**

Let $(\Omega,\mathcal{A},P)$ be a probability space. A sequence of random measures $\{\mu_{n}\}_{n\in\mathbb{N}}$ on a metric space $(E,\mathcal{B}(E))$ dependent on a random variable $\xi$ on $(\Omega,\mathcal{A},P)$ is said to weakly converge in $L^{1}(\Omega)$ to a random measure $\mu$ on the same metric space if for all bounded continuous functions $f\in C^{0}_{B}(E)$ we have

[TABLE]

In this case we write $\mu_{n}\xrightharpoonup{L^{1}}\mu$.

We can now state our second main result, whose proof is the main focus of Section 3.2.

Theorem 2**.**

Let the hypotheses of Theorem 1 be satisfied. Then the sequence of random measures $\{\mu^{\varepsilon}-\mu_{h}^{0}\}_{\varepsilon,h}$, where $\mu^{\varepsilon}$ and $\mu_{h}^{0}$ are defined in (18), satisfies

[TABLE]

Remark 2*.*

It is possible to verify that in both Theorem 1 and Theorem 2 the limits with respect to $\varepsilon$ and $h$ can be interchanged.

3 Convergence analysis

In this section we prove Theorem 1 and Theorem 2, the main results if this work. As announced above, the analysis is carried out in the finite dimensional case $X=\mathbb{R}^{M}$ and $Y=\mathbb{R}^{L}$, but it can be generalized to the infinite dimensional setting. For the purpose of the analysis, we introduce on top of the forward maps $\mathcal{G}^{\varepsilon}$ and $\mathcal{G}^{0}_{h}$, which have been introduced in Section 2.2, the operator $\mathcal{G}^{0}=\mathcal{O}\circ\mathcal{S}^{0}$, where $\mathcal{S}^{0}\colon X\to H^{1}_{0}(\Omega)$ is the exact solution operator associated with the homogenized PDE (15).

3.1 Convergence of the point estimate

We now focus on Theorem 1. It is clear from the desired bound that the effects of homogenization and discretization can be analysed separately. In particular, we first show the convergence of the ensemble generated employing the forward operator $\mathcal{G}^{\varepsilon}$ to the one generated employing the exact homogenized operator $\mathcal{G}^{0}$ for $\varepsilon\to 0$. Then, in an analogous fashion, we prove the convergence of the ensemble generated with $\mathcal{G}^{0}_{h}$ to the ensemble generated employing $\mathcal{G}^{0}$. In order to introduce a compact notation, we denote by $\mathcal{U}_{J,M}$ the set of ensembles of dimension $J$ with elements in $\mathbb{R}^{M}$ and we consider the homogenization error function $e\colon\mathbb{R}\times\mathcal{U}_{J,M}\to\mathbb{R}$, which is defined for a generic ensemble $u$ as

[TABLE]

and a discretization error function $\tilde{e}\colon\mathbb{R}\times\mathcal{U}_{J,M}\to\mathbb{R}$ as

[TABLE]

Before proving the main theorem, we introduce some preliminary results.

Let us first consider a generic forward operator involving an elliptic PDE and show that the associated forward map is Lipschitz continuous.

Lemma 1**.**

Let $\mathcal{G}\colon\mathbb{R}^{M}\to\mathbb{R}^{L}$, $\mathcal{G}=\mathcal{O}\circ\mathcal{S}$ be a forward operator such that $\mathcal{O}\colon H^{1}_{0}(\Omega)\to\mathbb{R}^{L}$ is Lipschitz and $\mathcal{S}\colon\mathbb{R}^{M}\to H^{1}_{0}(\Omega)$, $\mathcal{S}\colon u\mapsto p$ is defined by the solution of

[TABLE]

where $\Omega\subset\mathbb{R}^{d}$ is an open bounded set, the right-hand side $f\in L^{2}(\Omega)$ and the tensor $A_{u}$ satisfies Assumption 1. Then $\mathcal{G}$ is Lipschitz with a constant depending only on the Poincaré constant of $\Omega$, on the constants $M$ and $\alpha$ appearing in Assumption 1, on the right-hand side $f$ and on the Lipschitz constant of the operator $\mathcal{O}$.

The proof of Lemma 1 is given in the Appendix. In the following Lemma, whose proof is also given in the Appendix, we consider the homogenization error defined in (19) and shows that it vanishes in the limit $\varepsilon\to 0$.

Lemma 2**.**

Let $e$ be defined as (19). Under Assumption 2, we have for all $u\in\mathcal{U}_{J,M}$

[TABLE]

Moreover, if the solution of the homogenized problem (15) is in $H^{2}(\Omega)$ independently of $u$, then there exists $K>0$ independent of $\varepsilon$ and $u$ such that

[TABLE]

Finally, we consider the particle empirical covariances of ensembles given by the EnKF algorithm, thus proving their boundedness and Lipschitz continuity. The proof of this Lemma can be found in the Appendix.

Lemma 3**.**

Let $C^{up}(u)\in\mathbb{R}^{M\times L}$ and $C^{pp}(u)\in\mathbb{R}^{L\times L}$ be defined as

[TABLE]

where $\bar{u}\in\mathbb{R}^{M}$ and $\bar{\mathcal{G}}\in\mathbb{R}^{L}$ are the empirical averages

[TABLE]

and let $\mathcal{G}\colon\mathbb{R}^{M}\to\mathbb{R}^{L}$ be Lipschitz with constant $C_{\mathcal{G}}$. Then, there exist four constants $C_{i}>0$, $i=1,\ldots,4$, such that

1. (i)

$\left\|C^{up}(u)\right\|_{2}\leq C_{1}$, 2. (ii)

$\left\|C^{pp}(u)\right\|_{2}\leq C_{2}$, 3. (iii)

$\left\|C^{up}(u_{1})-C^{up}(u_{2})\right\|_{2}\leq C_{3}\left\|u_{1}-u_{2}\right\|$, 4. (iv)

$\left\|C^{pp}(u_{1})-C^{pp}(u_{2})\right\|_{2}\leq C_{4}\left\|u_{1}-u_{2}\right\|$,

for all ensembles $u,u_{1},u_{2}\in\mathcal{U}_{J,M}$ which are stable in the sense of Assumption 3.

In order to clarify the exposition, we first consider the amplification the error over one step between the EnKF algorithms employing the multiscale and the homogenized forward operators respectively, which is summarized in the following lemma.

Lemma 4**.**

For all $n=0,\dots,N-1$, let $u_{n}^{0}=\{u_{n}^{0,(j)}\}_{j=1}^{J}$,$u_{n}^{\varepsilon}=\{u_{n}^{\varepsilon,(j)}\}_{j=1}^{J}$ be the ensembles of particles at the $n$-th iteration of the EnKF for the forward operators $\mathcal{G}^{0}$ and $\mathcal{G}^{\varepsilon}$ respectively. Then, under Assumption 1, Assumption 2 and Assumption 3, there exist positive constants $\alpha$ and $\gamma$ such that

[TABLE]

where $e(\varepsilon,u)$ is given in (19).

Proof.

First, due to Assumption 2 and the Poincaré inequality with constant $C_{p}$ we have

[TABLE]

which shows that $\mathcal{O}$ is Lipschitz with constant $C_{\mathcal{O}}C_{p}$. Therefore, applying Lemma 1, we deduce that both $\mathcal{G}^{0}$ and $\mathcal{G}^{\varepsilon}$ are Lipschitz with constant $C_{\mathcal{G}}$ independent of $\varepsilon$. The Kalman update formulae (8) restricted to the $u$ variable read (see [11])

[TABLE]

Combining (23) and (24), we have

[TABLE]

and using the triangle inequality we obtain

[TABLE]

where

[TABLE]

Let us introduce two useful inequalities which will be employed in the following. Given $A$ and $B$ square invertible matrices of the same size, it holds

[TABLE]

Moreover, if $A$ is positive semidefinite and $B$ is positive definite, it holds

[TABLE]

Let us first consider $S_{1}$. Applying Lemma 3 and (30) to the first two factors gives

[TABLE]

Moreover, since $y_{n+1}^{(j)}=y+\eta_{n+1}^{(j)}$ and since $y=\mathcal{G}^{\varepsilon}(u^{*})+\eta$, where $u^{*}$ is the true value of the unknown and $\eta$ is the true realization of the noise, the triangle inequality yields

[TABLE]

which, since $\mathcal{G}^{\varepsilon}$ is Lipschitz and due to Assumption 3, implies

[TABLE]

Hence, we get

[TABLE]

Finally, the random variables $\zeta_{n+1}^{(j)}\coloneqq\eta_{n+1}^{(j)}+\eta$ are i.i.d., distributed as $\zeta\sim\mathcal{N}(0,2\Gamma)$ and independent of $u_{n}^{\varepsilon}$ and $u_{n}^{0}$, which implies first

[TABLE]

and second, defining $\alpha_{1}\coloneqq C_{3}\left\|\Gamma^{-1}\right\|_{2}(C_{\mathcal{G}}R+\sqrt{2\mathrm{tr}(\Gamma)})$, yields the final bound

[TABLE]

Let us now consider the second term $S_{2}$. We apply Lemma 3 to the norm of $C^{up}(u_{n}^{0})$. Moreover, applying the inequalities (29), (30) and Lemma 3 gives

[TABLE]

Reasoning as for $S_{1}$ for the third factor appearing in (27) finally yields

[TABLE]

where $\alpha_{2}\coloneqq C_{1}C_{4}\left\|\Gamma^{-1}\right\|_{2}^{2}(C_{\mathcal{G}}R+\sqrt{2\mathrm{tr}(\Gamma)})$. We now consider the last term $S_{3}$. The first factor appearing in (28) can be bounded by Lemma 3 and for the second factor we use (30), thus obtaining

[TABLE]

Regarding the third factor of (28), we apply the triangle inequality and the Lipschitz continuity of the forward operator $\mathcal{G}^{\varepsilon}$, which yield

[TABLE]

Substituting back into $S_{3}$ and by definition of $e(\varepsilon,u_{n}^{0})$ and of the ensemble norm we obtain

[TABLE]

Therefore, defining $\alpha_{3}=C_{1}\left\|\Gamma^{-1}\right\|_{2}C_{\mathcal{G}}$ and $\gamma=C_{1}\left\|\Gamma^{-1}\right\|_{2}$ we have the bound

[TABLE]

Finally, defining $\alpha\coloneqq 1+\alpha_{1}+\alpha_{2}+\alpha_{3}$, and using the results (25), (34), (36) and (38), we obtain the desired result. ∎

We now present the main result about global multiscale convergence of the EnKF algorithm.

Proposition 1**.**

Under the notation and assumptions of Lemma 4, letting $u_{0}^{\varepsilon}=u_{0}^{0}$ be the same initial ensemble, we have

[TABLE]

Moreover, if the solution of the homogenized problem (15) is sufficiently regular, namely $p^{0}\in H^{2}(\Omega)$, then there exists $K_{1}>0$ independent of $\varepsilon$ such that

[TABLE]

Proof.

Since $u_{0}^{\varepsilon}=u_{0}^{0}$, iterating the estimate of Lemma 4 yields

[TABLE]

Applying Lemma 2, we have $e(\varepsilon,u_{i}^{0})\to 0$ for all $i=0,\dots,N-1$, hence as $\varepsilon\to 0$

[TABLE]

Moreover, if $p^{0}$ belongs to $H^{2}(\Omega)$, applying Lemma 2 gives

[TABLE]

where $K_{1}=\gamma(\alpha^{N}-1)K/(\alpha-1)$, which is the desired result. ∎

We now consider convergence with respect to the FEM discretization of the homogenized problem. First, we introduce a preliminary result, which plays the role of Lemma 2 in the context of numerical convergence and whose proof is given in the Appendix.

Lemma 5**.**

Let $\tilde{e}$ be defined in (20) and let Assumption 2 hold. If the exact solution $p^{0}$ of the homogenized problem (21) is in $H^{q+1}(\Omega)$, the right-hand side $f$ is in $H^{q-1}(\Omega)$ and we employ polynomials of degree $r$ for the finite element basis, then

[TABLE]

where $s=\min\{r,q\}$.

We can now state the main result concerning convergence with respect to the numerical discretization of the homogenized problem.

Proposition 2**.**

Let $u_{N}^{0}=\{u_{N}^{0,(j)}\}_{j=1}^{J}$, $u_{N,h}^{0}=\{u_{N,h}^{0,(j)}\}_{j=1}^{J}$ be the ensembles of particles at the last iteration of the iterative ensemble Kalman filter for the forward operators $\mathcal{G}^{0}$ and $\mathcal{G}_{h}^{0}$ respectively. Then, under Assumption 1, Assumption 2, Assumption 3 and if the exact solution $p^{0}$ of the homogenized problem (21) is in $H^{q+1}(\Omega)$ and we use polynomials of degree $r$ for the finite element basis, we have

[TABLE]

where $s=\min\{r,q\}$ and $K_{2}$ is a positive constant independent of $h$.

Proof.

The proof of Proposition 2 is identical to the proof of Proposition 1, except that all the ensembles $\{u_{n}^{\varepsilon}\}_{n=1}^{N}$ obtained by the multiscale operator $\mathcal{G}^{\varepsilon}$ have to be replaced by the ensembles $\{u_{n,h}^{0}\}_{n=1}^{N}$ obtained by the finite element discretization of the homogenized operator $\mathcal{G}^{0}_{h}$. Moreover Lemma 2 for the error $e$ has to be replaced by Lemma 5 for the error $\tilde{e}$. ∎

Applying Proposition 1 and Proposition 2, we finally prove Theorem 1.

Proof of Theorem 1.

An application of the triangle inequality yields

[TABLE]

The two addends can be bounded applying Proposition 1 and Proposition 2, thus obtaining the desired result for $C=\max\{K_{1},K_{2}\}$. ∎

3.2 Convergence of the posterior distributions

In this section, we give the proof of Theorem 2, i.e., the convergence of the discrete posterior measures $\mu^{\varepsilon}$ to $\mu^{0}_{h}$ introduced in (18) as $\varepsilon,h\to 0$. Let $u^{*}\in\mathbb{R}^{M}$ and let $B_{R}(u^{*})$ be the ball of radius $R$ centered in $u^{*}$ with respect to the norm $\left\|\cdot\right\|_{s}$ with $s\in[1,\infty]$. Due to the discrete nature of these distributions, we study convergence with respect to the Wasserstein metrics, for which we report its standard definition in the metric spaces $(B_{R}(u^{*}),\left\|\cdot\right\|_{s})$, which can be found, e.g., in [16].

Definition 2**.**

Let $\mu$ and $\nu$ be two probability measures on the metric space $(B_{R}(u^{*}),\left\|\cdot\right\|_{s})$. The Wasserstein distance between $\mu$ and $\nu$ is defined for all $p\in[1,\infty)$ as

[TABLE]

where $\Gamma(\mu,\nu)$ denotes the collection of all joint distributions on $B_{R}(u^{*})\times B_{R}(u^{*})$ with marginals $\mu$ and $\nu$ on the first and second factors respectively.

Remark 3*.*

If $\mu$ and $\nu$ are two discrete distributions on finite state spaces, respectively $\Omega_{1}=\{u_{1},\dots,u_{K_{1}}\}$ and $\Omega_{2}=\{v_{1},\dots,v_{K_{2}}\}$ included in $B_{R}(u^{*})$, then (40) can be written as

[TABLE]

where the matrix $\gamma$ has to satisfy the following constraints

[TABLE]

We now show that the distance $W_{1,2}$ is bounded by the distance induced by the ensemble norm defined in (17). This result will be crucial later to prove Theorem 1.

Lemma 6**.**

Let $u_{1}=\{u_{1}^{(j)}\}_{j=1}^{J}$, $u_{2}=\{u_{2}^{(j)}\}_{j=1}^{J}$ be two ensembles of particles and let $\mu_{1},\mu_{2}$ be the corresponding distributions defined as sum of Dirac masses

[TABLE]

Then for all $s\in[1,\infty]$ and $p\in[1,\infty)$ it holds

[TABLE]

and, in particular,

[TABLE]

Proof.

Take $\gamma^{*}$ defined as

[TABLE]

which satisfies the constraints (42), and note that

[TABLE]

Therefore, by definition of Wasserstein distance for discrete distributions on finite spaces (41), we deduce that

[TABLE]

which is the desired result. Finally, taking $p=1$ and $s=2$ and recalling the ensemble norm defined in (17), we obtain the second inequality. ∎

We now analyze the relationship between the weak $L^{1}$ convergence introduce in Definition 1 and the convergence with respect to the expectation of the Wasserstein distance for random probability measures. In particular, we prove that the latter implies the former, which was already proved in [16] for non-random measures. Here, we extend the result to random probability measures. The proof of the following Lemma is given in the Appendix.

Lemma 7**.**

Let $(\Omega,\mathcal{A},P)$ be a probability space. Let the sequence $\{\mu_{n}\}_{n\in\mathbb{N}}$ and $\mu$ be random probability measures on the metric space $(B_{R}(u^{*}),\left\|\cdot\right\|_{s})$ dependent on the random variable $\xi$ on $(\Omega,\mathcal{A},P)$. If

[TABLE]

then $\mu_{n}\xrightharpoonup{L^{1}}\mu$.

We can now complete the proof of Theorem 2.

Proof of Theorem 2.

Applying Lemma 6 and due to Theorem 1, we deduce that for $\varepsilon,h\to 0$ it holds

[TABLE]

Note that the only difference in the update step of the EnKF when used for a point estimate and in the Bayesian framework is that $\Gamma$ is replaced by $\Delta^{-1}\Gamma$ where $\Delta=1/N$. The constants of the proof of Theorem 1 depend on $\left\|\Gamma^{-1}\right\|_{2}$, which is now replaced by $\left\|(\Delta^{-1}\Gamma)^{-1}\right\|_{2}$, which can be bounded by $\left\|\Gamma^{-1}\right\|_{2}$ as

[TABLE]

Finally, applying Lemma 7, we obtain the desired result. ∎

4 Modelling error

In this section, we consider the effects of model misspecification due to the homogenization and discretization error. All the results presented in Section 3 deal with the asymptotic case $h,\varepsilon\to 0$, which is unrealistic in applications. Let us recall that the original inverse problem involves predicting the exact unknown $u^{*}$ from observations originated by the model

[TABLE]

where $\eta\sim\mathcal{N}(0,\Gamma)$ is the noise. Since evaluating $\mathcal{G}^{\varepsilon}$ is too expensive and in many applications unfeasible, we wish to employ the cheaper forward operator $\mathcal{G}^{0}_{h}$. Hence, we rewrite (43) as

[TABLE]

where

[TABLE]

The quantity $\mathcal{E}(u^{*})$ represents the error introduced by misspecification of the forward model. Equation (44) shows that the observed data $y$ can be seen as data originating by the discrete homogenized model which is affected by two sources of errors, the original noise and the modelling error. This formulation of modelling error was originally presented in [6], and then applied to multiscale inverse problems in [3]. Following [6, 3], we assume that the modelling error is a Gaussian random variable independent of the noise $\eta$, so that $\mathcal{E}\sim\mathcal{N}(m,\Sigma)$ for all $u$, and write

[TABLE]

where $\zeta\sim\mathcal{N}(0,\Sigma)$. There is no theoretical guarantee for the modelling error to be distributed as a Gaussian in this framework. Nevertheless, it has been shown in [13] that in the one-dimensional case a Gaussian assumption can be employed effectively for the modelling error, thus partially justifying our choice. Then we define

[TABLE]

and, from (45), we obtain

[TABLE]

Therefore, if the mean $m$ and covariance $\Sigma$ of the modelling error are known, a more reliable approximation of the unknown $u^{*}$ can be obtained applying the EnKF to (46). The modelling error distribution, by assumption fully determined by its mean and covariance, is approximated offline. We sample $N_{\mathcal{E}}$ unknowns $\{u_{i}\}_{i=1}^{N_{\mathcal{E}}}$ from $\mu_{0}$ and, for all $i=1,\dots,N_{\mathcal{E}}$, we apply both the forward operators $\mathcal{G}^{\varepsilon}(u_{i})$ and $\mathcal{G}^{0}_{h}(u_{i})$. Then we compute

[TABLE]

and the mean $m$ and the covariance $\Sigma$ are obtained as the empirical mean and covariance of the sample $\{\mathcal{E}_{i}\}_{i=1}^{N_{\mathcal{E}}}$. This procedure is computationally involved due to the multiple evaluations of $\mathcal{G}^{\varepsilon}$, but it has to be performed only once and can then be applied to different sets of observations and true values $u^{*}$. Let us also remark that on the one hand, due to the theory of homogenization, the modelling error can be considered negligible when $\varepsilon$ is very small, and the expensive estimation of $\mathcal{E}$ may not be necessary. On the other hand, when $\varepsilon$ is larger, the homogenized equation does not provide with a good approximation of the multiscale problem, and an estimation of $\mathcal{E}$ is required. One may rightfully argue that in case $\varepsilon=\mathcal{O}(1)$, it is possible to evaluate the forward operator $\mathcal{G}^{\varepsilon}$ without a large computational effort. Hence, the techniques presented in this section are relevant for mid-range values of $\varepsilon$, for which $\mathcal{E}$ is significant with respect to the noise $\eta$. Moreover, we remarked in practice via numerical experiments that a small number $N_{\mathcal{E}}$ can be employed to obtain a satisfactory approximation of the modelling error. A theoretical justification of this property is provided by Theorem 3 and Theorem 4.

In order to obtain a more reliable approximation of the distribution of the modelling error, we can follow a dynamic approach based on the estimation of the mean $m$ and the covariance $\Sigma$ online, i.e., during the run of the EnKF algorithm. This methodology has been developed in [5]. In particular, we sequentially apply the ensemble Kalman method for $\mathcal{L}$ levels and, at each level $\ell=1,\dots,\mathcal{L}$, we update the distribution of the modelling error, which is denoted by $\nu^{\ell}=\mathcal{N}(m^{\ell},\Sigma^{\ell})$. Letting

[TABLE]

be the approximation of the distribution of the particles at iteration $n$ at level $\ell$, $\mu_{0}^{\ell+1}=\mu_{N^{\ell}}^{\ell}$ and $\mu_{0}^{1}=\mu_{0}$, where $N^{\ell}$ is the number of iterations at level $\ell$, then the mean $m^{\ell}$ and the covariance $\Sigma^{\ell}$ are approximated as in the offline approach with the only difference that $\mu_{0}$ is replaced by $\mu_{0}^{\ell}$. This approach provides indeed a better approximation of the modelling error as instead of taking the samples from the prior distribution, they are drawn from distributions which are progressively closer to the true posterior. On the other hand, this procedure has to be done online and it is computationally expensive because it requires the resolution of $N_{\mathcal{E}}=\sum_{\ell=1}^{\mathcal{L}}N_{\mathcal{E}}^{\ell}$ full multiscale problems.

Finally, we are interested in studying whether the simple offline method for estimating the modelling error provides indeed a good approximation. In this direction, we give in Theorem 3 and Theorem 4 a criterion on how to choose the number $N_{\mathcal{E}}$ of full multiscale problems which has to be solved in order to have a reliable approximation of the true mean $m^{*}$ and covariance $\Sigma^{*}$ of the modelling error with respect to $\varepsilon$ and $h$. Before stating Theorem 3 and Theorem 4, let us recall the Hoeffding’s and McDiarmid’s inequalites, which will be used in the proofs. Let $\{Y_{i}\}_{i=1}^{N}$ be independent random variables with values in $[a,b]$, and let $\bar{Y}$ be the sample average of $\{Y_{i}\}_{i=1}^{N}$. Then, the Hoeffding’s inequality states that for all $\eta\in\mathbb{R}$ it holds

[TABLE]

Moreover, let $\{X_{i}\}_{i=1}^{N}$ be independent random variables with values in the space $\mathcal{X}$, and let $\varphi\colon\mathcal{X}^{N}\to\mathbb{R}$ satisfy for all $i=1,\dots,N$

[TABLE]

then the McDiarmid’s inequality states that for all $\eta\in\mathbb{R}$ it holds

[TABLE]

Theorem 3**.**

Let $\alpha\in(0,1)$, $\eta>0$ and $C_{\mathcal{E}}=\max\{K,\tilde{K}\}$, where $K$ and $\tilde{K}$ are the constants of Lemma 2 and Lemma 5. Let $\{\mathcal{E}_{i}\}_{i=1}^{N_{\mathcal{E}}}\subset\mathbb{R}^{L}$ be given by

[TABLE]

for a sample of realizations $\{u_{i}\}_{i=1}^{N_{\mathcal{E}}}$ from the standard normal distribution $\mathcal{N}(0,I)$, let $m$ be the sample mean of $\{\mathcal{E}_{i}\}_{i=1}^{N_{\mathcal{E}}}$ and $m^{*}=\mathbb{E}[\mathcal{E}_{i}]$. If

[TABLE]

where $s$ is given by Lemma 5, then

[TABLE]

Proof.

First, note that the modelling error is bounded, indeed by Lemma 2 and Lemma 5, we have for each $i=1,\dots,N_{\mathcal{E}}$

[TABLE]

so each component $(\mathcal{E}_{i})_{l}$, for $l=1,\dots,L$, is bounded by the same constant

[TABLE]

Observe that if

[TABLE]

then

[TABLE]

which implies that

[TABLE]

Using (49) and applying Hoeffding’s inequality we have

[TABLE]

Define the events $A_{l}=\left\{\left\lvert m_{l}-m^{*}_{l}\right\rvert\leq\frac{\eta}{\sqrt{L}}\right\}$ for each $l=1,\dots,L$, then we have

[TABLE]

and, applying the De Morgan’s laws and the union bound, we obtain

[TABLE]

Therefore, thanks to (52), (53) and (54), we have

[TABLE]

and if $N_{\mathcal{E}}$ satisfies the hypothesis we obtain the desired result. ∎

Theorem 4**.**

Let $\alpha\in(0,1)$, $\eta>0$ and $C_{\mathcal{E}}=\max\{K,\tilde{K}\}$, where $K$ and $\tilde{K}$ are the constants of Lemma 2 and Lemma 5. Let $\{\mathcal{E}_{i}\}_{i=1}^{N_{\mathcal{E}}}\subset\mathbb{R}^{L}$ be given by

[TABLE]

for a sample of realizations $\{u_{i}\}_{i=1}^{N_{\mathcal{E}}}$ from the standard normal distribution $\mathcal{N}(0,I)$, let $m$ and $\Sigma$ be the sample mean and covariance of $\{\mathcal{E}_{i}\}_{i=1}^{N_{\mathcal{E}}}$ and $m^{*}=\mathbb{E}[\mathcal{E}_{i}]$ and $\Sigma^{*}=\mathbb{E}[(\mathcal{E}_{i}-m)(\mathcal{E}_{i}-m)^{T}]$. If

[TABLE]

where $s$ is given by Lemma 5 and $\widehat{C}$ is specified in the proof, then

[TABLE]

Proof.

First, repeating verbatim the first part of the proof of Theorem 3 we have

[TABLE]

Observe that if

[TABLE]

then denoting by $\left\|\cdot\right\|_{F}$ the Frobenius norm we have

[TABLE]

which implies that

[TABLE]

For all $j,k=1,\dots,L$ define the functions $\varphi_{j,k}\colon(\mathbb{R}^{L})^{N}\to\mathbb{R}$ as

[TABLE]

where

[TABLE]

and the function $\Phi\colon(\mathbb{R}^{L})^{N}\to\mathbb{R}^{L\times L}$ whose component $(j,k)$ is given by $\varphi_{j,k}$. Observe that $\Sigma=\Phi(\mathcal{E}_{1},\dots,\mathcal{E}_{N})$ and $\Sigma^{*}=\mathbb{E}[\Phi(\mathcal{E}_{1},\dots,\mathcal{E}_{N})]$. Since the modelling error is bounded, we can restrict the functions $\varphi_{j,k}$ to the ball of radius $C_{\mathcal{E}}(\varepsilon+h^{s+1})$ centred in [math], $\varphi_{j,k}\colon(B_{C_{\mathcal{E}}(\varepsilon+h^{s+1})})^{N}\to\mathbb{R}$, allowing us to prove the following bound

[TABLE]

where

[TABLE]

In fact we have

[TABLE]

where

[TABLE]

Now we bound the two terms separately. First, we have

[TABLE]

Let $x,y\in B_{C_{\mathcal{E}}(\varepsilon+h^{s+1})}$ and note that for all $j=1,\dots,L$ we have

[TABLE]

By equations (62) and (64), it holds for all $j=1,\dots,L$

[TABLE]

therefore we obtain

[TABLE]

Moreover, by (64) we also get

[TABLE]

witch, together with (61), (65) and the fact that $(2N-1)/(N-1)\leq 3$ for all $N\geq 2$, implies (60)

[TABLE]

Therefore, applying McDiarmid’s inequality we have

[TABLE]

where $\widehat{C}=2304$. Finally, we define the events $A_{j,k}=\left\{\left\lvert\Sigma_{j,k}-\Sigma^{*}_{j,k}\right\rvert\leq\frac{t}{L}\right\}$ for each $j,k=1,\dots,L$ and we repeat the same argument as in the last part of the proof of Theorem 3. Hence, due to (59) and (68) we have

[TABLE]

and if $N_{\mathcal{E}}$ satisfies the hypothesis we obtain the desired result. ∎

Remark 4*.*

Note that, in Theorem 3 and Theorem 4, as expected, the number $N_{\mathcal{E}}$ of full multiscale problems tends to infinity if we require no error between the sample and the true mean and covariance ($\eta\to 0$) or certainty that the error is below a certain value ($\alpha\to 0$). Moreover, observe that for any given accuracy the number of samples required $N_{\mathcal{E}}$ is a increasing function of $\varepsilon$ and $h$, so that if the model $\mathcal{G}^{0}_{h}$ is a good approximation of $\mathcal{G}$, thus computationally expensive, then only few samples are needed. In particular, notice that in order to obtain a good approximation of the true mean, the number of full multiscale problems is

[TABLE]

while to have a reliable approximation of the covariance matrix it is required that

[TABLE]

5 Numerical experiments

In this section, using the setting of [3], we present some numerical experiments to illustrate the iterative ensemble Kalman method to solve multiscale inverse problems.

Let $\Omega$ be a bounded open domain. We consider a class of parametrized multiscale locally periodic tensors of the type $A^{\varepsilon}_{\sigma^{*}}(x)=A(\sigma^{*}(x),x/\varepsilon)$, where $\sigma^{*}\colon\Omega\to\mathbb{R}$. We assume to know the map $(t,x)\to A(t,x/\varepsilon)$ for all $x\in\Omega$ and $t\in\mathbb{R}$ and we want to estimate the function $\sigma^{*}$ given measurements computed from the model

[TABLE]

Remark 5*.*

Note that the theory has been developed for Dirichlet homogeneous boundary conditions, but it can be applied to the non-homogeneous case by considering an extension of the function at the boundary and slightly modifying the PDE. For more details we refer to [15, Remark 8.10].

For the unknown $\sigma^{*}$ we consider the following admissible set

[TABLE]

where $\sigma^{-}$ and $\sigma^{+}$ are two given values.

The measurements, which we take into account, are the integrals of the normal flux multiplied by some functions with compact support in a portion of the boundary of the domain. More precisely, we consider $I\in\mathbb{N}$ disjoint portions of $\Omega$, which we denote by $\Gamma_{i}\in\partial\Omega$, $i=1,\dots,I$, $\Gamma_{i}\cap\Gamma_{j}=\emptyset$ for $i\neq j$, and $I$ functions $\varphi_{i}\in H^{1/2}(\partial\Omega)$ with compact support $\mathrm{supp}\;(\varphi_{i})\subset\Gamma_{i}$ for all $i=1,\dots,I$. Moreover, we solve (72) for $K\in\mathbb{N}$ Dirichlet data $g_{k}$, $k=1,\dots,K$, and we denote by $p_{k}^{\varepsilon}$ the solution of the problem. Let $\Lambda_{A^{\varepsilon}_{\sigma}}\colon H^{1/2}(\partial\Omega)\to H^{-1/2}(\partial\Omega)$ be the operator which maps the Dirichlet data $g$ to the normal flux of the solution $p^{\varepsilon}$ of (72)

[TABLE]

where $\nu$ is the exterior unit normal vector to $\partial\Omega$. Then we define the multiscale operator $\mathcal{F}^{\varepsilon}\colon\Sigma\to\mathbb{R}^{L}$ where $L=IK$ by components

[TABLE]

which, with an abuse of notation, can be written

[TABLE]

The final vector of observations $y$ is given by the sum of the operator $\mathcal{F}^{\varepsilon}$ and a noise

[TABLE]

where $\eta\sim\mathcal{N}(0,\Gamma)$ and $\Gamma$ is a given symmetric positive definite covariance matrix, which, in our experiments, is a multiple of the identity $\Gamma=\gamma^{2}I$ and $\gamma$ is a given value. Observations are computed with a refined Finite Element Method (FEM) with mesh size $h_{\mathrm{obs}}\ll\varepsilon$, while the homogenized version of problem (72) is solved using a macro mesh size $h\gg h_{\mathrm{obs}}$. We call $\mathcal{T}_{h}$ the macro triangulation and $N_{h}$ the total number of nodes defining $\mathcal{T}_{h}$. We assume that the prior distribution for the discretization of the unknown $\sigma^{*}$ on the macro triangulation $\mathcal{T}_{h}$ is given by $\mathcal{N}(\sigma_{0},C)$, where $\sigma_{0}$ is a given discretization of a function in $\Sigma$ and $C\in\mathbb{R}^{N_{h}\times N_{h}}$ is defined by

[TABLE]

where $\delta,\lambda\in\mathbb{R}^{+}$ and $\{x_{i}\}_{i=1}^{N_{h}}$ are the nodes of the macro triangulation $\mathcal{T}_{h}$. The parameter $\lambda$ is a correlation length that describes how the values at different positions of the functions supported by the prior measure are related, while the parameter $\delta$ is an amplitude scaling factor. Regarding the prior modelling, we need to take into account that even if in the homogenized problem the coarse and fine scales have been separated, functions drawn from the prior distribution on the coarse scale can exhibit multiple scales, including the fine scale of our multiscale model, depending on the rate of decay of the prior covariance. This issue can thus be controlled by setting the parameters $\delta$ and $\lambda$. Even though this does not ensure a clear separation between coarse and fine scales, our numerical results illustrate that it is sufficient in practice.

In order to reduce the dimensionality of the unknown we use a truncated Karhunen-Loève expansion. Any sample from the prior distribution $\mathcal{N}(\sigma_{0},C)$ can be represented as

[TABLE]

where $\{\psi_{m}\}_{m=1}^{N_{h}}$ is an orthonormal set of eigenvectors of $C$ with corresponding eigenvalues $\{\lambda_{m}\}_{m=1}^{N_{h}}$ in decreasing order, and $\{u_{m}\}_{m=1}^{N_{h}}$ is an i.i.d sequence with $u_{m}\sim\mathcal{N}(0,1)$. Note that the Karhunen-Loève expansion works also in the infinite dimensional setting, where $\sigma_{0}\in\Sigma$, $C$ is a covariance operator and $\{\lambda_{m},\psi_{m}\}_{m=1}^{\infty}$ is an orthonormal set of eigenvalues-eigenfunctions with respect to the scalar product in $L^{2}(\Omega)$. Then the truncated Karhunen-Loève expansion of the discretization of $\sigma$ consists of taking the first $M$ components of the series in (76)

[TABLE]

and the actual unknown becomes the vector $u\in\mathbb{R}^{M}$, whose components are the coefficients $u_{m}$ in (77). Then we define the multiscale forward operator $\mathcal{G}^{\varepsilon}\colon\mathbb{R}^{M}\to\mathbb{R}^{L}$ as the composition of $\mathcal{F}^{\varepsilon}$ with the truncated Karhunen-Loève expansion

[TABLE]

In the iterative ensemble Kalman method we do not compute the exact solution of problem (72), but we solve its homogenized version numerically using the macro triangulation $\mathcal{T}_{h}$, therefore we obtain the homogenized discrete solution $p^{0}_{h}$. The problem is solved applying the finite element heterogeneous multiscale method (FE-HMM), which is described in [1, 4]. Hence, analogously to the multiscale case, we define the discrete homogenized operator $\mathcal{F}_{h}^{0}\colon\Sigma\to\mathbb{R}^{L}$ with an abuse of notation as

[TABLE]

and the discrete homogenized forward operator $\mathcal{G}^{0}_{h}\colon\mathbb{R}^{M}\to\mathbb{R}^{L}$, which is actually used in the algorithm, as

[TABLE]

Finally, we call $u_{\mathrm{EnKF}}$ the solution of the iterative ensemble Kalman algorithm and the estimated $\sigma_{\mathrm{EnKF}}$ is obtained from the truncated Karhunen-Loève expansion

[TABLE]

5.1 Data

In the numerical results presented in the following section the computational domain is the unit square

[TABLE]

For the discretization parameters we set $\varepsilon=1/64$ and $h_{\mathrm{obs}}=1/4096$ and for the forward homogenized problem we use a macro mesh size $h=1/32$, which is much larger than $h_{\mathrm{obs}}$ and reduces the computational cost significantly. We solve the problem for $K=3$ Dirichlet conditions $\{g_{k}\}_{k=1}^{3}$ and $g_{k}=\sqrt{\mu_{k}}\vartheta_{k}$ where $\{(\mu_{k},\vartheta_{k})\}_{k=1}^{3}$ are couples of eigenvalues and eigenfunctions of the one dimensional discrete Laplacian operator corresponding to the first $K=3$ smallest eigenvalues. For each $g_{k}$ we consider its restriction to the boundary $\partial\Omega$ in order to obtain a Dirichlet condition. These functions are orthonormal with respect to the scalar product in $L^{2}(\Omega)$ and this ensures that each function gives independent information.

To compute the boundary integrals in (74) and (78), we consider $I=12$ boundary portions, three for each side of the square $\Omega$. In particular, for each side, all $\Gamma_{i}$ have length equal to $0.2$ and they consist of the intervals $(0.1,0.3),(0.4,0.6)$ and $(0.7,0.9)$. The functions $\{\varphi_{i}\}_{i=1}^{12}$ are hat functions with $\mathrm{supp}\;(\varphi_{i})=\Gamma_{i}$, which take value one at the midpoint and value [math] at the extremes of $\Gamma_{i}$. Then the parameter of the noise, which perturbs the observations, is $\gamma=0.01$.

Moreover, regarding the prior distribution for the unknown, we consider $\sigma_{0}=0$ and the parameters of the covariance matrices are $\delta=0.05$ and $\lambda=0.5$. In the truncated Karhunen-Loève expansion we take $M=100$. Finally, about the ensemble Kalman method, we consider $J=1000$ particles for each ensemble and $500$ iterations.

The exact tensor $A_{\sigma^{*}}^{\varepsilon}$ is given by

[TABLE]

where

[TABLE]

and

[TABLE]

Figure 1 shows the exact unknown $\sigma^{*}$. Note that $\sigma^{*}$ is a non-continuous function, but, in order to approximate it, we are using a truncated Karhunen-Loève expansion, where the eigenfunctions are smooth.

One can verify that the tensor $A^{\varepsilon}_{\sigma}$ satisfies Assumption 1. In particular, for $\xi\in\mathbb{R}^{2}$ we have

[TABLE]

Moreover, since the $\mathrm{EnKF}$ algorithm estimates the coefficients $\{u_{m}\}_{m=1}^{M}$ of the truncated Karhunen-Loève expansion, we show that $A^{\varepsilon}(u)\colon\mathbb{R}^{M}\to L^{\infty}(\Omega,\mathbb{R}^{d\times d})$, which maps $u$ into $A^{\varepsilon}_{\sigma_{u}}$, is Lipschitz. In fact we first have

[TABLE]

then using the truncated Karhunen-Loève expansion and the Cauchy-Schwarz inequality we obtain

[TABLE]

which shows that $A^{\varepsilon}(u)$ is Lipschitz with constant equal to $\sqrt{13}e^{\sigma^{+}}\sup_{x\in\Omega}\left(\sum_{m=1}^{M}\lambda_{m}\psi_{m}^{2}(x)\right)^{1/2}$.

5.2 Results

We first fix the multiscale parameter $\varepsilon=1/32$ and the ensemble size $J=500$ and study the evolution with respect to the number of steps. In Figure 2 we plot the estimation $\sigma_{\mathrm{EnKF}}$ after $10,50,250$ and $500$ iterations of the ensemble Kalman algorithm. We clearly see that the approximation gets better as the number of iterations increases and that convergence has been reached. In particular, already after $N=250$ iterations the algorithm seem to have reached convergence. We point out that we obtain a quite good approximation of the real unknown $\sigma^{*}$ , indeed we are trying to recover a non-continuous function in the whole domain given only some observations at the boundary.

We now perform a sensitivity analysis with respect to the ensemble size. In Figure 3 we vary the number of particles $J$ and we compare the results obtained at the end of the algorithm after $500$ iterations for $\varepsilon=1/32$. As expected, the approximation becomes better when the ensemble contains more particles. In particular, note that if the number of particles is too small, e.g. $J=10$, then the approximation is not satisfying.

Further, we fix the ensemble size $J=500$ and we perform $N=500$ iterations of the EnKF for different values of the multiscale parameter. Results, shown in Figure 4, highlight how the approximation becomes worse when $\varepsilon$ is bigger, indeed the homogenized problem becomes too different with respect to the multiscale one and, if $\varepsilon$ is too big, the solution does not approximate the true unknown.

Moreover, in order to obtain good results even in case $\varepsilon$ is not close to the asymptotic limit $\varepsilon\to 0$, in Figure 5 we apply offline modelling error estimation with $N_{\mathcal{E}}=20$ and we plot the solution of the inverse problem (46) for different values of the multiscale parameter $\varepsilon$. Comparing these plots with the ones in Figure 4, in particular for $\varepsilon=1/4$, we observe that the modelling error estimation significantly improves the results.

Finally, in Figure 6 we show the results obtained by applying the ensemble Kalman method with dynamic updating of the modelling error distribution with $\mathcal{L}=5$ levels, $N_{\mathcal{E}}^{\ell}=4$ samples and $N^{\ell}=100$ iterations at each level $\ell=1,\dots,\mathcal{L}$. The number of resolutions of the full multiscale problem is $20$ and the total number of iterations is $500$, which are equal to the previous approach, where the distribution of the modelling error was approximated offline. Comparing these plots with the ones in Figure 5, we note that updating the distribution of the modelling error dynamically still improves the results.

6 Conclusion

In this paper we analyzed the ensemble Kalman inversion methodology in the context of inverse problems for multiscale elliptic PDEs with tensors highly oscillatory at a scale $\varepsilon\ll 1$. The multiscale algorithm we propose relies on the EnKF, on a surrogate homogenized forward operator and on numerical homogenization techniques such as the FE-HMM. It guarantees a significant reduction in computational cost for problems which would be otherwise computationally involved or unfeasible. In Theorem 1 we have shown that the ensemble of particles approximating the unknown parameter generated by our multiscale algorithm converges to the ensemble generated by the true model as the small scale parameter $\varepsilon$ and the numerical discretization parameter $h$ go to zero. Furthermore in a Bayesian framework, we have shown in Theorem 2 that the discrete probability measure based on the ensemble originating from our multiscale algorithm converges to the measure generated by the true model, again as $\varepsilon$ and $h$ go to zero. Hence when $\varepsilon\ll 1$ and the full model is expensive to solve, the multiscale numerical method we propose is both accurate and efficient to recover an unknown parameter in multiscale elliptic PDEs. Moreover, we equipped our method with a technique which allows to account for the discrepancy between the artificial homogenized surrogate forward model and the true multiscale data, thus alleviating the effects of model misspecification. This technique requires additional offline or online computations involving the numerical solution of the full multiscale problem. The optimal number of such additional solves is quantified in Theorem 3 and Theorem 4. In particular, we have proved that the number of solves needed to reach any required accuracy tends to zero when the small scale parameter $\varepsilon$ and the numerical discretization parameter $h$ vanish. Hence, we can conclude that accounting for model misspecification is particularly beneficial for mid-range values of $\varepsilon$, when a small number of full solves should be computationally affordable. The efficiency and usefulness of the multiscale algorithm have been further demonstrated through a series of numerical experiments.

Appendix

Proof of Lemma 1

Let $u_{1},u_{2}\in\mathbb{R}^{M}$, and $p_{1}=\mathcal{S}(u_{1})$, $p_{2}=\mathcal{S}(u_{2})$. From the weak formulations of (21) we get that

[TABLE]

which yields

[TABLE]

Then choosing $v=p_{1}-p_{2}$, by the hypotheses on $A_{u}$ and applying the Hölder inequality we obtain

[TABLE]

which due a standard coercivity argument implies

[TABLE]

where $C_{p}$ is the Poincaré constant associated to the domain $\Omega$. Hence (79) shows that $\mathcal{S}$ is Lipschitz with constant

[TABLE]

Finally, since $\mathcal{G}$ is the composition of two Lipschitz operators, we deduce that it is also Lipschitz with constant $L_{\mathcal{G}}=L_{\mathcal{O}}L_{\mathcal{S}}$. ∎

Proof of Lemma 2

Let us consider an ensemble $u\in\mathcal{U}_{J,M}$ with particles $u^{(j)}\in\mathbb{R}^{M}$, for $j=1,\ldots,J$. For each particle we have

[TABLE]

where we write explicitly the dependence of the solutions $p^{\varepsilon}$ and $p^{0}$ on the particle they are generated by. Due to homogenization theory, we have that $p^{\varepsilon}(u^{(j)})\rightharpoonup p^{0}(u^{(j)})$ in $H^{1}_{0}(\Omega)$ for all $j=1,\dots,J$, and therefore $p^{\varepsilon}(u^{(j)})\to p^{0}(u^{(j)})$ in $L^{2}(\Omega)$, which implies

[TABLE]

Moreover, if the solution of the homogenized problem $p^{0}$ is sufficiently smooth independently of $u$, namely $p^{0}\in H^{2}(\Omega)$, letting $C>0$ be a constant independent of $\varepsilon$, we have by [12] for all $j=1,\dots,J$

[TABLE]

which implies

[TABLE]

and defining $K=C_{\mathcal{O}}C$ gives the desired result. ∎

Proof of Lemma 3

First, for all $x\in B_{R}(u^{*})$ we have

[TABLE]

We can also deduce the same bounds for the mean values

[TABLE]

Then by (81) and (82) we get

[TABLE]

and defining $C_{1}=4Mm$ we get $(i)$. The argument is similar for the matrix $C^{pp}(u)$, for which we have

[TABLE]

and defining $C_{2}=4M^{2}$ we get $(ii)$. Before proving $(iii)$ and $(iv)$, we need the following estimates for two ensemble of particles $u_{1}$ and $u_{2}$

[TABLE]

Then we have

[TABLE]

and since $\mathcal{G}$ is Lipschitz and due to (81), (82), (83), we obtain

[TABLE]

and defining $C_{3}=2(J+1)(mC_{\mathcal{G}}+M)$ we get $(iii)$. The argument is similar for the matrix $C^{pp}(u)$, for which we have

[TABLE]

and defining $C_{4}=4(J+1)MC_{\mathcal{G}}$ we get $(iv)$, which concludes the proof. ∎

Proof of Lemma 5

Let us consider an ensemble $u\in\mathcal{U}_{J,M}$ with particles $u^{(j)}\in\mathbb{R}^{M}$, for $j=1,\ldots,J$. For each particle we have

[TABLE]

where we write explicitly the dependence of the solutions $p^{0}$ and $p^{0}_{h}$ on the particle they are generated by. Then due to standard a priori error estimates of FEM (see e.g. [7, Theorem 3.2.5]) and higher order boundary regularity results for elliptic partial differential equations (see e.g. [9, Theorem 6.3.5]) we have for all $j=1,\ldots,J$

[TABLE]

where $C>0$ is a constant independent of $h$. Therefore, we obtain

[TABLE]

and defining $\tilde{K}=C_{\mathcal{O}}C\left\|f\right\|_{H^{q-1}(\Omega)}$ gives the desired result. ∎

Proof of Lemma 7

We follow the same steps of the proof of Theorem 5.9 in [16]. Let us first recall the duality formula for the Wasserstein distance with $p=1$

[TABLE]

where $\Phi$ is the set of all globally Lipschitz continuous functions $\varphi\colon B_{R}(u^{*})\to\mathbb{R}$ with Lipschitz constant $C_{\mathrm{Lip}}\leq 1$. Note that if $\varphi\in\Phi$, then also $-\varphi\in\Phi$. Hence we deduce that

[TABLE]

Then we have

[TABLE]

where the right hand side vanishes by hypothesis. Therefore we obtain

[TABLE]

for all $\varphi\in\Phi$. Finally, we extend (90) to all Lipschitz functions by linearity and to all bounded continuous functions by density, thus proving the desired result. ∎