Nonlinear reduced models for state and parameter estimation

Parametrized partial differential equations are of common used to model complex physical systems. Such equations can generally be written in abstract form as

$${\cal P}(u,y)=0,$$

(1.1)

where $y=(y_{1},\dots,y_{d})$ is a vector of scalar parameters ranging in some domain $Y\subset\mathbb{R}^{d}$. We assume well-posedness, that is, for any $y\in Y$ the problem admits a unique solution $u=u(y)$ in some Hilbert space $V$. We may therefore consider the parameter to solution map

$$y\mapsto u(y),$$

(1.2)

from $Y$ to $V$, which is typically nonlinear, as well as the solution manifold

$${\cal M}:=\{u(y)\,:\,y\in Y\}\subset V$$

(1.3)

that describes the collection of all admissible solutions. Throughout this paper, we assume that $Y$ is compact in $\mathbb{R}^{d}$ and that the map (1.2) is continuous. Therefore ${\cal M}$ is a compact set of $V$. We sometimes refer to the solution $u(y)$ as the state of the system for the given parameter vector $y$.

The parameters are used to represent physical quantities such as diffusivity, viscosity, velocity, source terms, or the geometry of the physical domain in which the PDE is posed. In several relevant instances, $y$ may be high or even countably infinite dimensional, that is, $d\gg 1$ or $d=\infty$.

In this paper, we are interested in inverse problems which occur when only a vector of linear measurements

$$z=(z_{1},\dots,z_{m})\in\mathbb{R}^{m},\quad z_{i}=\ell_{i}(u),\quad i=1,\dots,m,$$

(1.4)

is observed, where each $\ell_{i}\in V^{\prime}$ is a known continuous linear functional on $V$. We also sometimes use the notation

$$z=\ell(u),\quad\ell=(\ell_{1},\dots,\ell_{m}).$$

(1.5)

One wishes to recover from $z$ the unknown state $u\in{\cal M}$ or even the underlying parameter vector $y\in Y$ for which $u=u(y)$. Therefore, in an idealized setting, one partially observes the result of the composition map

$$y\in Y\mapsto u\in{\cal M}\mapsto z\in\mathbb{R}^{m}.$$

(1.6)

for the unknown $y$. More realistically, the measurements may be affected by additive noise

$$z_{i}=\ell_{i}(u)+\eta_{i},$$

(1.7)

and the model itself might be biased, meaning that the true state $u$ deviates from the solution manifold ${\cal M}$ by some amount. Thus, two types of inverse problems may be considered:

1. (i)

   State estimation: recover an approximation $u^{*}$ of the state $u$ from the observation $z=\ell(u)$. This is a linear inverse problem, in which the prior information on $u$ is given by the manifold ${\cal M}$ which has a complex geometry and is not explicitly known.

2. (ii)

   Parameter estimation: recover an approximation $y^{*}$ of the parameter $y$ from the observation $z=\ell(u)$ when $u=u(y)$. This is a nonlinear inverse problem, for which the prior information available on $y$ is given by the domain $Y$.

These problems become severely ill-posed when ${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}Y}$ has dimension $d>m$. For this reason, they are often addressed through Bayesian approaches [17, 25]: a prior probability distribution $P_{y}$ being assumed on $y\in Y$ (thus inducing a push forward distribution $P_{u}$ for $u\in{\cal M}$), the objective is to understand the posterior distributions of $y$ or $u$ conditioned by the observations $z$ in order to compute plausible solutions $y^{*}$ or $u^{*}$ under such probabilistic priors. The accuracy of these solutions should therefore be assessed in some average sense.

In this paper, we do not follow this avenue: the only priors made on $y$ and $u$ are their membership to $Y$ and ${\cal M}$. We are interested in developping practical estimation methods that offer uniform recovery guarantees under such deterministic priors in the form of upper bounds on the worst case error for the estimators over all $y\in Y$ or $u\in{\cal M}$. We also aim to understand whether our error bounds are optimal in some sense. Our primary focus will actually be on state estimation (i). Nevertheless we present in §4 several implications on parameter estimation (ii), which to our knowledge are new. For state estimation, error bounds have recently been established for a class of methods based on linear reduced modeling, as we recall next.

In several relevant instances, the particular parametrized PDE structure allows one to construct linear spaces $V_{n}$ of moderate dimension $n$ that are specifically tailored to the approximation of the solution manifold ${\cal M}$, in the sense that

$${\rm dist}({\cal M},V_{n})=\max_{u\in{\cal M}}\min_{v\in V_{n}}\|u-v\|\leq\varepsilon_{n},$$

(1.8)

where $\varepsilon_{n}$ is a certified bound that decays with $n$ significanly faster than when using for $V_{n}$ classical approximation spaces of dimension $n$ such as finite elements, algebraic or trigonometric polynomials, or spline functions. Throughout this paper

$${\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\langle\cdot,\cdot\rangle_{V}^{1/2}=:}\|\cdot\|=\|\cdot\|_{V},$$

(1.9)

denotes the norm of the Hilbert space $V$. The natural benchmark for such approximation spaces is the Kolmogorov $n$-width

$$d_{n}({\cal M}):=\min_{\dim(E)=n}{\rm dist}({\cal M},E).$$

