# Reduced order modeling for elliptic problems with high contrast diffusion coefficients

Albert Cohen, Matthieu Dolbeault, Agustin Somacal, Wolfgang Dahmen

► To cite this version:

Albert Cohen, Matthieu Dolbeault, Agustin Somacal, Wolfgang Dahmen. Reduced order modeling for elliptic problems with high contrast diffusion coefficients. ESAIM: Mathematical Modelling and Numerical Analysis, In press. hal-03549810v2

HAL Id: hal-03549810

https://hal.science/hal-03549810v2

Submitted on 20 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Reduced order modeling for elliptic problems with high contrast diffusion coefficients

Albert Cohen*, Wolfgang Dahmen†, Matthieu Dolbeault, Agustin Somacal*

April 20, 2023

## Abstract

We consider a parametric elliptic PDE with a scalar piecewise constant diffusion coefficient taking arbitrary positive values on fixed subdomains. This problem is not uniformly elliptic, as the contrast can be arbitrarily high, contrarily to the Uniform Ellipticity Assumption (UEA) that is commonly made on parametric elliptic PDEs. We construct reduced model spaces that approximate uniformly well all solutions with estimates in relative error that are independent of the contrast level. These estimates are sub-exponential in the reduced model dimension, yet exhibiting the curse of dimensionality as the number of subdomains grows. Similar estimates are obtained for the Galerkin projection, as well as for the state estimation and parameter estimation inverse problems. A key ingredient in our construction and analysis is the study of the convergence towards limit solutions of stiff problems when diffusion tends to infinity in certain domains.

# 1 Introduction

## 1.1 Reduced models for parametrized PDEs

Parametric PDE's are commonly used to describe complex physical phenomena. With $y = (y_1, \dots, y_d)$ denoting a parameter vector ranging in some domain $Y \subset \mathbb{R}^d$, and $u(y)$ the corresponding solution to the PDE of interest, assumed to be well defined in some Hilbert space $V$, we denote by

$$ \mathcal{M} := \{u(y) : y \in Y\}, \qquad (1.1) $$

the collection of all solutions, called the *solution manifold*.

There are two main ranges of problems associated to parametric PDEs:

<footnote>*Laboratoire Jacques-Louis Lions, Sorbonne Université, 4 place Jussieu, 75005 Paris, France (albert.cohen@sorbonne-universite.fr, matthieu.dolbeault@sorbonne-universite.fr and agustin.somacal@sorbonne-universite.fr)</footnote>

<footnote>†Department of Mathematics, University of South Carolina, 1523 Greene St, Columbia, SC 29208, USA (wolfgang.anton.dahmen@googlemail.com)</footnote>
1. Forward modeling: in applications where many queries of the parameter to solution map $y \mapsto u(y)$ are required, one needs numerical forward solvers that efficiently compute approximations $\tilde{u}(y)$ with a prescribed accuracy.

2. Inverse problems: when the exact value of the parameter $y$ is unknown, one is interested in either recovering an approximation to $u(y)$ (state estimation) or to $y$ (parameter estimation), from a limited number of observations $z_i = \ell_i(u(y))$, possibly corrupted by noise.

*Reduced order modeling* is widely used for tackling both problems. In its most common form, its aim is to construct linear spaces $V_n$ of moderate dimension $n$ that approximate all solutions $u(y)$ with best possible certified accuracy. The natural benchmark for measuring the performance of such linear reduced models is provided by the *Kolmogorov n-width* of the solution manifold

$$d_n(\mathcal{M})_V := \inf_{\dim(V_n)=n} \mathrm{dist}(\mathcal{M}, V_n)_V \quad (1.2)$$

that describes the performance of an optimal space. Here

$$\mathrm{dist}(\mathcal{M}, V_n)_V := \sup_{u \in \mathcal{M}} \inf_{v \in V_n} \|u - v\|_V = \sup_{u \in \mathcal{M}} \|u - P_{V_n} u\|_V,$$

where $P_{V_n}$ is the $V$-orthogonal projector onto $V_n$. We refer the reader to [27] for a general treatment of $n$-widths.

While an optimal space achieving the above infimum is usually out of reach, there exist two main approaches aiming to construct “sub-optimal yet good” spaces. The first one consists in building expansions of the parameter to solution map, for example by polynomials

$$u_n(y) := \sum_{\nu \in \Lambda_n} u_{\nu} y^{\nu}, \quad y^{\nu} := y_1^{\nu_1} \dots y_d^{\nu_d}, \quad (1.3)$$

where $\Lambda_n \subset \mathbb{N}^d$ is a set of cardinality $n$. The coefficients $u_{\nu}$ are elements of $V$ and therefore, for all $y \in Y$ the approximation $u_n(y)$ is picked from the space

$$V_n := \mathrm{span}\{u_{\nu} : \nu \in \Lambda_n\}.$$

Notice that $u_n(y)$ is not the orthogonal projection $P_{V_n} u(y)$ in this case, but $u_n(y)$ is easy to compute for a given query $y$ once the $u_{\nu}$ have been constructed (usually through a high fidelity finite element solver). We refer to [5, 6, 7, 11, 17, 18, 32] for instances of this approach.

The second approach is the reduced basis method [21, 28, 29], that consists in taking

$$V_n := \mathrm{span}\{u^1, \dots, u^n\},$$

where the $u^j = u(y^j)$ are particular solution instances corresponding to a selection of parameter vectors $y^j \in Y$. A close variant is the proper orthogonal decomposition method [16, 34, 35], where the reduced spaces are obtained by principal component analysis applied to large training set of such instances. In
the reduced basis method, the parameter vectors $y^1, \dots, y^n$ can be selected by a greedy algorithm, introduced in [33] and originally studied in [15]. For such a selection process, it is proved in [12, 19] that if $d_n(M)_V$ has a certain algebraic or exponential rate of decay with $n$, then a similar rate is achieved by $\text{dist}(M, V_n)_V$ for the reduced basis spaces.

It follows that the reduced basis spaces constructed by the greedy algorithm are close to optimal. This is in contrast to the spaces $V_n$ spanned by the polynomial coefficients $u_\nu$ for which the approximation rate is not guaranteed to be optimal. We refer to [10] for instances where reduced basis methods can be proved to converge with a strictly higher rate than polynomial approximations. On the other hand, the polynomial constructions (1.3) have certain numerical advantages. Namely, for several relevant classes of parametrized PDEs, it can be shown that the parameter to solution mapping $y \mapsto u(y)$ has certain smoothness properties that can be used to obtain a-priori bounds on the $\|u_\nu\|_V$ without actually computing these norms. This allows an a priori selection of an appropriate set $\Lambda_n$ and the proof of concrete approximation estimates for the error $\sup_{y \in Y} \|u(y) - u_n(y)\|_V$. These estimates in turn provide an upper bound for $d_n(M)_V$, and therefore for reduced basis approximations.

## 1.2 Parametrized elliptic PDEs

One prototypal instance where the convergence analysis described above has been deeply studied is the parametrized second order elliptic equation

$$ - \operatorname{div}(a(y)\nabla u(y)) = f \quad \text{in } \Omega, \quad u|_{\partial\Omega} = 0 \quad \text{on } \partial\Omega, \qquad (1.4) $$

where $\Omega \subset \mathbb{R}^m$ is the spatial domain, $f \in H^{-1}(\Omega)$ is a source term, and $a(y)$ has the affine form

$$ a(y) = \bar{a} + \sum_{j=1}^{d} y_j \psi_j, \qquad (1.5) $$

with $\bar{a}$ and $(\psi_1, \dots, \psi_d)$ some fixed functions in $L^\infty(\Omega)$.

The corresponding solution $u(y) \in H_0^1(\Omega)$ is defined through the standard variational formulation in $H_0^1(\Omega)$ equipped with its usual norm. Up to renormalization, it is usually assumed that the $y_j$ range in $[-1, 1]$, or equivalently $Y = [-1, 1]^d$. To ensure existence and uniqueness of solutions, one typically assumes that the so-called *Uniform Ellipticity Assumption* (UEA) holds: for some fixed $0 < r \le R < \infty$,

$$ r \le a(x, y) \le R, \quad x \in \Omega, \quad y \in Y, \qquad (1.6) $$

where $a(x, y) := a(y)(x) = \bar{a}(x) + \sum_{j=1}^d y_j \psi_j(x)$, or in short $r \le a(y) \le R$ for all $y \in Y$. Under this assumption, Lax-Milgram theory ensures that the solution map $y \mapsto u(y)$ is well defined from $Y$ into $H_0^1(\Omega)$, with the uniform bound

$$ \|u(y)\|_{H_0^1} := \|\nabla u(y)\|_{L^2} \le \frac{C_f}{r}, \quad y \in Y. $$
Here and throughout this paper

$$C_f := \|f\|_{H^{-1}}. \qquad (1.7)$$

It was proved in [6, 32] that, under UEA, polynomial approximations (1.3) of given total degree converge sub-exponentially: for $\Lambda_n = \{|\nu| \le k\}$ with $n = \binom{k+d}{d}$, one has

$$\sup_{y \in Y} \|u(y) - u_n(y)\|_{H_0^1} \le C' \exp(-cn^{1/d}), \qquad (1.8)$$

Such sub-exponential rates show that the spaces $V_n$ based on polynomial expansions or reduced bases perform significantly better than standard finite element spaces, at least for a moderate number $d$ of parameters. It is possible to maintain a rate of convergence as $d$ grows, and even when $d = \infty$, when assuming some anisotropy in the variable $y_j$ through the decay of the size of $\psi_j$ as $j \to \infty$, see in particular [11, 17, 18] for results of this type.

## 1.3 High contrast problems

The Uniform Ellipticity Assumption (1.6) implies that there is a uniform control on the level of contrast in the diffusion function

$$\kappa(y) := \frac{\max_{x \in \Omega} a(x, y)}{\min_{x \in \Omega} a(x, y)} \le \frac{R}{r}, \quad y \in Y. \qquad (1.9)$$

This assumption also plays a key role in the derivation of the above approximation results, since it guarantees that the parameter to solution map has a holomorphic extension to a sufficiently large complex neighbourhood of $Y$. In this case, a good polynomial approximation $u_n$ may be defined by simply truncating the power series $\sum_{\nu \in \mathbb{N}^d} u_\nu y^\nu$, leading to the estimate (1.8).

On the other hand, there exist various situations where one would like to avoid such a strong restriction on the level of contrast. Perhaps the most representative setting is when the domain $\Omega$ is partitioned into disjoint subdomains $\{\Omega_1, \dots, \Omega_d\}$, each of them admitting a constant diffusivity level that could vary strongly between subdomains. This is typically the case when modeling diffusion in materials having multiple layers or inclusions that could have very different nature, for example air or liquid versus solid. This situation can be encountered in groundwater flow applications, where certain subdomains correspond to cavities, for which the diffusion function becomes nearly infinite, as opposed to subdomains containing sediments or other porous rocks.

In such a case, we do not want to limit the contrast level. To represent this setting, we let

$$a(y)|_{\Omega_j} = y_j, \quad y_j \in ]0, \infty[ \qquad (1.10)$$

or equivalently $a(y) = \sum_{j=1}^d y_j \chi_{\Omega_j}$, which corresponds to the affine form (1.5) with $\bar{a} = 0$ and $\psi_j = \chi_{\Omega_j}$, now with

$$Y := ]0, \infty[. \qquad (1.11)$$
We take (1.11) as the definition of the parameter domain Y for the remainder of this paper. The solution $u(y)$ satisfies the variational formulation

$$ \sum_{j=1}^{d} y_j \int_{\Omega_j} \nabla u(y) \cdot \nabla v \, dx = \langle f, v \rangle_{H^{-1}, H_0^1}, \quad v \in H_0^1(\Omega), \quad (1.12) $$

or equivalently $-y_j \Delta u(y) = f$ as elements of $H^{-1}(\Omega_j)$ on each $\Omega_j$, with the standard jump conditions $[a(y) \partial_n u(y)] = 0$ across the boundaries between sub-domains.

Let us observe that in this setting, it is hopeless to find spaces $V_n$ that approximate all solutions $u(y)$ uniformly well. Indeed, the following homogeneity property obviously holds: for any $y \in Y$ and $t > 0$, one has

$$ u(ty) = t^{-1}u(y). \quad (1.13) $$

This property implies in particular that $\|u(y)\|_{H_0^1}$ tends to infinity as $y \to 0$, and so does $\|u(y) - P_{V_n} u(y)\|_{H_0^1}$ in general. In fact, this also shows that the solution manifold $\mathcal{M}$ is not relatively compact and does not have finite n-widths.

In addition to this principal difficulty, let us remind that when using the spaces $V_n$ in forward modeling, we typically use the Galerkin method, that delivers the orthogonal projection onto $V_n$ however for the energy norm

$$ \|v\|_y^2 := \sum_{j=1}^{d} y_j \int_{\Omega_j} |\nabla v|^2 dx. \quad (1.14) $$

This approximation is thus optimal in $H_0^1(\Omega)$, however up to the constant $\kappa(y)^{1/2}$, which deteriorates with high contrast.

The main contribution of this paper is to treat these issues, and derive approximation estimates that are robust to high contrast, in the sense that they are independent of $y \in Y$.

Due to the main objection coming from the homogeneity property (1.13), it is natural to look for uniform approximation estimates in relative error, that is, estimates of the form

$$ \|u(y) - P_{V_n} u(y)\|_{H_0^1} \le \varepsilon_n \|u(y)\|_{H_0^1}, \quad y \in Y, \quad (1.15) $$

with $\lim_{n \to \infty} \varepsilon_n = 0$, and similarly for $P_{V_n}^y u(y)$. Our main results, Theorems 3.7 and 4.2, exhibit spaces $V_n$ ensuring the validity of such uniform estimates with $\varepsilon_n$ having sub-exponential decay with $n$, similar to the known results under UEA.

**Remark 1.1.** High contrast problems have been the object of intense investigation, in particular with the objective of developing techniques for multilevel or domain decomposition preconditioning [2, 3, 20] and a-posteriori error estimation [4, 9], that are provably robust with respect to the level of contrast. To our knowledge, the present work is the first in which this robustness is established for reduced modeling methods.
## 1.4 Outline

Throughout this paper, we consider the parametrized elliptic PDE (1.4) with $a(y)$ having piecewise constant form (1.10) over a fixed partition. In view of the homogeneity property (1.13), we are led to consider the subset

$$ Y' := [1, \infty]^d \qquad (1.16) $$

of parameters corresponding to the coercive regime. Any result on relative approximation error that is established for $Y'$ extends automatically to all of $Y$ because of the homogeneity property. Accordingly, we let

$$ \mathcal{B} := \{u(y) : y \in Y'\}. \qquad (1.17) $$

In §2, we start by proving that $\mathcal{B}$ is a precompact set of $H_0^1(\Omega)$. One crucial ingredient for this analysis are the *limit solutions* of the so-called *stiff problem*, obtained as $y_j \to \infty$ for certain $j \in \{1, \dots, d\}$.

In §3, we construct specific reduced model spaces for which the approximation estimate (1.15) holds with $\varepsilon_n$ decaying sub-exponentially. Our construction is based on partitioning the parametric domain $Y'$ into rectangular regions and using a different polynomial approximations on each region. This results in global reduced model space $V_n$ for which the accuracy bound remains sub-exponential, however in $\exp(-cn^{1/2d-2})$. A key ingredient for establishing these sub-exponential rates is the derivation of quantitative estimates on the convergence of $u(y)$ towards limit solutions defined in §2 as some $y_j$ tend to infinity. These estimates are established under an additional geometrical assumption on the partition, similar results for a general partition of $\Omega$ being an open problem.

In §4, we discuss the use of these reduced model spaces in forward modeling and inverse problems. Our main result relative to forward modeling is that the estimate (1.15) also holds for the Galerkin projection with the same exponential decay $\varepsilon_n$. We show that such a result is only possible if $V_n$ includes functions that have constant values over some subdomains. For the state estimation problem, we follow the Parametrized Background Data Weak (PBDW) method [13, 24], and obtain recovery bounds that are uniform over $y \in Y$ in relative error. For the parameter estimation problem, we introduce an ad-hoc strategy that specifically exploits the piecewise constant structure of the diffusion coefficient and obtain similar recovery bounds for the inverse diffusivity.

We conclude in §5 by presenting some numerical illustrations revealing the effectiveness of the reduced model spaces even in the high-contrast regime, as expressed by the approximation results.

**Acknowledgements:** We thank the anonymous reviewers for their constructive comments. We also thank François Murat for useful discussions in the understanding of the convergence process towards limit solutions, Hamza Maimoune for leading us to this work through his remarks during his master project, and Jules Pertinand for useful discussions.
## 2 Uniform approximation in relative error

In this section we work under no particular geometric assumption on the partition {$\Omega_1, \dots, \Omega_d$} of $\Omega$, and consider the solution manifold $\mathcal{M}$ defined by (1.1), where $u(y) \in H_0^1(\Omega)$ is solution to the elliptic boundary value problem with variational formulation (1.12). Our objective is to show the existence of spaces $V_n$ that uniformly approximate $\mathcal{M}$ in the relative error sense expressed by (1.15).

### 2.1 Limit solutions and the extended solution manifold

Our first observation is that this collection can be continuously extended when $y_j = \infty$ for some values of $j$, through limit solutions of stiff inclusions problems. Such limit solutions have for example been considered in the context of homogenization, see e.g. p.98 of [22].

For this purpose, to any $S \subset \{1, \dots, d\}$, we associate the space

$$V_S := \{v \in H_0^1(\Omega) : \nabla v|_{\Omega_j} = 0, j \in S\}. \qquad (2.18)$$

In other words, $V_S$ consists of the functions from $H_0^1(\Omega)$ that have constant values on the subdomains $\Omega_j$ for $j \in S$ (or on each of their connected components if these subdomains are not connected). It is a closed subspace of $H_0^1(\Omega)$. We decompose the parameter vector $y$ according to

$$y = (y_S, y_{S^c}), \quad y_S := (y_j)_{j \in S} \quad \text{and} \quad y_{S^c} := (y_j)_{j \in S^c}. \qquad (2.19)$$

For any finite and positive vector $y_{S^c}$, similar to the $\| \cdot \|_y$ norm (1.14), we may define

$$\|v\|_{y_{S^c}}^2 := \sum_{j \in S^c} y_j \int_{\Omega_j} |\nabla v|^2 dx, \qquad (2.20)$$

which is a semi-norm on $H_0^1(\Omega)$, and a full norm equivalent to the $H_0^1$-norm on $V_S$. Also note that when $y = (y_S, y_{S^c})$ is finite, one then has $\|v\|_{y_{S^c}} = \|v\|_y$ for any $v \in V_S$.

For any finite and positive vector $y_{S^c}$, we define the function $u_S(y_{S^c}) \in V_S$ solution to the following stiff inclusions problem:

$$\sum_{j \in S^c} y_j \int_{\Omega_j} \nabla u_S(y_{S^c}) \cdot \nabla v \, dx = \langle f, v \rangle_{H^{-1}, H_0^1}, \quad v \in V_S. \qquad (2.21)$$

The following result shows that this solution is well defined and is the limit of $u(y)$, when $y_{S^c}$ is fixed and $y_j \to \infty$ for $j \in S$. Note that the weak convergence is established in [22] (p. 98) and so we concentrate the proof on the strong convergence.

**Lemma 2.1.** There exists a unique $u_S(y_{S^c}) \in V_S$ solution to (2.21), which is the limit in $H_0^1(\Omega)$ of the solution $u(y_S, y_{S^c})$ as $y_j \to \infty$ for all $j \in S$.
*Proof.* Using the bilinear form $(u, v) \mapsto \sum_{j \in S^c} y_j \int_{\Omega_j} \nabla u \cdot \nabla v \, dx$ in the space $V_S$, Lax-Milgram theory implies the existence of a unique solution $u_S(y_{S^c}) \in V_S$ to (2.21).

Consider now a sequence $(y^n)_{n \ge 1} \in Y^N$, with $y_{S^c}^n = y_{S^c}^n$ and $y_j^n \to \infty$ for all $j \in S$. Denoting $u_n = u(y^n)$, it is readily seen that $(u_n)_{n \ge 1}$ is uniformly bounded in $H_0^1$ norm by $C = C_f c^{-1}$, where $c := \min_{n \ge 1} \min_{1 \le j \le d} y_j^n > 0$, and that any weak limit of a sequence extraction is solution to the variational equation (2.21). Therefore the whole sequence $(u_n)_{n \ge 1}$ weakly converges to $\bar{u} = u_S(y_{S^c})$.

We finally prove strong convergence by writing

$$ 
\begin{aligned}
c \|u_n - \bar{u}\|_{H_0^1}^2 & \le \int_{\Omega} a(y^n) |\nabla (u_n - \bar{u})|^2 dx \\
& = \langle f, u_n \rangle_{H^{-1}, H_0^1} - 2 \langle \bar{u}, u_n \rangle_{y_{S^c}} + \|\bar{u}\|_{y_{S^c}}^2 \\
& \xrightarrow{n \to \infty} \langle f, \bar{u} \rangle_{H^{-1}, H_0^1} - \|\bar{u}\|_{y_{S^c}}^2 = 0.
\end{aligned}
$$

□

The above lemma allows us to readily extend the solution manifold by introducing

$$ \tilde{Y} := ]0, \infty]^d, $$

and

$$ \overline{M} := \{u(y) : y \in \tilde{Y}\}, $$

where we have formally set

$$ u(y) := u_S(y_{S^c}), $$

when $y_j = \infty$ for $j \in S$ and $y_j < \infty$ for $j \in S^c$. Note that when $S = \{1, \dots, d\}$ the space $V_S$ is trivial and one has

$$ u(\infty, \dots, \infty) = 0. $$

**Remark 2.2.** Although we do not make explicit use of it, it can be checked that despite the fact that $y_j = 0$ is excluded in the definition of $\overline{M}$, it indeed coincides with the closure of $M$ in $H_0^1(\Omega)$ due to the fact that $\|u(y)\|_{H_0^1} \to \infty$ as $y \to 0$.

**Remark 2.3.** More precisely, when some $y_j$ tend to zero, $u(y)$ converges to the solution of the so-called soft inclusions problem (see [22], chapter 3), outside the corresponding subdomains $\Omega_j$. Here, due to the fact that the approximation estimates that we prove further are in relative error, these other limit solutions are of no use in our analysis.
## 2.2 A compactness result

As already observed in the introduction, the manifold $\overline{\mathcal{M}}$ is not bounded in $H_0^1(\Omega)$ due to the homogeneity property (1.13) and therefore not compact.

In order to treat this defect, we consider

$$ \tilde{Y}' := [1, \infty]^d, $$

and the submanifold