(1.10)

The space $E_{n}$ that achieves the above minimum is thus the best possible reduced model for approximating all of ${\cal M}$, however it is computationally out of reach.

One instance of computational reduced model spaces is generated by sparse polynomial approximations of the form

$$u_{n}(y)=\sum_{\nu\in\Lambda_{n}}u_{\nu}y^{\nu},\quad y^{\nu}:=\prod_{j\geq 1}y_{j}^{\nu_{j}},$$

(1.11)

where $\Lambda_{n}$ is a conveniently chosen set of multi-indices such that $\#(\Lambda_{n})=n$. Such approximations can be derived, for example, by best $n$-term truncations of infinite Taylor or orthogonal polynomial expansions. We refer to [11, 12] where convergence estimates of the form

$$\sup_{y\in Y}\|u(y)-u_{n}(y)\|\leq Cn^{-s},$$

(1.12)

are established for some $s>0$ even when $d=\infty$. Therefore, the space $V_{n}:={\rm span}\{u_{\nu}\,:\,\nu\in\Lambda_{n}\}$ approximates the solution manifold with accuracy $\varepsilon_{n}=Cn^{-s}$.

Another instance, known as reduced basis approximation, consists of using spaces of the form

$$V_{n}:={\rm span}\{u_{1},\dots,u_{n}\},$$

(1.13)

where $u_{i}=u(y^{i})\in{\cal M}$ are instances of solutions corresponding to a particular selection of parameter values $y^{i}\in Y$ (see [20, 23, 24]). One typical selection procedure is based on a greedy algorithm: one picks $y^{k}$ such that $u_{k}=u(y^{k})$ is furthest away from the previously constructed space $V_{k-1}$ in the sense of maximizing a computable and tight a-posteriori bound of the projection error $\|u(y)-P_{V_{k-1}}u(y)\|$ over a sufficiently fine discrete training set $\tilde{Y}\subset Y$. In turn, this method also delivers a computable upper estimate $\varepsilon_{k}$ for ${\rm dist}({\cal M},V_{k})$. It was proved in [1, 16] that the reduced basis spaces resulting from this greedy algorithm have near-optimal approximation property, in the sense that if $d_{n}({\cal M})$ has a certain polynomial or exponential rate of decay as $n\to\infty$, then the same rate is achieved by ${\rm dist}({\cal M},V_{n})$.

In both cases, these reduced models come in the form of a hierarchy $(V_{n})_{n\geq 1}$, with computable decreasing error bounds $(\varepsilon_{n})_{n\geq 1}$, where $n$ corresponds to the level of truncation in the first case and the step of the greedy algorithm in the second case. Given a reduced model $V_{n}$, one way of tackling the state estimation problem is to replace the complex solution manifold ${\cal M}$ by the simpler prior class described by the cylinder

$${\cal K}={\cal K}(V_{n},\varepsilon_{n})=\{v\in V\;:\;\mathop{\rm dist}(v,V_{n})\leq\varepsilon_{n}\}.$$

(1.14)

that contains ${\cal M}$. The set ${\cal K}$ therefore reflects the approximability of ${\cal M}$ by $V_{n}$. This point of view leads to the Parametrized Background Data Weak (PBDW) method introduced in [19], also termed as one space method and further analyzed in [2], that we recall below in a nutshell.

In the noiseless case, the knowledge of $z=(z_{i})_{i=1,\dots,m}$ is equivalent to that of the orthogonal projection $w=P_{W}u$, where

$$W:={\rm span}\{\omega_{1},\dots,\omega_{m}\}$$

(1.15)

and $\omega_{i}\in V$ are the Riesz representers of the linear functionals $\ell_{i}$, that is

$$\ell_{i}(v)=\langle\omega_{i},v\rangle,\quad v\in V.$$

(1.16)

Thus, the data indicates that $u$ belongs to the affine space

$$V_{w}:=w+W^{\perp}.$$

(1.17)

Combining this information with the prior class ${\cal K}$, the unknown state thus belongs to the ellipsoid

$${\cal K}_{w}:={\cal K}\cap V_{w}=\{v\in{\cal K}\;:\;P_{W}v=w\}.$$

(1.18)

For this posterior class ${\cal K}_{w}$, the optimal recovery estimator $u^{*}$ that minimizes the worst case error $\max_{u\in{\cal K}_{w}}\|u-u^{*}\|$ is therefore the center of the ellipsoid, which is equivalently given by

$$u^{*}=u^{*}(w):={\rm argmin}\{\|v-P_{V_{n}}v\|\;:\;P_{W}v=w\}.$$

(1.19)

It can be computed from the data $w$ in an elementary manner by solving a finite set of linear equations. The worst case performance for this estimator, both over ${\cal K}$ and ${\cal K}_{w}$, for any $w$, is thus given by the half-diameter of the ellipsoid which is the product of the width $\varepsilon_{n}$ of ${\cal K}$ and the quantity

$$\mu_{n}=\mu(V_{n},W):=\max_{v\in V_{n}}\frac{\|v\|}{\|P_{W}v\|}.$$

(1.20)

Note that $\mu_{n}$ is the inverse of the cosine of the angle between $V_{n}$ and $W$. For $n\geq 1$, this quantity can be computed as the inverse of the smallest singular value of the $n\times m$ cross-Gramian matrix with entries $\langle\phi_{i},\psi_{j}\rangle$ between any pair of orthonormal bases $(\phi_{i})_{i=1,\dots,n}$ and $(\psi_{j})_{j=1,\dots,m}$ of $V_{n}$ and $W$, respectively. It is readily seen that one also has

$$\mu_{n}=\max_{w\in W^{\perp}}\frac{\|w\|}{\|P_{V_{n}^{\perp}}w\|},$$

(1.21)

allowing us to extend the above definition to the case of the zero-dimensional space $V_{n}=\{0\}$ for which $\mu(\{0\},W)=1$.

Since ${\cal M}\subset{\cal K}$, the worst case error bound over ${\cal M}$ of the estimator, defined as

$$E_{wc}:=\max_{u\in{\cal M}}\|u-u^{*}(P_{W}u)\|,$$

(1.22)

satisfies the error bound

$$E_{wc}=\max_{u\in{\cal M}}\|u-u^{*}(P_{W}u)\|\leq\max_{u\in{\cal K}}\|u-u^{*}(P_{W}u)\|{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}=\mu_{n}\varepsilon_{n}}.$$

(1.23)

Note that $\mu_{n}\geq 1$ increases with $n$ and that its finiteness imposes that $\dim(V_{n})\leq\dim(W)$, that is $m\geq n$. Therefore, one natural way to decide which space $V_{n}$ to use is to take the value of $n\in\{0,\dots,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{m}}\}$ that minimizes the bound $\mu_{n}\varepsilon_{n}$. This choice is somehow crude since it might not be the value of $n$ that minimizes the true reconstruction error for a given $u\in{\cal M}$, and for this reason it was referred to as a poor man algorithm in [1].

The PBDW approach to state estimation can be improved in various ways:

• •

  One variant that is relevant to the present work is to use reduced models of affine form

  $$V_{n}=\overline{u}_{n}+\overline{V}_{n},$$

  (1.32)

  where $\overline{V}_{n}$ is a linear space and $\overline{u}$ is a given offset. The optimal recovery estimator is again defined by the minimization property (1.19). Its computation amounts to the same type of linear systems and the reconstruction map $w\mapsto u^{*}(w)$ is now affine. The error bound (1.23) remains valid with $\mu_{n}=\mu(\overline{V}_{n},W)$ and $\varepsilon_{n}$ a bound for ${\rm dist}({\cal M},V_{n})$. Note that $\varepsilon_{n}$ is also a bound for the distance of ${\cal M}$ to the linear space $\overline{V}_{n+1}:=\overline{V}_{n}\,\oplus\,\mathbb{R}\overline{u}_{n}$ of dimension $n+1$. However, using instead this linear space, could result in a stability constant $\mu_{n+1}=\mu(\overline{V}_{n+1},W)$ that is much larger than $\mu_{n}$, in particular, when the offset $\overline{u}_{n}$ is close to $W^{\perp}$.

• •

  Another variant proposed in [10] consists in using a large set ${\cal T}_{N}=\{u_{i}=u(y^{i})\,:\,i=1,\dots,N\}$ of precomputed solutions in order to train the reconstruction maps $w\mapsto u^{*}(w)$ by minimizing the least-square fit $\sum_{i=1}^{N}\|u_{i}-u^{*}(P_{W}u_{i})\|^{2}$ over all linear or affine map, which amounts to optimizing the choice of the space $V_{n}$ in the PBDW method.

• •

  Conversely, for a given reduced basis space $V_{n}$, it is also possible to optimize the choice of linear functionals $(\ell_{1},\dots,\ell_{m})$ giving rise to the data, among a dictionary ${\cal D}$ that represent a set of admissible measurement devices. The objective is to minimize the stability constant $\mu(V_{n},W)$ for the resulting space $W$, see in particular [3] where a greedy algorithm is proposed for selecting the $\ell_{i}$. We do not take this view in the present paper and think of the space $W$ as fixed once and for all: the measurement devices are given to us and cannot be modified.

The simplicity of the PBDW method and its above variants come together with a fundamental limitation of its performance: since the map $w\mapsto u^{*}(w)$ is linear or affine, the reconstruction necessarily belongs to an $m$ or $m+1$ dimensional space, and thefore the worst case performance is necessarily bounded from below by the Kolmogorov width $d_{m}({\cal M})$ or $d_{m+1}({\cal M})$. In view of this limitation, one principle objective of the present work is to develop nonlinear state estimation techniques which provably overcome the bottleneck of the Kolmogorov width $d_{m}({\cal M})$.

In §2, we introduce various benchmark quantities that describe the best possible performance of a recovery map in a worst case sense. We first consider an idealized setting where the state $u$ is assumed to exactly satisfy the theoretical model described by the parametric PDE, that is $u\in{\cal M}$. Then we introduce similar benchmarks quantities in the presence of model bias and measurement noise. All these quantities can be substantially smaller than $d_{m}({\cal M})$.