$$ \bar{\mathcal{B}} := \{u(y) : y \in \tilde{Y}'\}, $$

which is now bounded in $H_0^1(\Omega)$, from the standard a-priori estimate

$$ \|u(y)\|_{H_0^1} \le \frac{C_f}{\min y_j} \le C_f, $$

that is obtained by taking $v = u(y)$ in the variational formulation (1.12), with $C_f = \|f\|_{H^{-1}}$ as in (1.7). This estimate trivially extends to $u_S(y_{S^c})$ when the $y_j$ have infinite value for $j \in S$. In addition we have the following result.

**Theorem 2.4.** The set $\bar{\mathcal{B}}$ is compact in $H_0^1(\Omega)$.

*Proof.* Consider any sequence of vectors $y^n = (y_1^n, \dots, y_d^n) \in \tilde{Y}'$ for $n \ge 1$. We need to prove that the corresponding sequence of solutions $(u(y^n))_{n \ge 1}$ admits a converging subsequence. For this purpose, we observe that there exists a subset $S \in \{1, \dots, d\}$ such that, up to subsequence extraction,

$$ \lim_{n \to \infty} y_j^n = \infty, \quad j \in S, $$

and

$$ \lim_{n \to \infty} y_j^n = y_j < \infty, \quad j \in S^c. $$

Note that $S$ could be empty, for instance in the case where the $y_j^n$ are uniformly bounded for all $j$.

Let $\varepsilon > 0$. Using the strong convergence result in Lemma 2.1, for all $n \ge 1$ there exists an auxiliary vector $\tilde{y}^n$ such that $\tilde{y}_j^n = y_j^n$ when $y_j^n < \infty$, $\tilde{y}_j^n < \infty$ when $y_j^n = \infty$, such that by having picked $\tilde{y}_j^n$ large enough in the second case

$$ \|u(y^n) - u(\tilde{y}^n)\|_{H_0^1} \le \varepsilon/3. $$

In addition we may assume that $\tilde{y}_j^n \to \infty$ for $j \in S$. Next we introduce the vector $\tilde{y}^n$ such that $\tilde{y}_j^n = \tilde{y}_j^n$ when $j \in S$ and $\tilde{y}_j^n = y_j$ when $j \in S^c$. Applying again Lemma 2.1, we find that with $y_{S^c} = (y_j)_{j \in S^c}$, one has

$$ \|u(\tilde{y}^n) - u_S(y_{S^c})\|_{H_0^1} \le \varepsilon/3, $$

for $n$ sufficiently large. Finally we argue that

$$ \|u(\tilde{y}^n) - u(\bar{y}^n)\|_{H_0^1} \le \varepsilon/3, $$
for $n$ large enough. This is a consequence of the following variant of Strang first lemma (which proof is similar and left as an exercise to the reader) that says that for two diffusion functions $\bar{a}$ and $\tilde{a}$, the corresponding solution $\bar{u}$ and $\tilde{u}$ with the same data $f$ satisfy

$$ \|\bar{u} - \tilde{u}\|_{H_0^1} \le \frac{C_f \|\bar{a} - \tilde{a}\|_{L^\infty}}{\min\{\bar{a}_{\min}, \tilde{a}_{\min}\}^2}. $$

We then apply this to $\bar{a} := \bar{a}_n = a(\bar{y}^n)$ and $\tilde{a} := \tilde{a}_n = a(\tilde{y}^n)$, observing that from their definition, $\|\bar{a} - \tilde{a}\|_{L^\infty} = \max_{j \in S^c} |\bar{y}_j^n - y_j| \to 0$ as $n \to \infty$. Therefore $\|u(y^n) - u_S(y_{S^c})\|_{H_0^1} \le \varepsilon$ for $n$ sufficiently large, which concludes the proof. $\square$

We next observe that any $y \in Y$ can be rewritten as

$$ y = t\tilde{y}, $$

with $\tilde{y} \in Y'$ and normalization $\min \tilde{y}_j = 1$, for some $t > 0$, and from (1.13) one has $u(y) = t^{-1}u(\tilde{y})$. This motivates the study of the further reduced manifold

$$ \mathcal{N} := \{u(y) : y \in \tilde{Y}', \min y_j = 1\}, \qquad (2.22) $$

which is a subset of $\bar{B}$.

One important observation is that the solutions contained in $\mathcal{N}$ are also uniformly bounded from below, under mild assumptions on the data $f$.

**Lemma 2.5.** The set $\mathcal{N}$ is compact in $H_0^1(\Omega)$. Moreover, one has the framing

$$ \min_{1 \le j \le d} \|f\|_{H^{-1}(\Omega_j)} \le \|u(y)\|_{H_0^1} \le C_f, \qquad (2.23) $$

for all $u(y) \in \mathcal{N}$.

*Proof.* The compactness of $\mathcal{N}$ follows from that of $\bar{B}$, since $\mathcal{N}$ is a closed subset of $\bar{B}$. For the framing, as $a(y) \ge 1$ on $\Omega$,

$$ \|u\|_{H_0^1}^2 \le \sum_{j \in S^c} y_j \int_{\Omega_j} |\nabla u(y)|^2 dx = \langle f, u(y) \rangle_{H^{-1}, H_0^1} \le C_f \|u(y)\|_{H_0^1}, $$

so $\|u(y)\|_{H_0^1} \le C_f$. Now take $j \in \{1, \dots, d\}$ such that $y_j = 1$, and consider $\phi \in H_0^1(\Omega_j)$. Then

$$ \langle f, \phi \rangle_{H^{-1}, H_0^1} = \int_{\Omega_j} \nabla u(y) \cdot \nabla \phi \, dx \le \|u(y)\|_{H_0^1(\Omega)} \|\phi\|_{H_0^1(\Omega_j)}, $$

which gives the result. $\square$

<figure><img src="image_2.png" /></figure>

In the sequel of this paper, we always work under the condition that the lower bound in (2.23) is strictly positive

$$ c_f := \min_{1 \le j \le d} \|f\|_{H^{-1}(\Omega_j)} > 0. \qquad (2.24) $$
Let us observe that when $f$ is a function in $L^2(\Omega)$, this is ensured as soon as $f$ is not identically zero on one of the $\Omega_j$. We thus have

$$0 < c_f \le \|u(y)\|_{H_0^1} \le C_f, \qquad (2.25)$$

for all $u(y) \in \mathcal{N}$.

**Remark 2.6.** The condition $c_f > 0$ is in general necessary for controlling $\|u(y)\|_{H_0^1}$ from below. Indeed assume $\|f\|_{H^{-1}(\Omega_j)} = 0$ for some $j$ such that $\bar{\Omega} \setminus \bar{\Omega}_j$ is connected. Then taking $y_k = \infty$ for $k \ne j$ and $y_j = 1$, we find that $u(y) \in V_S$ with $S = \{j\}^c$, which is equivalent to $u(y) \in H_0^1(\Omega_j)$ since it vanishes on the other sub-domains. As $\|f\|_{H^{-1}(\Omega_j)} = 0$, we obtain $u(y) = 0$.

**Remark 2.7.** One also has the uniform framing in the $\| \cdot \|_y$ norm since

$$0 < c_f \le \|u(y)\|_{H_0^1} \le \|u(y)\|_y = \sqrt{\langle f, u \rangle_{H^{-1}, H_0^1}} \le C_f, \qquad (2.26)$$

for all $u(y) \in \mathcal{N}$ when all $y_j$ are finite.

The framing (2.25) has an implication on the existence of reduced model spaces that approximate uniformly well all solutions $u(y) \in \bar{\mathcal{M}}$ in relative error.

**Theorem 2.8.** There exists a sequence of linear spaces $(V_n)_{n \ge 1}$ such that $\dim(V_n) = n$, and a sequence $(\varepsilon_n)_{n \ge 1}$ that converges to zero such that

$$\|u(y) - P_{V_n} u(y)\|_{H_0^1} \le \varepsilon_n \|u(y)\|_{H_0^1} \qquad (2.27)$$

for all $y \in \tilde{Y}$, where $P_{V_n}$ is the $H_0^1(\Omega)$-orthogonal projector onto $V_n$.

*Proof.* Since $\mathcal{N}$ is compact, there exists a sequence of spaces $(V_n)_{n \ge 1}$ with $\dim(V_n) = n$ and a sequence $(\sigma_n)_{n \ge 1}$ that tends to 0, such that

$$\|v - P_{V_n} v\|_{H_0^1} \le \sigma_n, \quad v \in \mathcal{N}.$$

Now let $y \in \tilde{Y}$ differing from $(\infty, \dots, \infty)$, for which there is nothing to prove since $u(\infty, \dots, \infty) = 0$, and let $t^{-1} = \min_{1 \le j \le d} y_j < \infty$. By homogeneity, $t^{-1}u(y) = u(ty) \in \mathcal{N}$, and therefore

$$\|u(y) - P_{V_n} u(y)\|_{H_0^1} = t \|u(ty) - P_{V_n} u(ty)\|_{H_0^1(\Omega)} \le t \sigma_n.$$

On the other hand, $\|u(y)\|_{H_0^1(\Omega)} = t \|u(ty)\|_{H_0^1(\Omega)} \ge tc_f$ by framing (2.23), which proves Theorem 2.8 with $\varepsilon_n = \sigma_n/c_f$. $\square$

The above theorem tells us that we can achieve contrast-independent approximation in relative error. It is however still unsatisfactory from two perspectives:

1. It does not describe the rate of decay of $\varepsilon_n$ as the reduced dimension $n$ grows. In practice, one would like to construct reduced spaces $V_n$ such that this decay is fast, similar to the exponential decay obtained under UEA.
<figure><img src="image_3.png" /><figcaption>Figure 1: Partition of $[0, 1]^d$ by the inverse rectangles $R_\ell^{-1}$ in the case $d=2$.</figcaption></figure>

2. The approximation property is expressed in terms of the orthogonal projection $P_{V_n}$. In applications to forward modeling, we approximate the solution $u(y)$ in the space $V_n$ by the Galerkin projection $P_{V_n}^y u(y)$. We thus wish for uniform estimates also for such approximations.

These two problems are treated in §3 and §4 respectively.

### 3 Approximation rates

Our construction of efficient reduced model spaces is based on a certain partitioning of the parameter domain $\tilde{Y}'$ associated to the manifold $\tilde{B}$. To any $\ell = (\ell_1, \dots, \ell_d) \in \mathbb{N}_0^d$ we associate the dyadic rectangle

$$R_{\ell} = [2^{\ell_1}, 2^{\ell_1+1}] \times \dots \times [2^{\ell_d}, 2^{\ell_d+1}], \quad (3.28)$$

For a positive integer $L$ to be fixed further, we modify the definition of $R_\ell$ by replacing the interval $[2^{\ell_j}, 2^{\ell_j+1}]$ by $[2^{\ell_j}, \infty]$ when $\ell_j = L$ for some $j$. This leads to the partition

$$\tilde{Y}' = \bigcup_{\ell \in \{0, \dots, L\}^d} R_{\ell}. \quad (3.29)$$

This partition is best visualized in the inverse parameter domain by setting

$$z = (z_1, \dots, z_d) := (y_1^{-1}, \dots, y_d^{-1}) \in [0, 1]^d. \quad (3.30)$$

Then, the inverse rectangles $R_\ell^{-1}$ split the unit cube, as shown on Figure 1. In particular, the rectangles touching the axes correspond to rectangles $R_\ell$ of infinite size.
We build reduced model spaces through a piecewise polynomial approximation over this partition. In other words, for each $\ell \in \{0, \dots, L\}^d$, we use different polynomials

$$u_{\ell,k}(y) = \sum_{|\nu| \le k} u_{\ell,\nu} y^{\nu},$$

of total degree $k$ for approximating $u(y)$ when $y \in R_\ell$, leading to a family of local reduced model spaces

$$V_{\ell,k} = \operatorname{span}\{u_{\ell,\nu} : |\nu| \le k\}, \qquad (3.31)$$

that can be either used individually when approximating $u(y)$ if the rectangle $R_\ell$ containing $y$ is known, or summed up in order to obtain a global reduced model space.

In this section we show that this construction yields exponential convergence rates in (1.15), similar to those obtained under a Uniform Ellipticity Assumption. This requires a proper tuning between the total polynomial degree $k$ and the integer $L$ that determines the size of the partition. In the study of local polynomial approximation, we treat separately the inner rectangles for which $\ell \in \{0, \dots, L-1\}^d$ and the infinite rectangles for which one or several $\ell_j$ are equal to $L$. The estimates obtained in the latter case rely on the additional assumption that the partition has a geometry of disjoint inclusions.

## 3.1 Polynomial approximation on inner rectangles

Inner rectangles $R_\ell$ are particular cases of rectangles of the form

$$R = [a_1, 2a_1] \times \dots \times [a_d, 2a_d], \qquad (3.32)$$

for some $a_j \ge 1$. The following lemma, adapted from [10], shows that one can approximate the parameter to solution map in the $\| \cdot \|_y$ and $\| \cdot \|_{H_0^1}$ norms on such rectangles, with a rate that decreases exponentially in the total polynomial degree.

**Lemma 3.1.** Let $R$ be any rectangle of the form (3.32). Then, for each $k \ge 0$, there exists functions $u_\nu \in H_0^1(\Omega)$ such that

$$\|u(y) - \sum_{|\nu| \le k} u_\nu y^\nu\|_{y} \le C3^{-k}, \quad y \in R, \qquad (3.33)$$

where $C := \frac{1}{\sqrt{3}}C_f$, and

$$\|u(y) - \sum_{|\nu| \le k} u_\nu y^\nu\|_{H_0^1} \le C3^{-k}, \quad y \in R, \qquad (3.34)$$

where $C := \frac{1}{\sqrt{6}}C_f$.
Proof. The exponential rate is established in [10] for a single parameter domain with uniform ellipticity assumption. Here the difficulty lies in the fact that we want the same estimate for all parametric rectangles R and thus without control on the uniform ellipticity. Still the technique of proof, based on power series, is similar.

The elliptic equation $-\mathrm{div}(a(y)u(y)) = f$ may be written in operator form

$$A_y u(y) = f,$$

where the invertible operator $A_y : H_0^1(\Omega) \to H^{-1}(\Omega)$ is defined by

$$\langle A_y v, w \rangle_{H^{-1}, H_0^1} := \int a(y) \nabla v \cdot \nabla w \, dx = \langle v, w \rangle_y.$$

We introduce

$$\bar{y} := \frac{3}{2}(a_1, \dots, a_d),$$

the center of the rectangle, and write any $y \in R$ as

$$y = \bar{y} + \tilde{y},$$

where the components $\tilde{y}_j$ of $\tilde{y}$ vary in $[-a_j/2, a_j/2]$. We may write $A_y = A_{\bar{y}} + \sum_{j=1}^d \tilde{y}_j A_j$, where the operators $A_j : H_0^1(\Omega) \to H^{-1}(\Omega)$ are defined by

$$\langle A_j v, w \rangle_{H^{-1}, H_0^1} := \int_{\Omega_j} \nabla v \cdot \nabla w \, dx.$$

This allows us to rewrite the equation as

$$(I + B(\tilde{y}))u(y) = g,$$

where $g := A_{\bar{y}}^{-1} f \in H_0^1(\Omega)$ and $B(\tilde{y}) = \sum_{j=1}^d \tilde{y}_j A_{\bar{y}}^{-1} A_j$ acts in $H_0^1(\Omega)$. We then observe that

$$\langle B(\tilde{y})v, w \rangle_{\bar{y}} = \langle A_{\bar{y}} B(\tilde{y})v, w \rangle_{H^{-1}, H_0^1} = \sum_{j=1}^d \tilde{y}_j \langle A_j v, w \rangle_{H^{-1}, H_0^1} = \sum_{j=1}^d \tilde{y}_j \int_{\Omega_j} \nabla v \cdot \nabla w \, dx,$$

and therefore, since $|\tilde{y}_j| \le \frac{1}{3}\bar{y}_j$,

$$|\langle B(\tilde{y})v, w \rangle_{\bar{y}}| \le \frac{1}{3} \sum_{j=1}^{d} \bar{y}_j \left| \int_{\Omega_j} \nabla v \cdot \nabla w \, dx \right| \le \frac{1}{3} \|v\|_{\bar{y}} \|w\|_{\bar{y}},$$

which shows that $\|B(\tilde{y})\|_{\bar{y} \to \bar{y}} \le \frac{1}{3}$. We may thus approximate $(I + B(\tilde{y}))^{-1}$ by the partial Neumann series

$$\sum_{l=0}^{k} (-1)^l B(\tilde{y})^l,$$
which is a polynomial in $\tilde{y}$ of total degree $k$. The corresponding polynomial approximation to $u(y)$ is given by

$$N_k u(y) = \sum_{l=0}^{k} (-1)^l B(\tilde{y})^l g = \sum_{l=0}^{k} (-1)^l \left( \sum_{j=1}^{d} \tilde{y}_j A_{\tilde{y}}^{-1} A_j \right)^l g = \sum_{|\nu| \le k} v_\nu \tilde{y}^\nu,$$

and coincides with the truncated power series of $\tilde{u}(\tilde{y}) := u(\bar{y} + \tilde{y})$ at $\tilde{y} = 0$, that is,

$$v_\nu := \frac{1}{\nu!} \partial^\nu u(\bar{y}), \quad \nu! := \prod \nu_j!.$$

It can be rewritten in the form

$$N_k u(y) = \sum_{|\nu| \le k} u_\nu y^\nu.$$

One has

$$\|u(y) - N_k u(y)\|_{\bar{y}} \le \sum_{l>k} \|B(\tilde{y})^l g\|_{\bar{y}} \le \left(\sum_{l>k} 3^{-l}\right) \|A_{\bar{y}}^{-1} f\|_{\bar{y}} = \frac{3^{-k}}{2} \|A_{\bar{y}}^{-1} f\|_{\bar{y}},$$

and

$$\|A_{\bar{y}}^{-1} f\|_{\bar{y}}^{2} = \langle A_{\bar{y}} A_{\bar{y}}^{-1} f, A_{\bar{y}}^{-1} f \rangle_{H^{-1}, H_{0}^{1}} = \langle f, u(\bar{y}) \rangle_{H^{-1}, H_{0}^{1}} \le C_{f} \|u(\bar{y})\|_{H_{0}^{1}} \le C_{f}^{2},$$

where the last inequality follows from Lax-Milgram estimate since $a(\bar{y}) \ge 1$. This proves the estimate

$$\|u(y) - \sum_{|\nu| \le k} u_\nu y^\nu\|_{\bar{y}} \le C 3^{-k}, \quad y \in \mathbb{R}, \qquad (3.35)$$

with $C := \frac{1}{2}C_f$. Using the inequalities

$$\|v\|_y^2 \le \frac{4}{3} \|v\|_{\bar{y}}^2, \quad v \in H_0^1(\Omega), y \in \mathbb{R},$$

and

$$\|v\|_{H_0^1}^2 \le \frac{2}{3} \|v\|_{\bar{y}}^2, \quad v \in H_0^1(\Omega),$$

we obtain the estimate (3.33) and (3.34) with the modified multiplicative constants. $\square$

**Remark 3.2.** The above lemma shows that the set $M_R := \{u(y) : y \in \mathbb{R}\}$ can be approximated with accuracy $C3^{-k}$ by the space

$$V_R := \operatorname{span}\{u_\nu : |\nu| \le k\}. \qquad (3.36)$$

The dimension of $V_R$ is at most $\binom{k+d}{d}$, however, as noticed in [10], it can in fact be seen that

$$\dim(V_R) \le \binom{k+d-1}{d-1}. \qquad (3.37)$$
This stems from the fact that the operators defined in the above proof satisfy the
dependency relation

$$
A_{\bar{y}} = \sum_{j=1}^{d} \bar{y}_j A_j,
$$

and therefore, one can rewrite $A_y$ as

$$
A_y := (1 + \tilde{y}_d / \bar{y}_d) A_{\bar{y}} + \sum_{j=1}^{d-1} (\tilde{y}_j - \tilde{y}_d \bar{y}_j / \bar{y}_d) A_j.
$$

Using this form, the partial Neumann sum $N_k u(y)$ has at most $\binom{k+d-1}{d-1}$ independent terms.

We shall also make use of the following adaptation of the above lemma to the approximation of the limit solution map $y_{S^c} \mapsto u_S(y_{S^c})$, defined by (2.21). Its proof is an immediate adaptation of the previous one and is therefore omitted.

**Lemma 3.3.** Let $S \subset \{1, \dots, d\}$, and for some $a_j \ge 1$, let $R$ be a rectangle of the form

$$
R = \prod_{j \in S^c} [a_j, 2a_j]. \qquad (3.38)
$$

Then, there exists functions $u_\nu \in V_S$ such that

$$
\| u_S(y_{S^c}) - \sum_{|\nu| \le k} u_\nu y_{S^c} \|_{y_{S^c}} \le C 3^{-k}, \quad y_{S^c} \in R, \quad (3.39)
$$

where $C := \frac{1}{\sqrt{3}} C_f$, and

$$
\| u_S(y_{S^c}) - \sum_{|\nu| \le k} u_\nu y_{S^c} \|_{H_0^1} \le C 3^{-k}, \quad y_{S^c} \in R, \quad (3.40)
$$

where $C := \frac{1}{\sqrt{6}} C_f$.

3.2 Polynomial approximation on infinite rectangles

We now consider the infinite rectangles $R_\ell$, corresponding to the $\ell$ such that
some of the $\ell_j$ equal $L$. We define

$$
S := \{j : \ell_j = L\}, \tag{3.41}
$$

the set of such indices. When $y \in R_\ell$, we thus have

$y_j \ge 2^L, \quad j \in S,$

and so $u(y)$ should be close to $u_S(y_{S^c})$ as $L$ is large. On the other hand $y_{S^c}$
belongs to a rectangle of the form

$$
R_{\ell_{S^c}} = \prod_{j \in S^c} [2^{\ell_j}, 2^{\ell_j+1}].
$$
Therefore, by Lemma 3.3, we can approximate $u_S(y_{S^c})$ by a polynomial of total degree $k$ in these restricted variables.

In order to conclude that this polynomial is a good approximation to $u(y)$ on $R_\ell$, we need a quantitative estimate on the convergence of $u(y)$ towards $u_S(y_{S^c})$. Let us observe that since

$$ \sum_{j=1}^{d} y_j \int_{\Omega_j} \nabla u(y) \cdot \nabla v \, dx = \langle f, v \rangle_{H^{-1}, H_0^1} = \sum_{j \in S^c} y_j \int_{\Omega_j} \nabla u_S(y_{S^c}) \cdot \nabla v \, dx, \quad v \in V_S, $$

the function $u_S(y_{S^c})$ coincides with the orthogonal projection of $u(y)$ onto $V_S$ for the $y$-norm, as well as for the $y_{S^c}$-norm:

$$ u_S(y_{S^c}) = P_{V_S}^y u(y) = P_{V_S}^{y_{S^c}} u(y). \quad (3.42) $$

In addition, with

$$ \Omega_S := \bigcup_{j \in S} \Omega_j, \quad (3.43) $$

we have

$$ 2^L \|\nabla u(y)\|_{L^2(\Omega_S)}^2 \le \sum_{j \in S} y_j \int_{\Omega_j} |\nabla u(y)|^2 dx \le \langle f, u(y) \rangle_{H^{-1}, H_0^1} \le C_f^2, $$

since $\|u(y)\|_{H_0^1} \le C_f$, and therefore, since $\nabla u_S(y_{S^c}) = 0$ on $\Omega_S$, we find that

$$ \|\nabla u(y) - \nabla u_S(y_{S^c})\|_{L^2(\Omega_S)} \le C_f 2^{-L/2}. \quad (3.44) $$

Our objective is to obtain a similar error bound on the remaining domains $\Omega_j$ for $j \in S^c$. This turns out to be feasible, with an even better rate $2^{-L}$, when making certain geometric assumptions on the partition of the domain $\Omega$.

**Definition 3.4.** We say that $\{\Omega_1, \dots, \Omega_d\}$ is a Lipschitz partition if and only if for any subset $T \subset \{1, \dots, d\}$, the domain $\Omega_T = \bigcup_{j \in T} \Omega_j$ has Lipschitz boundaries.

<figure><img src="image_4.png" /><figcaption>Figure 2: A Lipschitz partition of $\Omega$ (left) and a counter-example (right) since $\Omega_1 \cup \Omega_4$ is not Lipschitz.</figcaption></figure>

Note that such a property is stronger than just saying that each domain is Lipschitz, see Figure 2 (right) for a counter-example. In a Lipschitz partition,
all subdomains $\Omega_j$ are Lipschitz, and the common boundary between two subdomains is either empty or a $(n-1)$-dimensional surface, as illustrated on Figure 2 (left). In particular, it is easily checked that partitions consisting of a background domain and well separated subdomains that have Lipschitz boundaries fall in this category. Similar to the $\Omega_T$, the individual $\Omega_j$ could have several connected components, that should then be well separated. Here by “well separated”, we mean that $\delta$-neighbourhoods of the subdomains remain disjoints for some $\delta > 0$.

For the inner domains $\Omega_T$ such that $\partial\Omega_T \cap \partial\Omega = \emptyset$, the classical Stein's extension theorem [30] guarantees the existence of continuous extension operators

$$E_T : H^1(\Omega_T) \to H^1(\Omega),$$

that satisfy $(E_T v)|_{\Omega_T} = v$ for all $v \in H^1(\Omega_T)$. We refer to chapter 5 of [1] for a relatively simple construction of the extension operator $E_j$ by local reflection after using a partitioning of unity along the boundary of $\Omega_T$ and local transformations mapping the boundary to the hyperplane $\mathbb{R}^{n-1}$.

For the domains $\Omega_T$ touching the boundary $\partial\Omega$, these operators are modified in order to take into account the homogeneous boundary condition, and we refer to [36] for such adaptations. Here, the relevant space is

$$\tilde{H}^1(\Omega_T) := R_T(H_0^1(\Omega)), \qquad (3.45)$$

where $R_T$ is the restriction to $\Omega_T$, over which $v \mapsto \|\nabla v\|_{L^2(\Omega_T)}$ is equivalent to the $H^1$ norm by Poincaré inequality. Then, there exists a continuous extension operator

$$E_T : \tilde{H}^1(\Omega_T) \to H_0^1(\Omega).$$

Note that the norm of all these operators depends on the geometry of the partition. These operators are instrumental in proving the following convergence estimate.

**Lemma 3.5.** Assume that $\{\Omega_1, \dots, \Omega_d\}$ is a Lipschitz partition of $\Omega$. Then there exists a constant $C_0$ that only depends on the geometry of the partition such that for any $S \subset \{1, \dots, d\}$ and $y = (y_S, y_{S^c}) \in Y'$, one has

$$\|u(y) - u_S(y_{S^c})\|_{H_0^1} \le C_0 C_f \max_{j \in S} y_j^{-1}. \qquad (3.46)$$

In particular, for the infinite rectangle $R_\ell$,

$$\|u(y) - u_S(y_{S^c})\|_{H_0^1} \le C_0 C_f 2^{-L}, \quad y \in R_\ell, \qquad (3.47)$$

with $S$ defined by (3.41).

*Proof.* We first note that it suffices to prove (3.46) in the particular case where the largest $y_j$ are those for which $j \in S$. Indeed, if this is not the case, we use the decomposition

$$u(y) - u_S(y_{S^c}) = (u(y) - u_{S'}(y_{S^c})) - (u(y') - u_{S'}(y_{S^c})) + (u(y') - u_S(y_{S^c})),$$
with $S' = \{i : y_i \ge \min_{j \in S} y_j\}$ and $y'$ defined by $y'_j = \max_{i=1,...,d} y_i$ if $j \in S$, $y'_j = y_j$ otherwise, so that each term falls in this particular case and will be bounded in $H_0^1$ norm by $C_0C_f \max_{j \in S} y_j^{-1}$. This leads to the same estimate (3.46) up to a factor 3 in constant $C_0$. In addition, up to reordering the subdomains $\Omega_j$, we may assume $y_1 \ge \dots \ge y_d$ and therefore $S = \{1, \dots, |S|\}$.

Fix $j \ge |S|$, and denote $u = u(y)$ and $u_S = u_S(y_{S^c})$ for simplicity. We define the Lipschitz domain $\Omega^j = \bar{\Omega}_1 \cup \dots \cup \bar{\Omega}_j$, remarking that

$$ \Omega_S = \bigcup_{j \in S} \Omega_j = \Omega^{|S|}. $$

Poincaré's inequality ensures that there exists a function $c$ on $\Omega^j$, constant on any connected component of $\Omega^j$, and null on $\partial\Omega \cap \Omega^j$, such that

$$ \|u - u_S - c\|_{H^1(\Omega^j)} \le C_P \|\nabla(u - u_S)\|_{L^2(\Omega^j)}, $$

with $C_P$ the maximal Poincaré constant of all unions of subdomains from the partition. Moreover, there is an extension $v \in H_0^1(\Omega)$ of $u - u_S - c \in \tilde{H}^1(\Omega^j)$ such that

$$ \|v\|_{H_0^1(\Omega)} \le C_E \|u - u_S - c\|_{H^1(\Omega^j)} \le C_E C_P \|\nabla(u - u_S)\|_{L^2(\Omega^j)}, $$

with $C_E$ the maximal norm of all extension operators $E_T$, $T \subset \{1, \dots, d\}$.

As $u - u_S - v = c$ on $\Omega_S \subset \Omega^j$, the function $u - u_S - v$ is in $V_S$, and therefore orthogonal to $u - u_S = u - P_{V_S}^y u$ for the $\| \cdot \|_y$ norm:

$$ \begin{aligned} 0 &= \langle u - u_S, u - u_S - v \rangle_y \\ &= \sum_{i=1}^{d} y_i \int_{\Omega_i} |\nabla(u - u_S)|^2 - \sum_{i=1}^{d} y_i \int_{\Omega_i} \nabla(u - u_S) \cdot \nabla v \\ &= \sum_{i>j} y_i \int_{\Omega_i} |\nabla(u - u_S)|^2 - \sum_{i>j} y_i \int_{\Omega_i} \nabla(u - u_S) \cdot \nabla v \end{aligned} $$

since $\nabla v = \nabla(u - u_S)$ on $\Omega^j$. In particular, we obtain

$$ \begin{aligned} \|\nabla(u - u_S)\|_{L^2(\Omega_{j+1})}^2 &\le \sum_{i>j} y_i \int_{\Omega_i} |\nabla(u - u_S)|^2 \\ &\le y_{j+1} \int_{\Omega \setminus \Omega^j} |\nabla(u - u_S) \cdot \nabla v| \\ &\le y_{j+1} \|u - u_S\|_{H_0^1(\Omega)} \|v\|_{H_0^1(\Omega)} \\ &\le y_{j+1} \|u - u_S\|_{H_0^1(\Omega)} C_P C_E \|\nabla(u - u_S)\|_{L^2(\Omega^j)}, \end{aligned} $$

and therefore

$$ \|\nabla(u - u_S)\|_{L^2(\Omega^{j+1})}^2 \le (1 + C_P C_E) \|\nabla(u - u_S)\|_{L^2(\Omega)} \|\nabla(u - u_S)\|_{L^2(\Omega^j)}. $$
Applying this inequality inductively for $j = d-1, \dots, d-k$, we get

$$ \|\nabla(u - u_S)\|_{L^2(\Omega)} \le (1 + C_P C_E)^{2^k - 1} \|\nabla(u - u_S)\|_{L^2(\Omega^{d-k})}, $$

for any $k = 1, \dots, d - |S|$. For $k = d - |S|$, this results in the bound

$$ \|\nabla(u - u_S)\|_{L^2(\Omega)}^2 \le C_0 \|\nabla(u - u_S)\|_{L^2(\Omega_S)}^2 = C_0 \|\nabla u\|_{L^2(\Omega_S)}^2, \quad (3.48) $$

for any non-empty $S$, with $C_0 = (1 + C_P C_E)^{2^{d-1}}$.

We now write

$$ \begin{aligned} \left(\min_{i \in S} y_i\right) \|\nabla(u - u_S)\|_{L^2(\Omega_S)}^2 &\le \|u - u_S\|_y^2 = \langle u, u - 2u_S \rangle_y + \langle u_S, u_S \rangle_{y_{S^c}} \\ &= \langle f, u - u_S \rangle_{H^{-1}, H_0^1} \le C_f \|\nabla(u - u_S)\|_{L^2(\Omega)}, \end{aligned} $$

which, combined to the previous estimate, gives

$$ \|u - u_S\|_{H_0^1} = \|\nabla(u - u_S)\|_{L^2(\Omega)} \le C_0 C_f \max_{i \in S} y_i^{-1}, $$

therefore proving (3.46). For (3.47), we simply notice that $\max_{j \in S} y_j^{-1} \le 2^{-L}$ for $y \in Y' \cap R_\ell$, and use a continuity argument when $y$ takes infinite values. $\square$

Combining the estimate (3.47) from the above lemma with (3.40) from Lemma 3.3, we obtain the following estimate for polynomial approximation on an infinite rectangle $R_\ell$:

$$ \left\| u(y) - \sum_{|\nu| \le k} u_\nu y_{S^c}^\nu \right\|_{H_0^1} \le \frac{C_f}{\sqrt{6}} 3^{-k} + C_0 C_f 2^{-L}, \quad y \in R_\ell, \quad (3.49) $$

where $C_0$ is the constant in (3.47). This estimate hints how the level $L$ in the partition should be tuned to the total polynomial degree $k$, so that the two contributions in the above estimate are of the same order.

**Remark 3.6.** Note that the constant $C_0 = (1 + C_P C_E)^{2^{d-1}}$ becomes prohibitive even for moderate values of $d$. However, under more restrictive geometric assumptions, for instance if the subdomains $\bar{\Omega}_2, \dots, \bar{\Omega}_d$ are disjoint inclusions in a background $\Omega_1$, better bounds can be obtained, with a constant $C_0$ that does not suffer a similar curse of dimensionality, by replacing the induction in the proof by a two-step procedure, consisting of extensions first from the high-diffusivity inclusions to the background, and then to the whole domain $\Omega$.

## 3.3 Approximation rates and n-widths

We are now in position to establish an approximation result for the reduced model spaces. For this purpose, we fix the smallest level $L = L_k \ge 1$ such that

$$ C_0 C_f 2^{-L} \le \frac{C_f}{\sqrt{3}} 3^{-k}. $$
In particular L scales linearly with k, with the bound $\alpha k + \beta \le L_k \le \alpha k + \gamma$, where

$$ \alpha := \frac{\ln 3}{\ln 2}, \quad \beta := \frac{\ln(\sqrt{3}C_0)}{\ln 2}, \quad \gamma := \frac{\ln(2\sqrt{3}C_0)}{\ln 2}. \qquad (3.50) $$

Then, the polynomial approximation estimates (3.34) and (3.49) show that for each $\ell \in \{0, \dots, L_k\}^d$, there exist functions $u_{\ell,\nu} \in H_0^1(\Omega)$ such that

$$ \|u(y) - \sum_{|\nu| \le k} u_{\ell,\nu} y^\nu\|_{H_0^1} \le \left( \frac{C_f}{\sqrt{6}} + \frac{C_f}{\sqrt{3}} \right) 3^{-k} \le C_f 3^{-k}, \quad y \in R_\ell. $$

Note that in the case of an infinite rectangle $R_\ell$, the $u_{\ell,\nu}$ are non trivial only for monomials of the form $y_{S_c}^\nu$ and they belong to $V_S$, where $S := \{j : \ell_j = L_k\}$.

Thus the solutions $u(y)$ for $y \in R_\ell$ are approximated with accuracy $C_f 3^{-k}$ in the space

$$ V_{\ell,k} := \operatorname{span}\{u_{\ell,\nu} : |\nu| \le k\}, $$

which in view of Remark 3.2 has dimension at most $\binom{k+d-1}{d-1}$.

Note also that approximating the reduced manifold $\mathcal{N}$ defined in (2.22) requires a smaller subset of rectangles, since

$$ \{y \in \tilde{Y}' : \min y_j = 1\} \subset \bigcup_{\ell \in E_k} R_\ell, \quad E_k := \{0, \dots, L_k\}^d \setminus \{1, \dots, L_k\}^d. $$

We thus introduce the reduced model space

$$ V_n := \bigoplus_{\ell \in E_k} V_{\ell,k}, \quad n = \dim(V_n) \le \#(E_k) \binom{k+d-1}{d-1}, \qquad (3.51) $$

and find that

$$ \|u(y) - P_{V_n} u(y)\|_{H_0^1} \le C_f 3^{-k}, \qquad (3.52) $$

for all $y \in \tilde{Y}'$ such that $\min y_j = 1$. In view of (3.50), there exists a constant C that depends on d and $C_0$, such that

$$ n \le ((L_k + 1)^d - L_k^d) \binom{k+d-1}{d-1} \le C(k+1)^{2d-2}. \qquad (3.53) $$

This leads to the following approximation theorem.

**Theorem 3.7.** Assume that the partition has the geometry of disjoint inclusions. The reduced basis space $V_n$ defined in (3.51) then satisfies

$$ \|u(y) - P_{V_n} u(y)\|_{H_0^1} \le C \exp\left(-cn^{\frac{1}{2d-2}}\right), \qquad (3.54) $$

for all $y \in \tilde{Y}' = [1, \infty]^d$ such that $\min y_j = 1$. The Kolmogorov n-width (1.2) of the reduced manifold $\mathcal{N}$ satisfies

$$ d_n(\mathcal{N})_{H_0^1} \le C \exp\left(-cn^{\frac{1}{2d-2}}\right). \qquad (3.55) $$
Over the full manifold $\overline{\mathcal{M}}$, one has the estimate in relative error

$$ \|u(y) - P_{V_n} u(y)\|_{H_0^1} \le C \exp\left(-cn^{\frac{1}{2d-2}}\right) \|u(y)\|_{H_0^1}, \quad (3.56) $$

for all $y \in \tilde{Y} = [0, \infty]^d$. The positive constants $c$ and $C$ only depend on $d$, $C_f$, and on the geometry of the partition through the constant $C_0$.

*Proof.* The estimate (3.54) follows directly by combining (3.52) and (3.53), and (3.55) is an immediate consequence. We then derive (3.56) by using the homogeneity property (1.13) and the lower inequality in (2.25), similar to the proof of (2.27) in Theorem 2.8. $\square$

**Remark 3.8.** In the above construction of $V_n$, the dimension $n$ only takes the values $n_k := \#(E_k)^{\binom{k+d-1}{d-1}}$ for $k \ge 0$. However it is easily seen that if we set $V_n = V_{n_k}$ for $n_k \le n < n_{k+1}$, then all the estimates in the above theorem remain valid up to a change in the constants $(c, C)$.

**Remark 3.9.** Note that the union of the $V_{\ell,k}$ for $\ell \in E_k$ would suffice to approximate $\mathcal{N}$ with uniform accuracy $C_f 3^{-k}$, their sum $V_n$ is an overkill. When $y$ is known, for example in forward modeling, it is therefore possible to first identify the proper space $V_{\ell,k}$ associated to the rectangle $R_\ell$ that contains $y$, and build the approximation to $u(y)$ from this space. This nonlinear reduced modeling strategy has been studied in [14] with similar local polynomial approximation under UEA, and in [26, 25, 23] with local reduced basis. The natural benchmark is given by the notion of library width introduced in [31], that is defined for any compact set $\mathcal{K}$ in a Banach space $V$ as

$$ d_{n,N}(\mathcal{K})_V := \inf_{\#(\mathcal{L}_n) \le N} \sup_{u \in \mathcal{K}} \min_{V_n \in \mathcal{L}_n} \min_{v \in V_n} \|u - v\|_V, \quad (3.57) $$

where the first infimum is taken over all libraries $\mathcal{L}_n$ of $n$-dimensional spaces with cardinality at most $N$. Our results thus show that

$$ d_{n,N}(\mathcal{N})_{H_0^1} \le C_f 3^{-k} \sim C \exp(-cn^{\frac{1}{d}}), \quad n := \binom{k+d-1}{d-1}, \quad N = (L_k+1)^d - L_k^d. $$

Note that the above sub-exponential rate can be misleading due to fact that the constant $c$ has a hidden dependence in $d$. As an example, up to the constant $C_f$, we find that taking $k = 4, 7, 9$ leads to error bounds $3^{-k}$ of order $10^{-2}, 10^{-3}, 10^{-4}$, with $n = 15, 36, 55$ for $d = 3$, and $n = 35, 120, 220$ for $d = 4$, which is far better than the value of $\exp(-n^{\frac{1}{d}})$.

**Remark 3.10.** In view of the results from [12] and [19], we are ensured that a proper selection of reduced basis elements in the manifold $\mathcal{N}$ should generate spaces $V_n$ that perform at least with the same exponential rates as those achieved by the spaces $V_n$ in Theorem 3.7. As explained in the introduction, reduced basis spaces may perform significantly better than reduced model spaces based on polynomial or piecewise polynomial approximation. This occurs in particular
when the polynomial coefficients have certain linear dependency, as established in [10] for the elliptic problem with piecewise constant coefficients in the low contrast regime, and recalled in Remark 3.2. There, it is shown that the rate $O(\exp(-cn^{1/d}))$ is at least improved to $O(\exp(-cn^{d-1}/d))$ and that further improvements in the rate may result from certain symmetry properties of the domain partition, however not circumventing the curse of dimensionality. While we do not pursue this analysis in the present high contrast setting, we expect similar results to hold.

# 4 Forward modeling and inverse problems

## 4.1 Galerkin projection

In the context of forward modeling, the reduced model space $V_n$ is used to approximate the parameter to solution map, by a map

$$y \mapsto u_n(y) \in V_n,$$

computed through the Galerkin method: $u_n(y) \in V_n$ is such that

$$\sum_{j=1}^{d} y_j \int_{\Omega_j} \nabla u_n(y) \cdot \nabla v \, dx = \langle f, v \rangle_{H^{-1}, H_0^1}, \quad v \in V_n.$$

Therefore $\langle u_n(y), v \rangle_y = \langle u(y), v \rangle_y$, that is

$$u_n(y) = P_{V_n}^y u(y),$$

where $P_{V_n}^y$ is the projection onto $V_n$ with respect to norm $\| \cdot \|_y$.

Hence, one would like to derive estimates on $\|u(y) - P_{V_n}^y u(y)\|_{H_0^1}$ in place of the estimates on $\|u(y) - P_{V_n} u(y)\|_{H_0^1}$ that we have obtained so far, since $P_{V_n} u(y)$ is not practically accessible. As explained in the introduction, we cannot be satisfied with combining the latter estimates with the bound

$$\|u(y) - P_{V_n}^y u(y)\|_{H_0^1} \le \kappa(y)^{1/2} \|u(y) - P_{V_n} u(y)\|_{H_0^1}$$

derived from Cea's lemma, since the multiplicative constant $\kappa(y)$ from (1.9) is not uniformly bounded over the manifolds $M$, $B$ or $N$. Here, we shall employ another approach to derive the same rates of convergence for $\|u(y) - P_{V_n}^y u(y)\|_{H_0^1}$.

One first observation is that in order for Galerkin projection $P_{V_n}^y$ onto a reduced model space $V_n$ to satisfy a convergence bound in relative error, it is critical that this space contains some functions from the limit spaces $V_S$. This is expressed by the following result.

**Proposition 4.1.** Assume that there exists $S \subset \{1, \dots, d\}$ such that $V_n \cap V_S = \{0\}$. Then for any $C \in [0, 1[$, there exists $y \in Y'$ such that

$$\|u(y) - P_{V_n}^y u(y)\|_{H_0^1} \ge C \|u(y)\|_{H_0^1}. \qquad (4.58)$$
*Proof.* Since $V_n \cap V_S = \{0\}$, the quantity $\|\nabla v\|_{L^2(\Omega_S)}$ is a norm on $V_n$ and one can define

$$ \alpha = \min_{v \in V_n} \frac{\|\nabla v\|_{L^2(\Omega_S)}}{\|v\|_{H_0^1}} > 0. $$

For any $\varepsilon > 0$, take $y_j = \varepsilon^{-2}$ for $j \in S$ and $y_j = 1$ for $j \in S^c$. Then, for $v = P_{V_n}^y u(y)$,

$$ \frac{\alpha}{\varepsilon} \|v\|_{H_0^1} \le \frac{1}{\varepsilon} \|\nabla v\|_{L^2(\Omega_S)} \le \|v\|_y \le \|u(y)\|_y \le C_f \le \frac{C_f}{c_f} \|u(y)\|_{H_0^1}, $$

where we have used the framings (2.25) and (2.26). Therefore, taking $\varepsilon = \frac{c_f}{C_f} \alpha(1-C)$ implies $\|v\|_{H_0^1} \le (1-C)\|u(y)\|_{H_0^1}$, and (4.58) follows. $\square$

However, in the construction of $V_n$ in §3, each space $V_{\ell,k}$ is a subset of $V_S$ for $S = \{j : \ell_j = L_k\}$. This prevents the phenomenon described in the previous proposition from occurring. Instead, we obtain similar convergence bounds as those obtained for $P_{V_n}$, as expressed in the following result.

**Theorem 4.2.** Assume that the partition of $\Omega$ has the geometry of disjoint inclusions. On the rectangles $R_\ell$ for $\ell \in \{0, \dots, L\}^d$, the following uniform convergence estimates hold:

$$ \|u(y) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1} \le \frac{C_f}{\sqrt{3}} 3^{-k}, \quad y \in R_\ell, \qquad (4.59) $$

if $\|\ell\|_\infty < L$, and

$$ \|u(y) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1} \le \frac{C_f}{\sqrt{3}} 3^{-k} + C_0 C_f 2^{-L}, \quad y \in R_\ell, \qquad (4.60) $$

if $\|\ell\|_\infty = L$. As a consequence, with $L = L_k$ and $V_n$ defined as in §3.3, one has the estimates

$$ \|u(y) - P_{V_n}^y u(y)\|_{H_0^1} \le C \exp \left( -cn^{\frac{1}{2d-2}} \right), \qquad (4.61) $$

for all $y \in \tilde{Y}'$ such that $\min y_j = 1$, and

$$ \|u(y) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1} \le C \exp \left( -cn^{1/(2d-2)} \right) \|u(y)\|_{H_0^1}, \qquad (4.62) $$

for all $y \in \tilde{Y}$, with constants $c$ and $C$ that only depend on $d$, $C_f$, and on the geometry of the partition through the constant $C_0$.

*Proof.* For bounded rectangles $R_\ell$ with $\|\ell\|_\infty < L$, we know from Lemma 3.1, and more precisely from (3.33), that

$$ \|u(y) - P_{V_{\ell,k}}^y u(y)\|_y = \min_{v \in V_{\ell,k}} \|u(y) - v\|_y \le \left\| \|u(y) - \sum_{|\nu| \le k} u_\nu y^\nu \|_y \le \frac{C_f}{\sqrt{3}} 3^{-k} \right\| $$
for any $y \in R_{\ell}$. Since all the $y_j$ are greater or equal to 1, one has $\|v\|_{H_0^1} \le \|v\|_y$ for all $y$ and therefore (4.59) follows.

For infinite rectangles $R_{\ell}$ such that $\|\ell\|_{\infty} = L$, we again introduce $S = \{j : \ell_j = L\}$. Then, using (3.47),

$$ \begin{aligned} \|u(y) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1} &\le \|u(y) - u_S(y_{S^c})\|_{H_0^1} + \|u_S(y_{S^c}) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1} \\ &\le C_0 C_f 2^{-L} + \|u_S(y_{S^c}) - P_{V_{\ell,k}}^y u(y)\|_{H_0^1}. \end{aligned} $$

Since $V_{\ell,k} \subset V_S$, we have

$$ P_{V_{\ell,k}}^y u(y) = P_{V_{\ell,k}}^y P_{V_S}^y u(y) = P_{V_{\ell,k}}^y u_S(y_{S^c}) = P_{V_{\ell,k}}^{y_{S^c}} u_S(y_{S^c}), $$

Similarly to the previous case, we apply (3.39) from Lemma 3.3:

$$ \|u_S(y_{S^c}) - P_{V_{\ell,k}}^y u_S(y_{S^c})\|_{H_0^1} \le \|u_S(y_{S^c}) - P_{V_{\ell,k}}^y u_S(y_{S^c})\|_y \le \frac{C_f}{\sqrt{3}} 3^{-k}, $$

and we thus obtain (4.60).

After taking $L = L_k$ and defining $V_n$ as the sum of the $V_{\ell,k}$ for $\ell \in E_k$, the derivation of (4.61) and (4.62) is exactly the same as for (3.54) and (3.56). $\square$

**Remark 4.3.** As in Remark 3.10, it is expected that the same rate of convergence is attained if $V_n$ is a reduced basis space generated by solutions $u(y^i)$, $i = 1, \dots, n$, as long as there are $O((k+d-1)/(d-1))$ samples $y^i$ in each rectangle, however with samples forced to be of the form $u_S(y_{S^c}^i) \in V_S$ in the case of infinite rectangles.

## 4.2 State and parameter estimation

The state estimation problem consists in retrieving the solution $\bar{u} = u(\bar{y})$ when the parameter $\bar{y}$ is unknown, and one observes $m$ linear measurements

$$ w_i = \ell_i(\bar{u}), \quad i = 1, \dots, m, $$

where the $\ell_i$ are continuous linear functional on the Hilbert space $V$ that contains the solution manifold. These linear functionals may thus be written in terms of Riesz representers

$$ \ell_i(v) = \langle \omega_i, v \rangle_V. $$

The Parametrized Background Data Weak (PBDW) method, introduced in [24] and further studied in [13], exploits the fact that all potential solutions are well approximated by reduced model spaces $V_n$. It is based on a simple recovery algorithm that consists in solving the problem

$$ \min_{u \in V_w} \min_{v \in V_n} \|u - v\|_V, \qquad (4.63) $$
where, for $w = (w_1, \dots, w_m) \in \mathbb{R}^m$,

$$V_w := \{u \in V : \ell_i(u) = w_i, i = 1, \dots, m\},$$

is the affine space of functions that agree with the measurements.

The analysis of this problem is governed by the quantity

$$\mu_n = \mu(V_n, W) := \sup_{v \in V_n} \frac{\|v\|_V}{\|P_W w\|_V}, \quad (4.64)$$

where $W := \operatorname{span}\{\omega_1, \dots, \omega_m\}$, which is finite if and only if $V_n \cap W^\perp = \{0\}$. Then, there exists a unique minimizing pair

$$(u^*, v^*) = (u^*(w), v^*(w)) \in V_w \times V_n$$

to (4.63), which satisfies the estimates

$$\|\bar{u} - v^*\|_V \le \mu_n \min_{v \in V_n} \|u - v\|_V, \quad (4.65)$$

and

$$\|\bar{u} - u^*\|_V \le \mu_n \min_{v \in V_n + (W \cap V_n^\perp)} \|u - v\|_V. \quad (4.66)$$

The computation of $(u^*, v^*)$ amounts to solving finite linear systems, and both solutions depend linearly on $w$.

Turning to our specific elliptic problem, and assuming that the $\ell_i$ belong to $H^{-1}(\Omega) = V'$ for $V = H_0^1(\Omega)$, we may apply the above PBDW method using the reduced basis spaces $V_n$ introduced in §3. As an immediate consequence of Theorem (3.7), we obtain a recovery estimate in relative error.

**Proposition 4.4.** Let $\bar{y} \in \tilde{Y}$ and $\bar{u} = u(\bar{y})$. Then both estimators $v^* \in V_n$ and $u^* \in V_w$ satisfy

$$\max\{\|\bar{u} - v^*\|_{H_0^1}, \|\bar{u} - u^*\|_{H_0^1}\} \le C\mu_n \exp\left(-cn^{\frac{1}{2d-2}}\right) \|\bar{u}\|_{H_0^1}. \quad (4.67)$$

The positive constants $c$ and $C$ only depend on $d$, $C_f$, and on the geometry of the partition through the constant $C_0$.

*Proof.* It follows readily by combining (3.56) applied to $y = \bar{y}$ with the recovery estimates (4.65) and (4.66). $\square$

We next turn to the problem of parameter estimation, namely recovering an approximation $y^*$ to $\bar{y}$ from the measurements $w$. In contrast to state estimation, this is a nonlinear inverse problem since the first mapping in

$$\bar{y} \mapsto \bar{u} \mapsto w$$

is typically nonlinear. One way of relaxing this problem into a linear one is by first using a recovery $u^*$ of the state $\bar{u}$, for example obtained by the PBDW method. One then defines $y^*$ as the minimizer over $\tilde{Y}$ of the residual

$$R(y) := \|\operatorname{div}(a(y)\nabla u^*) + f\|_{H^{-1}}.$$
This is a quadratic problem when $a(y)$ has an affine dependence in $y$, that can be solved by standard quadratic optimization methods. The rationale for this approach is the fact that

$$R(y) = \|A_y u^* - A_y u(y)\|_{H^{-1}} \sim \|u^* - u(y)\|_{H_0^1},$$

and therefore we should be close to finding the parameter $y$ that best explains the approximation $u^*$. Unfortunately, this approach is not much viable in the high-contrast regime since the equivalence $\|A_y v\|_{H^{-1}} \sim \|v\|_{H_0^1}$ has constants that are not uniform in $y$ and deteriorate with the level of contrast.

Instead, we propose a more specific approach that exploits the piecewise constant structure of $a(y)$, assuming that $V_n$ is a reduced space of the form

$$V_n = \operatorname{span}(u^1, \dots, u^n), \quad u^i = u(y^i),$$

for some properly selected parameter vectors

$$y^i = (y_1^i, \dots, y_d^i), \quad i = 1, \dots, n.$$

As mentioned, see Remark (3.10), these spaces satisfy the same exponential convergence bounds as the spaces constructed in §3.

The PBDW estimator $v^* = v^*(w) \in V_n$ thus has the form

$$v^* = \sum_{i=1}^{n} c_i u^i \in V_n$$

and satisfies a similar bound (4.67) as in the above proposition. Then, on the particular domain $\Omega_j$, one has

$$\frac{f}{\bar{y}_j} = -\Delta \bar{u}_{|\Omega_j} \approx -\sum_{i=1}^{n} c_i \Delta u^i = \sum_{i=1}^{n} c_i \frac{f}{y_j^i},$$

and therefore, a natural candidate for the parameter estimate is $y^* = (y_1^*, \dots, y_d^*)$ with

$$y_j^* := \left( \sum_{i=1}^{n} \frac{c_i}{y_j^i} \right)^{-1}. \qquad (4.68)$$

The following result gives a recovery bound in relative error for the inverse diffusivity.

**Proposition 4.5.** With the notation $1/y = (1/y_1, \dots, 1/y_d)$, the estimator $y^*$ defined by (4.68) satisfies the bound

$$\left\| \frac{1}{y^*} - \frac{1}{y} \right\|_{\infty} \le \frac{C_f}{c_f} C \mu_n \exp \left( -cn^{\frac{1}{2d-2}} \right) \left\| \frac{1}{y} \right\|_{\infty}, \qquad (4.69)$$

where $C_f$ and $c_f$ are as in (2.25), and the other constants as in (4.67).
Proof. For $1 \le j \le d$, take $\phi \in H_0^1(\Omega_j)$, then

$$
\begin{aligned}
\left| \frac{1}{y_j^*} - \frac{1}{\bar{y}_j} \right| |\langle f, \phi \rangle_{H^{-1}, H_0^1}| &= \left| \sum_{i=1}^n \frac{c_i}{y_j^i} \int_{\Omega_j} y_j^i \nabla u^i \cdot \nabla \phi \, dx - \frac{1}{\bar{y}_j} \int_{\Omega_j} \bar{y}_j \nabla \bar{u} \cdot \nabla \phi \, dx \right| \\
&= \left| \int_{\Omega_j} \nabla (v^* - \bar{u}) \cdot \nabla \phi \, dx \right| \\
&\le \|v^* - \bar{u}\|_{H_0^1(\Omega)} \|\phi\|_{H_0^1(\Omega_j)}.
\end{aligned}
$$

Optimizing over $\phi$ gives

$$
\left\| \frac{1}{y^*} - \frac{1}{\bar{y}} \right\|_{\infty} \le c_f^{-1} \|v^* - \bar{u}\|_{H_0^1},
$$

which combined with (4.67) gives

$$
\left\| \frac{1}{y^*} - \frac{1}{\bar{y}} \right\|_{\infty} \le c_f^{-1} C \mu_n \exp \left( -cn^{\frac{1}{2d-2}} \right) \|\bar{u}\|_{H_0^1}.
$$

Using the Lax-Milgram estimate

$$
\| \bar{u} \|_{H_0^1} \le C_f \| \frac{1}{\bar{y}} \|_{\infty},
$$

we reach (4.69).

<figure><img src="image_5.png" /></figure>

**Remark 4.6.** The bound (4.69) is not entirely satisfactory since the approximation error on $\bar{y}_j$ remains high when $\bar{y} \in \mathcal{N}$ with $\bar{y}_j \gg 1$. We do not know if a bound of the form

$$
\left| \frac{1}{y_j^*} - \frac{1}{\bar{y}_j} \right| \le \frac{\varepsilon_n}{\bar{y}_j}, \quad 1 \le j \le d,
$$

which would imply $|y_j^* - \bar{y}_j| \le \varepsilon_n/(1 - \varepsilon_n) \bar{y}_j$, holds uniformly over $\mathcal{N}$ with $\varepsilon_n \underset{n \to +\infty}{\longrightarrow} 0$.

**5 Numerical illustration**

The base model that will be used all along the numerical illustrations is the
diffusion equation (1.4) with data $f = 1$ set on the two-dimensional square
$\Omega = [-1, 1]^2$ with homogeneous Dirichlet boundary conditions. We consider a
piece-wise constant diffusion coefficient

$a|_{\Omega_j} = y_j, \quad 1 \le j \le d,$

on a partition of $\Omega$ into 16 squares of quarter side-length.

As such this partition does not satisfy the geometrical assumption of “Lipschitz partition” that was critical in our analysis for the application of Lemma
<figure><img src="image_6.png" /><figcaption>Figure 3: Lipschitz partition of Ω.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 4: Non-lipschitz partition of Ω.</figcaption></figure>

3.5. Therefore we consider sub-partitions that comply to the assumptions, such as illustrated on Figure 3, which amounts to equate the parameters $y_j$ of squares belonging to the same sub-domain. This way we can consider that $y = (y_A, y_B, y_C, y_D)$ consists of four parameters, one per each subdomain.

The numerical results that we next present aim to illustrate the robustness to high-contrast of the reduced basis method, and discuss in addition the effect of parameter selection, higher parametric dimensions, and inclusions that are not satisfying the geometric assumption as exemplified on Figure 4.

We construct different reduced bases {$u^1, \dots, u^n$} of moderate dimension $1 \le n \le 15$, where

$$u^k = u(y^k),$$

for certain parameter selections $y^1, \dots, y^n$. Each reduced basis element $u^k$ is numerically computed by the Galerkin method in a background finite element space $V_h$ of dimension 6241.

The reduced basis spaces are thus subspaces of $V_h$, thus strictly speaking spaces $V_{n,h}$ depending on $n$ and on the meshsize $h$. In our numerical computation, we always assess the error

$$P_{V_h}^y u(y) - P_{V_{n,h}}^y u(y).$$

We noticed that for the considered values of $n = 1, \dots, 15$ the error curves do not vary much when further reducing the mesh size $h$. In fact they are already essentially the same when the dimension of $V_h$ is four times smaller. Therefore, for simplicity of the presentation, we still write

$$u(y) - P_{V_h}^y u(y),$$

bearing in mind that the additional finite element error $u(y) - P_{V_h}^y u(y)$ depends on $h$ (with algebraic decay in the finite element dimension).

All the tests were done using Python 3.8. For more information and experiments not presented here we invite the reader to look into the github repository https://github.com/agussomacal/ROMHighContrast.
## 5.1 Parameter selection

We first study the case of a one parameter family: the diffusion coefficient $y_A$ of $\Omega_A$ in Figure 3 varies from 1 to $\infty$, while the other subdomains are considered as background with all coefficients equal to 1. Thus the $y^k$ are of the form $y^k = (y_A^k, 1, 1, 1)$.

<figure><img src="image_8.png" /><figcaption>Figure 5: Galerkin (left) and $H_0^1$ (right) projection error, both measured in $H_0^1$ relative error, maximized over the parameter domain, for different reduced bases, case $d=1$.</figcaption></figure>

In reduced basis constructions, two approaches for parameter selection are usually considered: random or greedy. Random selection usually performs well enough in many situations, however we shall see that it fails in the high contrast regime. This is in particular due to the fact that it does not capture the limit solutions, while we have observed in §4 that robust convergence of the Galerkin method in the high-contrast regime critically requires to include limit solutions in the space $V_n$. Here, there is only one limit solution $u_\infty = u(y_\infty)$ where $y_\infty = (\infty, 1, 1, 1)$, and this element is picked by the greedy method if initialized at any other point.

More precisely, we compare four strategies for selecting the $y_A^k \in [1, \infty]$:

* Random: the $y_A^k$ are drawn independently according to the uniform law for $\frac{1}{y_A} \in [0, 1]$.

* Random-$\infty$: First the limit solution corresponding to $y_A = \infty$ is put in the basis. The rest of the elements are randomly picked as in the previous case.
* Greedy $H_0^1$: The $y^k$ are picked incrementally, $y^{k+1}$ maximizing the relative $H_0^1$ projection error $\|u(y) - P_{V_k}u(y)\|_{H_0^1}/\|u(y)\|_{H_0^1}$.

* Greedy Galerkin: The $y^k$ are picked incrementally, $y^{k+1}$ maximizing the relative $H_0^1$ error of the Galerkin projection $\|u(y) - P_{V_k}^y u(y)\|_{H_0^1}/\|u(y)\|_{H_0^1}$.

Figure 5 displays on the left the evolution of the maximal relative error of the Galerkin projection

$$ \sup_{y_A \in [1, \infty]} \frac{\|u(y) - P_{V_n}^y u(y)\|_{H_0^1}}{\|u(y)\|_{H_0^1}}, $$

as a function of $n = \dim(V_n)$ for these various selection strategies. It reveals the superiority of the greedy selection that reaches machine precision after picking $n = 11$ reduced basis elements, and the gain in including the limit solution in the case of a random selection. As a comparison, we display on the right the decay of the relative $H_0^1$-orthogonal projection error

$$ \sup_{y_A \in [1, \infty]} \frac{\|u(y) - P_{V_n} u(y)\|_{H_0^1}}{\|u(y)\|_{H_0^1}} $$

for the same parameter selection strategies. Here, we notice that the inclusion of the limit solution $u_\infty$ is not anymore critical for reaching good accuracy. Nevertheless, these errors still decay faster for the greedy strategies.

**Remark 5.1.** As the diffusion coefficient is piecewise constant on the partition $\Omega_A \cup \Omega_A^c$, the parameter space dimension is $d = 2$ in this numerical example. The theoretical results thus provide a bound on the error of order $\exp(-c\sqrt{n})$. However, this bound is obtained with local reduced spaces $V_{\ell,k}$ on dyadic intervals, which does not perform as well as $V_n = \bigoplus_{\ell \in E_k} V_{\ell,k}$, for which one might expect a rate closer to $\exp(-cn)$. In Figure 5 for $n \le 11$, that is, until numerical precision issues arise, we even observe a faster than exponential convergence, that could be due to the superiority of reduced bases over polynomial approximations.

**Remark 5.2.** It is well known that the reduced basis can be very ill-conditioned, since $u^n$ becomes extremely close to $V_{n-1} = \operatorname{span}\{u^1, \dots, u^{n-1}\}$ as $n$ gets moderately large. In order to avoid numerical instabilities, prior to the computation of the Galerkin or $H_0^1$ projection onto $V_n$, we need to perform a change of basis, typically by some orthonormalization process. In our numerical test, we perform this orthonormalization with respect to the discrete $\ell^2$ inner product for the nodal values in the background finite element representation, using the QR decomposition, and obtain a satisfactory stable numerical behavior. However, this process is not invariant under permutations, and we observe that it behaves better in terms of numerical stability when sorting the reduced basis elements from higher contrast to lower contrast.

In this one parameter scenario, both greedy strategies behaved equally well. However, as we increase the dimensionality of the problem $d > 1$, Greedy
<figure><img src="image_9.png" /><figcaption>Figure 6: Galerkin and $H_0^1$ projection error (both measured in $H_0^1$ relative error maximized over the parameter domain) for different reduced bases, case $d = 2$.</figcaption></figure>

Galerkin appears to be the best selection procedure, as could be expected since it optimizes the error based on the approximation which is effectively computed in forward modeling. Figure 6 shows this effect when $d = 2$, where $y_A$ and $y_B$ are allowed to vary independently while $y_C$ and $y_D$ are taken as background always equal to 1.

## 5.2 Influence of dimensionality and geometry

In order to study the impact of dimensionality on the approximation rates, we compare the behavior of the Greedy Galerkin selection method, as we increase the number of freely varying parameters. As before, we will have for $y = (y_A, 1, 1, 1)$ when $d = 1$, then $y = (y_A, y_B, 1, 1)$ when $d = 2$, until having all four subdomains freely varying between 1 and $+\infty$.

In Figure 7 the degradation with respect to dimension is clearly observed as the approximation capabilities strongly decrease. Even thought the exponential decay rate is still conserved, the decay parameter shrinks from almost 3 down to 0.22 when $d = 4$.

Secondly, we study the case where the geometrical assumptions are not satisfied. We follow the same incremental subdomains unfreezing as in the previous case but using the geometry stated in Figure 4. We observe that the reduced basis approach still achieves exponential approximation rates, actually higher than in the previous example. This hints that the geometric assumptions which are needed in our proofs could be artificial, and leaves open the question of achieving such results without relying on these assumptions.
<figure><img src="image_10.png" /><figcaption>Figure 7: The Galerkin projection of Greedy Galerkin method for increasing dimensionality in geometries satisfying (left) or not (right) the assumptions.</figcaption></figure>

References

[1] R. A. Adam and J. F. Fournier, Sobolev spaces, Elsevier, 2003.

[2] B. Aksoyly and Z. Yelter, *Robust multigrid preconditioners for cell-centered finite volume discretization of the high-contrast diffusion equation*, Computing and Visualization in Science 13, 229-245, 2010.

[3] B. Aksoylu, I.G. Graham, H. Klie, and R. Scheichl, *Towards a rigorously justified algebraic preconditioner for high-contrast diffusion problems*, Computing and Visualization 11, 319-331, 2008.

[4] M. Ainsworth, *Robust a posteriori error estimation for nonconforming finite element approximation*, SIAM J. Num. Anal. 42-6, 2320-2341, 2005.

[5] I. Babuška, F. Nobile and R. Tempone, *A stochastic collocation method for elliptic partial differential equations with random input data*, SIAM J. Num. Anal. 45, 1005-1034, 2007.

[6] J. Beck, F. Nobile, L. Tamellini and R. Tempone, *Stochastic spectral Galerkin and collocation methods for PDEs with random coefficients: a numerical comparison*, Lecture Notes in Computational Science and Engineering 76, 43-62, 2010.

[7] J. Beck, F. Nobile, L. Tamellini and R. Tempone, *Implementation of optimal Galerkin and collocation approximations of PDEs with random coefficients*, ESAIM Proc 33, 10-21, 2011.
[8] J. Beck, F. Nobile, L. Tamellini and R. Tempone, *Convergence of quasi-optimal stochastic Galerkin methods for a class of PDES with random coefficients*, Computers & Mathematics with Applications 67(9), 732-751, 2014.

[9] C. Bernardi and R. Verfürth, *Adaptive finite element methods for elliptic equations with non-smooth coefficients*, Numerische Mathematik 85, 579-608, 2000.

[10] M. Bachmayr and A. Cohen, *Kolmogorov widths and low-rank approximations of parametric elliptic PDEs*, Mathematics of Computation, vol. 86, no. 304, 701-724, 2017.

[11] M. Bachmayr, A. Cohen, and G. Migliorati, *Sparse polynomial approximation of parametric elliptic pdes. part i: affine coefficients*, ESAIM:M2AN, 51(1), 321-339, 2017.

[12] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk, *Convergence rates for greedy algorithms in reduced basis methods*, SIAM Journal on Mathematical Analysis, 43(3), 1457-1472, 2011.

[13] P. Binev, A. Cohen, W. Dahmen, R. DeVore, G. Petrova, and P. Wojtaszczyk, *Data assimilation in reduced modeling*, SIAM/ASA Journal on Uncertainty Quantification, 5(1), 1-29, 2017.

[14] A. Bonito, A. Cohen, R. DeVore, D. Guignard, P. Jantsch and G. Petrova, *Nonlinear methods for model reduction*, ESAIM: Mathematical Modelling and Numerical Analysis, 55(2), 507-531, 2021.

[15] A. Buffa, Y. Maday, A.T. Patera, C. Prud'homme, and G. Turinici, *A Priori convergence of the greedy algorithm for the parameterized reduced basis*, Mathematical Modeling and Numerical Analysis, 46, 595-603, 2012.

[16] A. Chatterjee, *An introduction to the proper orthogonal decomposition*, Current Science 78, 539-575, 2000.

[17] A. Cohen and R. DeVore, *Approximation of high-dimensional parametric pdes*, Acta Numerica, 24, 1-159, 2015.

[18] A. Cohen, R. DeVore and C. Schwab, *Analytic regularity and polynomial approximation of parametric and stochastic PDEs*, Analysis and Applications, 9, 11-47, 2011.

[19] R. DeVore, G. Petrova, and P. Wojtaszczyk, *Greedy algorithms for reduced bases in Banach spaces*, Constructive Approximation, 37(3), 455-466, 2013.

[20] J. Galvis and Y. Efendiev, *Domain decomposition preconditioners for multiscale flows in high contrast media: reduced dimension coarse spaces*, SIAM Journal on Multiscale Modeling and Simulation 8(4), 1461-1483, 2010.
[21] B. Haasdonk, *Reduced basis methods for parametrized PDEs - a tutorial introduction for stationary and instationary problems*, in Model Reduction and Approximation - Theory and Algorithms, P. Benner, A. Cohen, M. Ohlberger, and K. Willcox eds, SIAM, 2017.

[22] V.V. Jikov, S.M. Kozlov, and O.A. Oleinik, *Homogeneization of differential operators and integral functionals*, Springer, 1994.

[23] D. Kouri, Z. Zou and W. Aquino, *An adaptive local reduced basis method for solving pdes with uncertain inputs and evaluating risk*, Computer Methods in Applied Mechanics and Engineering 345, 302–322, 2019.

[24] Y. Maday, A.T. Patera, J.D. Penn and M. Yano, *A parametrized-background data-weak approach to variational data assimilation: Formulation, analysis, and application to acoustics*, Int. J. Numer. Meth. Eng. 102, 933-965, 2015. DOI: 10.1002/nme.4747

[25] Y. Maday and B. Stamm, *Locally adaptive greedy approximations for anisotropic parameter reduced basis spaces*, SIAM Journal on Scientific Computing, 35(6), A2417-A2441, 2013.

[26] A.T. Patera, J.L. Eftang and E.M. Ronquist, *An "hp" certified reduced basis method for parametrized elliptic partial differential equations*, SIAM Journal on Scientific Computing, 32(6), 3170-3200, 2010.

[27] A. Pinkus, *N-widths in approximation theory*, Springer, 1985.

[28] G. Rozza, D.B.P. Huynh, and A.T. Patera, *Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations - application to transport and continuum mechanics*, Archive of Computational Methods in Engineering 15, 229-275, 2008. DOI: 10.1007/s00791-006-0044-7

[29] S. Sen, *Reduced-basis approximation and a posteriori error estimation for many-parameter heat conduction problems*, Numerical Heat Transfer B-Fund 54, 369-389, 2008. DOI: 10.1080/10407790802424204

[30] E.M. Stein, *Singular integrals and differentiability properties of functions*, Princeton University Press, 1970.

[31] V. Temlyakov, *Nonlinear Kolmogorov widths*, Math. Notes 63, 785-795, 1998.

[32] H. Tran, C.G. Webster, and G. Zhang, *Analysis of quasi-optimal polynomial approximations for parameterized PDEs with deterministic and stochastic coefficients*, Numerische Mathematik 137, 451-493, 2017.

[33] K. Veroy, C. Prudhomme, D.V. Rovas and T. Patera, *A Posteriori Error Bounds for Reduced-Basis Approximation of Parametrized Noncoercive and Nonlinear Elliptic Partial Differential Equations*, Proc. 16th AIAA Computational Fluid Dynamics Conference, Orlando, 2003.
[34] S. Volkwein, *Proper Orthogonal Decomposition: Theory and reduced order modeling.*, Lecture Notes, University of Konstanz, 2013.

[35] K. Willcox and J. Peraire, *Balanced model reduction via the proper orthogonal decomposition*, American Institute of Aeronautics and Astronautics 40, 2323-2330, 2022.

[36] A. Zenisek, *extensions from the Sobolev spaces $H^1$ satisfying prescribed Dirichlet boundary conditions*, Applications of Mathematics 49, 405-413, 2004.