In §3, we discuss a nonlinear recovery method, based on a family of affine reduced models $(V_{k})_{k=1,\dots,K}$, where each $V_{k}$ has dimension $n_{k}\leq m$ and serves as a local approximations to a portion ${\cal M}_{k}$ of the solution manifold. Applying the PBDW method with each such space, results in a collection of state estimators $u^{*}_{k}$. The value $k$ for which the true state $u$ belongs to ${\cal M}_{k}$ being unknown, we introduce a model selection procedure in order to pick a value $k^{*}$, and define the resulting estimator $u^{*}=u^{*}_{k^{*}}$. We show that this estimator has performance comparable to the benchmark introduced in §2. Such performances cannot be achieved by the standard PBDW method due to the above described limitations.

Model selection is a classical topic of mathematical statistics [22], with representative techniques such as complexity penalization or cross-validation in which the data are used to select a proper model. Our approach differs from these techniques in that it exploits (in the spirit of data assimilation) the PDE model which is available to us, by evaluating the distance to the manifold

$$\mathop{\rm dist}(v,{\cal M})=\min_{y\in Y}\|v-u(y)\|,$$

(1.33)

of the different estimators $v=u^{*}_{k}$ for $k=1,\dots,K$, and picking the value $k^{*}$ that minimizes it. In practice, the quantity (1.33) cannot be exactly computed and we instead rely on a computable surrogate quantity ${\cal S}(v,{\cal M})$ expressed in terms of the residual to the PDE, see § 3.4.

One typical instance where such a surrogate is available is when (1.1) has the form of a linear operator equation

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

(1.34)

where $A(y)$ is boundedly invertible from $V$ to $V^{\prime}$, or more generally, from $V\to Z^{\prime}$ for a test space $Z$ different from $V$, uniformly over $y\in Y$. Then ${\cal S}(v,{\cal M})$ is obtained by minimizing the residual

$${\cal R}(v,y)=\|A(y)v{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}f(y)\|_{Z^{\prime}},$$

(1.35)

over $y\in Y$. This task itself is greatly facilitated in the case where the operators $A(y)$ and source terms $f(y)$ have affine dependence in $y$. One relevant example that has been often considered in the literature is the second order elliptic diffusion equation with affine diffusion coefficient,

$$-{\rm div}(a\nabla u)=f(y),\quad a=a(y)=\overline{a}+\sum_{j=1}^{d}y_{j}\psi_{j}.$$

(1.36)

In §4, we discuss the more direct approach for both state and parameter estimation based on minimizing ${\cal R}(v,y)$ over both $y\in Y$ and $v\in w+W^{\perp}$. The associated alternating minimization algorithm amounts to a simple succession of quadratic problems in the particular case of linear PDE’s with affine parameter dependence. Such an algorithm is not guaranteed to converge to a global minimum (since the residual is not globally convex), and for this reason its limit may miss the optimaliy benchmark. On the other hand, using the estimator derived in §3 as a “good initialization point” to this mimimization algorithm leads to a limit state that has at least the same order of accuracy.

These various approaches are numerically tested in §5 for the elliptic equation (1.36), for both the overdetermined regime $m\geq d$, and the underdetermined regime $m<d$.

In the absence of model bias and when a noiseless measurement $w=P_{W}u$ is given, our knowledge on $u$ is that it belongs to the set

$${\cal M}_{w}:={\cal M}\cap V_{w}.$$

(2.2)

The best possible recovery map can be described through the following general notion.

In particular one has

$$\frac{1}{2}{\rm diam}(S)\leq{\rm rad}(S)\leq{\rm diam}(S),$$

(2.3)

where ${\rm diam}(S):={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sup}\{\|u-v\|\,:\,u,v\in S\}$ is the diameter of $S$. Therefore, the recovery map that minimizes the worst case error over ${\cal M}_{w}$ for any given $w$, and therefore over ${\cal M}$ is defined by

$$u^{*}(w)={\rm cen}({\cal M}_{w}).$$

(2.4)

Its worst case error is

$$E_{\rm wc}^{*}=\sup\{\mathop{\rm rad}({\cal M}_{w})\,:\,w\in W\}.$$

(2.5)

In view of the equivalence (2.3), we can relate $E^{*}_{\rm wc}$ to the quantity

$$\delta_{0}=\delta_{0}({\cal M},W):={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\sup}\{{\rm diam}({\cal M}_{w})\,:\,w\in W\}=\sup\{\|u-v\|\;:\;u,v\in{\cal M},\;u-v\in W^{\bot}\},$$

(2.6)

by the equivalence

$$\frac{1}{2}\delta_{0}\leq E^{*}_{\rm wc}\leq\delta_{0}.$$

(2.7)

Note that injectivity of the measurement map $P_{W}$ over ${\cal M}$ is equivalent to $\delta_{0}=0$. We provide in Figure (1(a)) an illustration the above benchmark concepts.

If $w=P_{W}u$ for some $u\in{\cal M}$, then any $u^{*}\in{\cal M}$ such that $P_{W}u^{*}=w$, meets the ideal benchmark $\|u-u^{*}\|\leq\delta_{0}$. Therefore, one way for finding such a $u^{*}$ would be to minimize the distance to the manifold over all functions such that $P_{W}v=w$, that is, solve

$$\min_{v\in V_{w}}{\rm dist}(v,{\cal M})=\min_{v\in V_{w}}\min_{y\in Y}\|u(y)-v\|.$$

(2.8)

This problem is computationally out of reach since it amounts to the nested minimization of two non-convex functions in high dimension.

Computationally feasible algorithms such as the PBDW methods are based on a simplification of the manifold ${\cal M}$ which induces an approximation error. We introduce next a somewhat relaxed benchmark that takes this error into account.

In order to account for manifold simplification as well as model bias, for any given accucary $\sigma>0$, we introduce the $\sigma$-offset of ${\cal M}$,

$${\cal M}_{\sigma}:=\{v\in V\;:\;\mathop{\rm dist}(v,{\cal M})\leq\sigma\}=\bigcup_{u\in{\cal M}}B(u,\sigma).$$

(2.9)

Likewise, we introduce the perturbed set

$${\cal M}_{\sigma,w}={\cal M}_{\sigma}\cap V_{w},$$

(2.10)

which, however, still excludes uncertainties in $w$. Our benchmark for the worst case error is now defined as (see Figure (1(b)) for an illustration)

$$\delta_{\sigma}:=\max_{w\in W}{\rm diam}({\cal M}_{\sigma,w})=\max\{\|u-v\|\;:\;u,v\in{\cal M}_{\sigma},\;u-v\in W^{\bot}\}.$$

(2.11)

The map $\sigma\mapsto\delta_{\sigma}$ satisfies some elementary properties:

• •

  Monotonicity and continuity: it is obviously non-decreasing

  $$\sigma\leq\tilde{\sigma}\implies\delta_{\sigma}\leq\delta_{\tilde{\sigma}}.$$

  (2.12)

  Simple finite dimensional examples show that this map may have jump discontinuities. Take for example a compact set ${\cal M}\subset\mathbb{R}^{2}$ consisting of the two points $(0,0)$ and $(1/2,1)$, and $W=\mathbb{R}e_{1}$ where $e_{1}=(1,0)$. Then $\delta_{\sigma}=2\sigma$ for $0\leq\sigma\leq\frac{1}{4}$, while $\delta_{\frac{1}{4}}({\cal M},W)=1$. Using the compactness of ${\cal M}$, it is possible to check that $\sigma\mapsto\delta_{\sigma}$ is continuous from the right and in particular $\lim_{\sigma\to 0}\delta_{\sigma}({\cal M},W)=\delta_{0}$.

• •

  Bounds from below and above: for any $u,v\in{\cal M}_{\sigma,w}$, and for any $\tilde{\sigma}\geq 0$, let $\tilde{u}=u+\tilde{\sigma}g$ and $\tilde{v}=v-\tilde{\sigma}g$ with $g=(u-v)/\|u-v\|$. Then, $\|\tilde{u}-\tilde{v}\|=\|u-v\|+2{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\tilde{\sigma}}}$ and $\tilde{u}-\tilde{v}\in W^{\bot}$, which shows that $\tilde{u},\,\tilde{v}\in{\cal M}_{\sigma+\tilde{\sigma},w}$, and

  $$\delta_{\sigma+\tilde{\sigma}}\geq\delta_{\sigma}+2\tilde{\sigma}.$$

  (2.13)

  In particular,

  $$\delta_{\sigma}\geq\delta_{0}+2\sigma\geq 2\sigma.$$

  (2.14)

  On the other hand, we obviously have the upper bound $\delta_{\sigma}\leq{\rm diam}({\cal M}_{\sigma})\leq{\rm diam}({\cal M})+2\sigma$.

• •

  The quantity

  $$\mu({\cal M},W):=\frac{1}{2}\sup_{\sigma>0}\frac{\delta_{\sigma}-\delta_{0}}{\sigma},$$

  (2.15)

  may be viewed as a general stability constant inherent to the recovery problem, similar to $\mu(V_{n},W)$ that is more specific to the particular PBDW method: in the special case where ${\cal M}=V_{n}$ and $V_{n}\cap W^{\perp}=\{0\}$, one has $\delta_{0}=0$ and $\frac{\delta_{\sigma}}{\sigma}=\mu(V_{n},W)$ for all $\sigma>0$. Note that $\mu({\cal M},W)\geq 1$ in view of (2.14).

Regarding measurement noise, it suggests to introduce the quantity

$$\tilde{\delta}_{\sigma}:=\max\{\|u-v\|\;:\;u,v\in{\cal M},\;\|P_{W}u-P_{W}v\|\leq\sigma\}.$$

(2.16)

The two quantities $\delta_{\sigma}$ and $\tilde{\delta}_{\sigma}$ are not equivalent, however one has the framing

$$\delta_{\sigma}-2\sigma\leq\tilde{\delta}_{2\sigma}\leq\delta_{\sigma}+2\sigma.$$

(2.17)

To prove the upper inequality, we note that for any $u,v\in{\cal M}_{\sigma}$ such that $u-v\in W^{\bot}$, there exists $\tilde{u},\tilde{v}\in{\cal M}$ at distance $\sigma$ from $u,v$, respectively, and therefore such that $\|P_{W}(\tilde{u}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}\tilde{v})\|\leq 2\sigma$. Conversely, for any $\tilde{u},\tilde{v}\in{\cal M}$ such that $\|P_{W}(\tilde{u}{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}\tilde{v})\|\leq 2\sigma$, there exists $u,v$ at distance $\sigma$ from $\tilde{u},\tilde{v}$, respectively, such that $u-v\in W^{\bot}$, which gives the lower inequality. Note that the upper inequality in (2.17) combined with (2.14) implies that $\tilde{\delta}_{2\sigma}\leq 2\delta_{\sigma}$.

In the following analysis of reconstruction methods, we use the quantity $\delta_{\sigma}$ as a benchmark which, in view of this last observation, also accounts for the lack of accuracy in the measurement of $P_{W}u$. Our objective is therefore to design an algorithm that, for a given tolerance $\sigma>0$, recovers from the measurement $w=P_{W}u$ an approximation to $u$ with accuracy comparable to $\delta_{\sigma}$. Such an algorithm requires that we are able to capture the solution manifold up to some tolerance $\varepsilon\leq\sigma$ by some reduced model.

Linear or affine reduced models, as used in the PBDW algorithm, are not suitable for approximating the solution manifold when the required tolerance $\varepsilon$ is too small. In particular, when $\varepsilon{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}<}d_{m}({\cal M})$ one would then need to use a linear space $V_{n}$ of dimension $n{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}>}m$, therefore making $\mu(V_{n},W)$ infinite.

One way out is to replace the single space $V_{n}$ by a family of affine spaces

$$V_{k}=\overline{u}_{k}+\overline{V}_{k},\quad k=1,\dots,K,$$

(3.1)

each of them having dimension

$$\dim(V_{k})=n_{k}\leq m,$$

(3.2)

such that the manifold is well captured by the union of these spaces, in the sense that

$${\rm dist}\Bigl{(}{\cal M},\bigcup_{k=1}^{K}V_{k}\Bigr{)}\leq\varepsilon$$

(3.3)

for some prescribed tolerance $\varepsilon>0$. This is equivalent to saying that there exists a partition of the solution manifold

$${\cal M}=\bigcup_{k=1}^{K}{\cal M}_{k},$$

(3.4)

such that we have local certified bounds

$${\rm dist}({\cal M}_{k},V_{k})\leq\varepsilon_{k}\leq\varepsilon,\quad k=1,\dots,K.$$

(3.5)

We may thus think of the family $(V_{k})_{k=1,\dots,K}$ as a piecewise affine approximation to ${\cal M}$. We stress that, in contrast to the hierarchies $(V_{n})_{n=0,\dots,m}$ of reduced models discussed in §1.2, the spaces $V_{k}$ do not have dimension $k$ and are not nested. Most importantly, $K$ is not limited by $m$ while each $n_{k}$ is.

The objective of using a piecewise reduced model in the context of state estimation is to have a joint control on the local accuracy $\varepsilon_{k}$ as expressed by (3.5) and on the stability of the PBDW when using any individual $V_{k}$. This means that, for some prescribed $\mu>1$, we ask that

$$\mu_{k}=\mu(\overline{V}_{k},W)\leq\mu,\quad k=1,\dots,K.$$

(3.6)

According to (1.23), the worst case error bound over ${\cal M}_{k}$ when using the PBDW method with a space $V_{k}$ is given by the product $\mu_{k}\varepsilon_{k}$. This suggests to alternatively require from the collection $(V_{k})_{k=1,\dots,K}$, that for some prescribed $\sigma>0$, one has

$$\sigma_{k}:=\mu_{k}\varepsilon_{k}\leq\sigma,{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}{\quad k=1,\dots,K.}}$$

(3.7)

This leads us to the following definitions.

Obviously, any $(\varepsilon,\mu)$-admissible family is $\sigma$-admissible with $\sigma:=\mu\varepsilon$. In this sense the notion of $(\varepsilon,\mu)$-admissibility is thus more restrictive than that of $\sigma$-admissibility. The benefit of the first notion is in the uniform control on the size of $\mu$ which is critical in the presence of noise, as hinted at by Remark 1.1.

If $u\in{\cal M}$ is our unknown state and $w=P_{W}u$ is its observation, we may apply the PBDW method for the different $V_{k}$ in the given family, which yields a corresponding family of estimators

$$u_{k}^{*}=u_{k}^{*}(w)={\rm argmin}\{\mathop{\rm dist}(v,V_{k})\,:\,v\in V_{w}\},\quad k=1,\dots,K.$$

(3.8)

If $(V_{k})_{k=1,\dots,K}$ is $\sigma$-admissible, we find that the accuracy bound

$$\|u-u_{k}^{*}\|\leq\mu_{k}{\rm dist}(u,V_{k})\leq\mu_{k}\varepsilon_{k}=\sigma_{k}\leq\sigma,$$

(3.9)

holds whenever $u\in{\cal M}_{k}$.

Therefore, if in addition to the observed data $w$ one had an oracle giving the information on which portion ${\cal M}_{k}$ of the manifold the unknown state sits, we could derive an estimator with worst case error

$$E_{wc}\leq\sigma.$$

(3.10)

This information is, however, not available and such a worst case error estimate cannot be hoped for, even with an additional multiplicative constant. Indeed, as we shall see below, $\sigma$ can be fixed arbitrarily small by the user when building the family $(V_{k})_{k=1,\dots,K}$, while we know from §2.1 that the worst case error is bounded from below by $E_{wc}^{*}\geq\frac{1}{2}\delta_{0}$ which could be non-zero. We will thus need to replace the ideal choice of $k$ by a model selection procedure only based on the data $w$, that is, a map

$$w\mapsto k^{*}(w),$$

(3.11)

leading to a choice of estimator $u^{*}=u^{*}_{k^{*}}$. We shall prove further that such an estimator is able to achieve the accuracy

$$E_{wc}\leq\delta_{\sigma},$$

(3.12)

that is, the benchmark introduced in §2.2. Before discussing this model selection, we discuss the existence and construction of $\sigma$-admissible or $(\varepsilon,\mu)$-admissible families.

For any arbitrary choice of $\varepsilon>0$ and $\mu\geq 1$, the existence of an $(\varepsilon,\mu)$-admissible family results from the following observation: since the manifold ${\cal M}$ is a compact set of $V$, there exists a finite $\varepsilon$-cover of ${\cal M}$, that is, a family $\overline{u}_{1},\dots,\overline{u}_{K}\in V$ such that

$${\cal M}\subset\bigcup_{k=1}^{K}B(\overline{u}_{k},\varepsilon),$$

(3.13)

or equivalently, for all $v\in{\cal M}$, there exists a $k$ such that $\|v-u_{k}\|\leq\varepsilon$. With such an $\varepsilon$ cover, we consider the family of trivial affine spaces defined by

$$V_{k}=\{\overline{u}_{k}\}=\overline{u}_{k}+\overline{V}_{k},\quad\overline{V}_{k}=\{0\},$$

(3.14)

thus with $n_{k}=0$ for all $k$. The covering property implies that (3.5) holds. On the other hand, for the $0$ dimensional space, one has

$$\mu(\{0\},W)=1,$$

(3.15)

and therefore (3.6) also holds. The family $(V_{k})_{k=1,\dots,K}$ is therefore $(\varepsilon,\mu)$-admissible, and also $\sigma$-admissible with $\sigma=\varepsilon$.

This family is however not satisfactory for algorithmic purposes for two main reasons. First, the manifold is not explicely given to us and the construction of the centers $\overline{u}_{k}$ is by no means trivial. Second, asking for an $\varepsilon$-cover, would typically require that $K$ becomes extremely large as $\varepsilon$ goes to $0$. For example, assuming that the parameter to solution $y\mapsto u(y)$ has Lipschitz constant $L$,

$$\|u(y)-u(\tilde{y})\|\leq L|y-\tilde{y}|,\quad y,\tilde{y}\in Y,$$

(3.16)

for some norm $|\cdot|$ of $\mathbb{R}^{d}$, then an $\varepsilon$ cover for ${\cal M}$ would be induced by an $L^{-1}\varepsilon$ cover for $Y$ which has cardinality $K$ growing like $\varepsilon^{-d}$ as $\varepsilon\to 0$. Having a family of moderate size $K$ is important for the estimation procedure since we intend to apply the PBDW method for all $k=1,\dots,K$.

In order to construct $(\varepsilon,\mu)$-admissible or $\sigma$-admissible families of better controlled size, we need to split the manifold in a more economical manner than through an $\varepsilon$-cover, and use spaces $V_{k}$ of general dimensions $n_{k}\in\{0,\dots,m\}$ for the various manifold portions ${\cal M}_{k}$. To this end, we combine standard constructions of linear reduced model spaces with an iterative splitting procedure operating on the parameter domain $Y$. Let us mention that various ways of splitting the parameter domain have already been considered in order to produce local reduced bases having both, controlled cardinality and prescribed accuracy [18, 21, 5]. Here our goal is slightly different since we want to control both the accuracy $\varepsilon$ and the stability $\mu$ with respect to the measurement space $W$.

We describe the greedy algorithm for constructing $\sigma$-admissible families, and explain how it should be modified for $(\varepsilon,\mu)$-admissible families. For simplicity we consider the case where $Y$ is a rectangular domain with sides parallel to the main axes, the extension to a more general bounded domain $Y$ being done by embedding it in such a hyper-rectangle. We are given a prescribed target value $\sigma>0$ and the splitting procedure starts from $Y$.

At step $j$, a disjoint partition of $Y$ into rectangles $(Y_{k})_{k=1,\dots,K_{j}}$ with sides parallel to the main axes has been generated. It induces a partition of ${\cal M}$ given by

$${\cal M}_{k}:=\{u(y)\,:\,y\in Y_{k}\},\quad k=1,\dots,K_{j}.$$

(3.17)

To each $k\in\{1,\dots,K_{j}\}$ we associate a hierarchy of affine reduced basis spaces

$$V_{n,k}=\overline{u}_{k}+\overline{V}_{n,k},\quad n=0,\dots,m.$$

(3.18)

where $\overline{u}_{k}=u(\overline{y}_{k})$ with $\overline{y}_{k}$ the vector defined as the center of the rectangle $Y_{k}$. The nested linear spaces

$$\overline{V}_{0,k}\subset\overline{V}_{1,k}\subset\cdots\subset\overline{V}_{m,k},\quad\dim({\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}\overline{V}}_{n,k})=n,$$

(3.19)

are meant to approximate the translated portion of the manifold ${\cal M}_{k}-\overline{u}_{k}$. For example, they could be reduced basis spaces obtained by applying the greedy algorithm to ${\cal M}_{k}-\overline{u}_{k}$, or spaces resulting from local $n$-term polynomial approximations of $u(y)$ on the rectangle $Y_{k}$. Each space $V_{n,k}$ has a given accuracy bound and stability constant

$${\rm dist}({\cal M}_{k},V_{n,k})\leq\varepsilon_{n,k}\quad{\rm and}\quad\mu_{n,k}:=\mu(\overline{V}_{n,k},W).$$

(3.20)

We define the test quantity

$$\tau_{k}=\min_{n=0,\dots,m}\mu_{n,k}\varepsilon_{n,k}.$$

(3.21)

If $\tau_{k}\leq\sigma$, the rectangle $Y_{k}$ is not split and becomes a member of the final partition. The affine space associated to ${\cal M}_{k}$ is

$$V_{k}=\overline{u}_{k}+\overline{V}_{k},$$

(3.22)

where $V_{k}=V_{n,k}$ for the value of $n$ that minimizes $\mu_{n,k}\varepsilon_{n,k}$. The rectangles $Y_{k}$ with $\tau_{k}>\sigma$ are, on the other hand, split into a finite number of sub-rectangles in a way that we discuss below. This results in the new larger partition $(Y_{k})_{k=1,\dots,K_{j+1}}$ after relabelling the $Y_{k}$. The algorithm terminates at the step $j$ as soon as $\tau_{k}\leq\sigma$ for all $k=1,\dots,K_{j}=K$, and the family $(V_{k})_{k=1,\dots,K}$ is $\sigma$-admissible. In order to obtain an $(\varepsilon,\mu)$-admissible family, we simply modify the test quantity $\tau_{k}$ by defining it instead as

$$\tau_{k}:=\min_{n=0,\dots,m}\max\Big{\{}\frac{\mu_{n,k}}{\mu},\frac{\varepsilon_{n,k}}{\varepsilon}\Big{\}}$$

(3.23)

and splitting the cells for which $\tau_{k}>1$.

The splitting of one single rectangle $Y_{k}$ can be performed in various ways. When the parameter dimension $d$ is moderate, we may subdivide each side-length at the mid-point, resulting into $2^{d}$ sub-rectangles of equal size. This splitting becomes too costly as $d$ gets large, in which case it is preferable to make a choice of $i\in\{1,\dots,d\}$ and subdivide $Y_{k}$ at the mid-point of the side-length in the $i$-coordinate, resulting in only $2$ sub-rectangles. In order to decide which coordinate to pick, we consider the $d$ possibilities and take the value of $i$ that minimizes the quantity

$$\tau_{k,i}=\max\{\tau_{k,i}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}},\tau_{k,i}^{+}\},$$

(3.24)

where $(\tau_{k,i}^{{\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}-}},\tau_{k,i}^{+})$ are the values of $\tau_{k}$ for the two subrectangles obtained by splitting along the $i$-coordinate. In other words, we split in the direction that decreases $\tau_{k}$ most effectively. In order to be certain that all sidelength are eventually split, we can mitigate the greedy choice of $i$ in the following way: if $Y_{k}$ has been generated by $l$ consecutive refinements, and therefore has volume $|Y_{k}|=2^{-l}|Y|$, and if $l$ is even, we choose $i={\color[rgb]{0,0,0}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,0}\pgfsys@color@gray@stroke{0}\pgfsys@color@gray@fill{0}(l/2\,{\rm mod}\,d)}$. This means that at each even level we split in a cyclic manner in the coordinates $i\in\{1,\dots,d\}$.

Using such elementary splitting rules, we are ensured that the algorithm must terminate. Indeed, we are guaranteed that for any $\eta>0$, there exists a level $l=l(\eta)$ such that any rectangle $Y_{k}$ generated by $l$ consecutive refinements has side-length smaller than $2\eta$ in each direction. Since the parameter-to-solution map is continuous, for any $\varepsilon>0$, we can pick $\eta>0$ such that

$$\|y-\tilde{y}\|_{\ell^{\infty}}\leq\eta\implies\|u(y)-u(\tilde{y})\|\leq\varepsilon,\quad y,\tilde{y}\in Y.$$

(3.25)

Applying this to $y\in Y_{k}$ and $\tilde{y}=\overline{y}_{k}$, we find that for $\overline{u}_{k}=u(\overline{y}_{k})$

$$\|u-\overline{u}_{k}\|\leq\varepsilon,\quad u\in{\cal M}_{k}.$$

(3.26)

Therefore, for any rectangle $Y_{k}$ of generation $l$, we find that the trivial affine space $V_{k}=\overline{u}_{k}$ has local accuracy $\varepsilon_{k}\leq\varepsilon$ and $\mu_{k}=\mu(\{0\},W)=1\leq\mu$, which implies that such a rectangle would not anymore be refined by the algorithm.