Adaptive PBDW approach to state estimation: noisy observations; user-defined update spaces
Yvon Maday Tommaso Taddei
math.NA 62-07, 93E24



[5]







patterns

arrows

 




 
     
 






thefnmark⋆footnotetext



Authors
  

  
Addresses
  

  
  
  
   
  




plain
	
	[C]submitted to SIAM Journal of Numerical Analysis, December 2017




  
  
  


 
  
DefinitionDefinition[section]
theoremTheorem[section]
theorem1Theorem[section]
theorem2Theorem[section]
proposition[theorem]Proposition
lemma[theorem1]Lemma
remark[theorem2]Remark
corollary[theorem]Corollary
example[theorem]Example
hypothesisHypothesis  
  
 
  
  
   








plain


 
 
 Adaptive PBDW approach to state estimation: noisy observations; user-defined update spaces
    Ralph Brinks
Hiller Research Center for Rheumatology
University Duesseldorf
    
==========================================================================================


 
 

Yvon Maday^1,2,
Tommaso Taddei^1



^1   Sorbonne Universités,
Laboratoire Jacques-Louis Lions, France       taddei@ljll.math.upmc.fr,
maday@ann.jussieu.fr 
 
^2 Brown University, Division of Applied Mathematics, USA
       yvon_jean_maday@brown.edu 




 We provide a number 
of extensions and further interpretations of  the  Parameterized-Background Data-Weak (PBDW) formulation,  a real-time and in-situ Data Assimilation (DA) framework for physical systems modeled by parametrized Partial Differential Equations (PDEs),   proposed in [Y Maday, AT Patera, JD Penn, M Yano, Int J Numer Meth Eng, 102(5), 933-965].
  Given M noisy measurements of the state,
  PBDW seeks an approximation of the form
  u^⋆ = z^⋆ + η^⋆, where 
  the background z^⋆  
  belongs to a N-dimensional background space informed by a parameterized mathematical model, 
and the update η^⋆  belongs to a M-dimensional update space
informed by the experimental observations.
The contributions of the present work are threefold:
first, we extend the adaptive formulation proposed in 
[T Taddei, M2AN, 51(5), 1827-1858] to general linear observation functionals, to effectively deal with noisy observations;
second, we consider an user-defined choice of the update space,  to improve convergence with respect to the number of measurements; third, we propose an a priori error analysis for general linear functionals in the presence of noise,
to identify the different sources of state estimation error and ultimately motivate the adaptive procedure. We present results for two  synthetic model problems  
in Acoustics, to illustrate the elements of the methodology and to prove its effectiveness.
We further present results for a synthetic problem in Fluid Mechanics to demonstrate the applicability of the approach to vector-valued fields.


Keywords:
variational data assimilation; parametrized partial differential equations; model order reduction; design of experiment.

 

§ INTRODUCTION


Data Assimilation (DA) refers to the process of integrating information coming from a (possibly parameterized) mathematical model with experimental observations, for prediction.
State estimation is a particular DA task in which the Quantity of Interest (QOI) is the state u^ true of a physical system over a domain of interest Ω⊂ℝ^d. 
In this work, we propose a number of extensions, and further interpretations, of  the  Parameterized-Background Data-Weak (PBDW) approach to state estimation,
first presented  in <cit.>.


As in <cit.>, we denote by { y_m }_m=1^M the set of experimental measurements, and we denote by u^ bk(μ) ∈𝒳 the solution to the  parameterized best-knowledge (bk) mathematical model for the parameter value μ∈𝒫^ bk,
G^ bk,μ(u^ bk(μ) ) = 0.
 Here, the space 𝒳=𝒳(Ω) is a suitable Hilbert space defined over Ω, 
 G^ bk,μ(· ) denotes 
the  parameterized bk mathematical model associated with the physical system, and
𝒫^ bk⊂ℝ^P reflects the uncertainty in the value of the parameters of the model.
We here consider  measurements  of the form
y_m = ℓ_m^o(u^ true) + ϵ_m, where 
ℓ_m^o is a  local average of the state about the location
x_m^ obs in Ω, and ϵ_m reflects the observational noise.
On the other hand, the uncertainty in the parameters of the model leads to the definition of the  bk  manifold ℳ^ bk := {u^ bk(μ)|_Ω:  μ∈𝒫^ bk}⊂𝒳 which collects the solution to the parameterized bk model for all values of the parameter in 𝒫^ bk, restricted to the domain of interest Ω.
Since in practice the model is affected by non-parameterized (unanticipated) uncertainty, 
the true field does not in general belong to ℳ^ bk.

The key idea of  the PBDW formulation is to  seek an approximation u^⋆= z^⋆ + η^⋆ to the true field u^ true employing projection-by-data. The first contribution to u^⋆, z^⋆∈𝒵_N, is the deduced background estimate.
The linear N-dimensional space 𝒵_N ⊂𝒳 is  informed by the bk manifold ℳ^ bk:
it thus encodes — in a mathematically-convenient way — the available prior knowledge about the system coming from the model.
The second contribution to  u^⋆, η^⋆∈𝒰_M, is the update estimate.
The linear M-dimensional space 𝒰_M is the span of the Riesz representers of the M observation functionals {ℓ_m^o }_m=1^M: 𝒰_M  improves the approximation properties of the search space associated with the state estimation procedure.
From a modelling perspective, the deduced background z^⋆ addresses the parameterized uncertainty in the bk   model,
while the update  addresses the non-parametric or unanticipated uncertainty.
In <cit.>, for  perfect measurements (ϵ_m = 0 for m=1,…,M),  
the pair (z^⋆, η^⋆) ∈𝒵_N ×𝒰_M is computed by searching for η^⋆ of minimum norm subject to the observation constraints ℓ_m^o(z^⋆ + η^⋆) = y_m for m=1,…,M.
In <cit.>,
for  pointwise noisy measurements,
the pair (z^⋆, η^⋆) is computed by solving a suitable Tikhonov regularization of the original constrained minimization statement proposed in <cit.>.









PBDW provides some new contributions.
First, the variational formulation facilitates the construction of a priori error estimates, which might guide the optimal choice of the experimental observations.
Second, the background space 𝒵_N accommodates anticipated uncertainty associated with the parameters of the model in a computationally-convenient way;
the construction of 𝒵_N based on ℳ^ bk relies on the application of  parametric Model Order Reduction (pMOR) techniques.  
Third, unlike standard least-squares methods, PBDW provides a mechanism — the update η^⋆ — to correct the deficiencies of the bk model. 
Finally,  projection-by-data, as opposed to projection-by-model, simplifies the treatment of uncertainty in boundary conditions, particularly
when the mathematical model is defined over a domain Ω^ bk that strictly contains Ω (<cit.>). 
We remark that several of these ingredients have appeared in different contexts: we refer to 
<cit.>
 for a thorough overview of the links between the PBDW formulation and other DA techniques, and to 
 <cit.> for a thorough introduction to Data Assimilation.

In this work, we propose 
three  contributions to the original PBDW formulation.



  * We extend the adaptive formulation proposed in 
<cit.> for pointwise measurements to general linear observation functionals.
The approach reads as a Tikhonov regularization of the  constrained minimization statement  in
<cit.>, and can also be interpreted as a convex relaxation of the Partial Spline Model (<cit.>).

  * In <cit.>, the update space 𝒰_M is induced by the particular inner product chosen for 𝒳
(variational update); in this work, we propose to first choose the space 𝒰_M, and then recover the variational formulation through the definition of a suitable inner product (user-defined update).
We demonstrate that our choice reduces the offline costs, and improves the approximation properties of the update space for smooth fields. We emphasize that the idea of selecting the approximation space before introducing the inner products is widely used in the kernel methods' literature
(see, e.g., <cit.>)  for pointwise measurements, and has already been exploited in <cit.>.

  * We propose an a priori error analysis for general linear functionals in the presence of noise.
First, we present a bound for the deterministic error between u^ true and the PBDW solution u^ opt, fed with perfect measurements.
Then, we present a bound for the stochastic error between the possibly regularized PBDW solution u^⋆ fed with imperfect (noisy) measurements,  and the PBDW solution u^ opt, fed with perfect measurements. 
The idea of identifying different sources of error and deriving distinct bounds for each source is the same exploited in <cit.>, for the estimate of the state estimation error based on local measurements of the state.


The outline of the paper is as follows.
In section <ref>, we present the methodology: 
we review  the derivation of the adaptive PBDW formulation for noisy measurements as discussed in <cit.>; 
we introduce the user-defined update space, and we discuss its variational interpretation; 
we discuss the selection of the observation functionals; 
we present an adaptive procedure for the selection of the hyper-parameter associated with the Tikhonov regularizer;
and we discuss the extension to vector-valued problems.
In section <ref>, we present an a priori error analysis for noisy measurements for the state estimation error.
In section <ref>   we present results for two  synthetic model problems  in Acoustics, to  illustrate the elements of the methodology and to prove its effectiveness.
 Finally, in section <ref>, we consider  a synthetic Fluid Mechanics model problem to demonstrate the applicability of PBDW to vector-valued problems.







§ FORMULATION




 §.§ Preliminaries



Given the Lipschitz domain Ω⊂ℝ^d,   we introduce the Hilbert space 𝒳
defined over  Ω; we 
 endow  𝒳 with the inner product (·, ·) and the induced norm · = √( (·, ·  )).
 We further denote by 𝒳' the dual space of 𝒳.
 For any closed linear  subspace 𝒬⊂𝒳, we denote by Π_𝒬 : 𝒳→𝒬 the orthogonal projection operator onto 𝒬, and we denote by 𝒬^⊥ its orthogonal complement.
Given  the linear functional ℓ∈𝒳', we denote by R_𝒳ℓ∈𝒳 the corresponding Riesz element,
(R_𝒳ℓ, f) = ℓ(f) for all f ∈𝒳.
 
 Given a random variable X, we denote by 𝔼[X] and by 𝕍[X] the mean and the variance, where 𝔼 denotes expectation.
We denote by  X ∼𝒩(m, σ^2) a Gaussian random variable with mean   m and variance  σ^2. Similarly, we  denote by  X ∼ Uniform(Ω) an uniform random variable over Ω.
Furthermore, we refer to an arbitrary random variable ε such that 𝔼[ε]=0 and 𝕍[ε] = σ^2 using the notation ε∼ (0,σ^2).

We state upfront that in this section,  we limit ourselves to real-valued problems; however, the discussion can be trivially extended to the complex-valued case. Furthermore, in sections <ref> — <ref>, we consider scalar fields; then, in section <ref>, we discuss the extension to vector-valued problems.
 


 §.§ Derivation of the PBDW statement



As explained in the introduction, we aim to estimate the deterministic state u^ true∈𝒳 over the domain of interest Ω⊂ℝ^d.
We shall afford ourselves two sources of information:
a bk mathematical model 

    G^ bk,μ(u^ bk(μ)) =0,
      μ∈𝒫^ bk

defined over a domain Ω^ bk that contains Ω;
and M experimental observations 
y_1, …, y_M such that

    y_m= ℓ_m^o(u^ true) + ϵ_m,     
     m=1,…,M,

where ℓ_1^o, …, ℓ_M^o∈𝒳' are suitable observation functionals,
 and 
{ϵ_m }_m=1^M are unknown disturbances caused by either systematic error in the data acquisition system or experimental random noise.
Here, 
𝒫^ bk⊂ℝ^P is a confidence region for the true values of the parameters of the model.
We further introduce the bk manifold ℳ^ bk = { u^ bk (μ)|_Ω: μ∈𝒫^ bk} associated with the solution to the parameterized model for all values of the parameters in 𝒫^ bk, restricted to the domain of interest Ω.

In order to combine the parameterized model with the experimental observations, 
Wahba proposed the following generalization  of the 3D-VAR 
(<cit.>)
statement for parameterized background, known as Partial Spline Model  (<cit.>):
find the state estimate
u_ξ^⋆ = u^ bk(μ_ξ^⋆)  + η_ξ^⋆
such that

    (μ_ξ^⋆, η_ξ^⋆ )
     :=
      argmin_(μ,η) ∈𝒫^ bk×𝒰 ξη^2
      
     +
      V_M (
     ℒ_M(u^ bk(μ) + η) 
     - 𝐲);

 where 
ξ >0 is a  tuning hyper-parameter,
 𝐲=[y_1,…,y_M] is the vector of experimental observations, ℒ_M(u)=[ℓ_1^o(u), …,ℓ_M^o(u)], and
 V_M: ℝ^M →ℝ_+ is a suitable strictly convex loss function with minimum in 0 that will be specified later.
 Finally, 𝒰 is 
the search space:  𝒰 is  a closed linear subspace 
contained in 𝒳
that will be specified in the next section.
We observe that    (<ref>) is non-convex in μ;
furthermore, evaluations of the map μ↦ u^ bk(μ) involve the solution to the bk model. Therefore, (<ref>) is
not suitable for real-time computations.

If we introduce the rank-N approximation (<cit.>) of the bk field u^ bk(μ),

    u_N^ bk|_Ω(x,μ)
    =
    ∑_n=1^N 
     
    z_n(μ)  ζ_n(x),
        
    x ∈Ω,
      μ∈𝒫^ bk,

we can approximate statement (<ref>) as

    (μ_ξ^⋆, η_ξ^⋆ )
     :=
      argmin_(μ,η) ∈𝒫^ bk×𝒰 ξη^2
      
     +
      
    V_M
     (
    ℒ_M
    ( ∑_n=1^N  z_n(μ) ζ_n   +  η)
    -
    𝐲).

Assuming that we are not interested in the estimate of the parameter μ_ξ^⋆, (<ref>) is equivalent to:

    (z_ξ^⋆, η_ξ^⋆)
    =
     argmin_  ( z, η  ) ∈Φ_N ×𝒰 ξη^2
      
     +
      
    V_M
     (
    ℒ_M
    ( ∑_n=1^N  z_n ζ_n   +  η)
    -
    𝐲)

where Φ_N = { [z_1(μ),…, z_N(μ)]: μ∈𝒫^ bk}⊂ℝ^N  is the image of the parameter-dependent coefficients, and
the corresponding state estimate is given by 
u_ξ^⋆ = ∑_n=1^N
(𝐳_ξ^⋆)_n ζ_n + η_ξ^⋆.

The set Φ_N in (<ref>) is non-convex; therefore, we convexify 
(<ref>) by minimizing the objective over the convex approximation of Φ_N, 
Φ_N:

    (z_ξ^⋆, η_ξ^⋆)
    =
     argmin_  ( z, η  ) ∈Φ_N ×𝒰 ξη^2
      
     +
      
    V_M
     (
    ℒ_M
    ( ∑_n=1^N  z_n ζ_n   +  η)
    -
    𝐲)

We refer to (<ref>) as to the PBDW formulation: the formulation generalizes the statements introduced in 
<cit.>.

In this work, we study the following instance of the general PBDW formulation:

    (z_ξ^⋆, η_ξ^⋆)
    =
     argmin_  ( z, η  ) ∈ℝ^N  ×𝒰 ξη^2
      
     +
      1/M ℒ_M
    ( ∑_n=1^N  z_n ζ_n   +  η)
    -
    𝐲_2^2,

which can also be reformulated as follows: 

    (z_ξ^⋆, η_ξ^⋆) :=
     arginf_  (z,η)  ∈𝒵_N ×𝒰 
    J_ξ(z, η):=
    ξη^2
     
    +
     1/M ℒ_M
    (z   +  η)
    -
    𝐲_2^2,

where 𝒵_N = {ζ_n}_n=1^N ⊂𝒳 is called the background space.
In the limit ξ→ 0^+, the formulation reduces to 
(<cit.>): 

    (z_ξ^⋆, η_ξ^⋆) :=
     arginf_  (z,η)  ∈𝒵_N ×𝒰 η,
        subject   to ℒ_M
    (z   +  η)
    = 
    𝐲;

which corresponds to the original formulation proposed in 
<cit.>.

Some comments are in order.
The choice of the loss function V_M should be based on the expected experimental disturbances: the particular form considered here is appropriate for homoscedastic noise
(ϵ_1,…,ϵ_M are independent realizations of the random variable ϵ∼ (0,σ^2)).
Furthermore, as observed in 
<cit.> in a related framework, the relaxation Φ_N = ℝ^N might lead to stability issues for N ≃ M. We refer to a future work for the analysis of more general losses and for tighter relaxations Φ_N of Φ_N.
The hyper-parameter ξ is selected adaptively, based on hold-out validation as in <cit.>;
for completeness, in section <ref>,
we summarize the validation procedure.




 §.§ Choice of 𝒰: variational and user-defined update spaces



We consider two choices for the search space 𝒰: 
(i) 𝒰= 𝒳; and
(ii) 𝒰  =  span{ψ_m }_m=1^M where 
ψ_1,…,ψ_M ∈𝒳 are linearly independent user-defined functions satisfying the unisolvency condition:

    ψ∈𝒰_M,
     ℓ_m^o(ψ) = 0,    m=1,…,M
      ⇔  ψ≡ 0.

The first choice is considered in the original PBDW papers, while the second choice is proposed here for the first time.

As shown in <cit.> for (<ref>), for 𝒰= 𝒳, the update η_ξ^⋆ belongs to 
the update space 𝒰_M =  span{ R_𝒳ℓ_m^o}_m=1^M;
we refer to this choice of the  update space 𝒰_M (or equivalently of 𝒰) as to variational update.
If ℓ_1^o,…,ℓ_M^o are linearly independent in 𝒳', it is easy to verify that  
the space 𝒰_M satisfies (<ref>).
Computation of the  update space 𝒰_M requires the solution to M
 variational problems, and depends on the choice of the inner product (·, ·) of 𝒳,
 and on the observation functionals ℓ_1^o,…,ℓ_M^o.

We refer to the second choice 
 𝒰_M = 𝒰 =  span{ψ_m    }_m=1^M  as user-defined update.
The functions ψ_1,…,ψ_M are chosen based on approximation considerations, to guarantee fast convergence with respect to the number of measurements.
We observe that 𝒰_M mildly depends on the observation functionals — it should satisfy (<ref>) — but it does not depend on the choice of the inner product (·, ·) of 𝒳. 
In particular, computation of the user-defined update does not in principle require the  solution to  any variational problem: 
we might thus rely on user-defined updates to simplify the implementation of the method.
It is possible to recover a variational interpretation for this update space for a suitable inner product of 𝒳: we discuss this point in  section <ref>.



  §.§.§ Practical choice of the update space


As explained in the introduction, we are interested in experimental observations 
that can be modelled by linear functionals 
of the form

    ℓ_m^o(u)
    =
    ℓ(
    u, x_m^ obs, r_ w)
    = 
    C(x_m^ obs)
    ∫_Ω ω(
     1/r_ w x - x_m^ obs_2   
    ) 
        u(x)   dx,

where r_ w reflects the filter width of the transducer,  
x_m^ obs denotes the  transducer location, and ω describes the local averaging process performed by the experimental device.
For this class of observation functionals, we consider updates of the form

    𝒰_M =  span{ϕ(·, x_m^ obs )
    }_m=1^M.

where ϕ(·, x_m^ obs) = R_𝒳ℓ(·,   x_m^ obs, r_ w ) (variational update), or ϕ might be an user-defined function.
We refer to ϕ as to update generator.

If 𝒳 is the free space H^s(Ω) for some s ≥ 1, we might consider the generator

    ϕ(·, x) =
    Φ(  · - x _2 ),

where  Φ: ℝ_+ →ℝ_+ is a positive definite kernel Φ (see, e.g., <cit.>).
Computation of (<ref>)
does not require the solution to any variational procedure, and for small values of r_ w it leads to superior approximation properties  compared to the standard H^1 variational update, as we will demonstrate in the numerical results.
Another potential choice is to consider:

    ϕ(·, x) = R_𝒳ℓ(·,   x, R_ w ),

for some R_ w > r_ w.
The generator (<ref>)  requires the solution to M variational problems during the offline stage (as opposed to (<ref>));
on the other hand, it permits the imposition of homogeneous Dirichlet boundary conditions[
Homogeneous Dirichlet boundary conditions could  be imposed using the generator (<ref>) by adding  measurements on the Dirichlet boundary:
this approach increases the dimensionality M of the update space, and thus it increases the online cost associated with the computation of the state estimate.
].



 §.§ Variational interpretation of the user-defined space, and well-posedness analysis



Given the update space 𝒰_M ⊂𝒳, we introduce the interpolation operator
ℐ_M: 𝒳→𝒰_M s.t. ℓ_m^o(  ℐ_M(u)   ) = ℓ_m^o(u) for all u ∈𝒳, m=1,…,M.
Then, we define the symmetric bilinear form


    ((u,v))
    =
    (  u - ℐ_M(u),   v- ℐ_M(v) ) 
     
    +
     
    (   ℐ_M(u),  ℐ_M(v) ).

The bilinear form (<ref>) and the induced (semi-)norm 
· = √((·, ·)). satisfy the following properties.



Let 𝒰_M =  span{ψ_m }_m=1^M be an M-dimensional space satisfying (<ref>),
where ψ_1, …, ψ_M are orthonormal. Then, the following hold.


  * u 
=
  u       ∀  u ∈𝒰_M.

  * Given the orthonormal basis ψ_1,…,ψ_M of 𝒰_M, we define 
𝕃_η∈ℝ^M,M such that
(𝕃_η)_m,m'= ℓ_m(ψ_m'). Then, we can rewrite (<ref>) as follows:

    ((u,v))
    =
    (  u - ℐ_M(u),   v- ℐ_M(v) ) 
     
    +
     ∑_m,m'=1^M  ℓ_m^o(u)  ℓ_m'^o(v)  𝕎_m,m',

where 𝕎 := 𝕃_η^-T𝕃_η^-1.

  * The bilinear form ((·, · ))  induces a norm over 𝒳. More precisely, the following estimate holds: 

    1/2  u  ^2
    ≤u^2
    ≤(
    2 + 3 ℐ_M _ℒ(𝒳)^2
    )   u  ^2,
      ℐ_M _ℒ(𝒳)
    =
    sup_v ∈𝒳 ℐ_M(v)   / v .


  * Let ϕ_m = ∑_p=1^M (𝕃_η)_m,pψ_p. Then, (( ϕ_m, v  )) = ℓ_m^o(v) for all v ∈𝒳.

  * If 𝒰_M =  span{  R_𝒳ℓ_m^o }_m=1^M,
then   ((·, ·)) = (·, ·).

  * Π_𝒰_M^· u = ℐ_M(u) for all u ∈𝒳.




If u ∈𝒰_M, ℐ_M(u)=u.
Therefore, 

    u^2
    =
    (u - ℐ_M(u), u - ℐ_M(u)) + 
    (ℐ_M(u), ℐ_M(u))
    =
    (u,u)
    =
      u ^2,

which proves the first statement.


Observing that ℐ_M(u) = ∑_m=1^M  ( 𝕃_η^-1ℒ_M(u) )_m ψ_m, and recalling that 
ψ_1,…,ψ_M are orthonormal, we find 

    [                                                                                  (ℐ_M(u), ℐ_M(v))                           =
    ∑_m,m'=1^M  
     ( 𝕃_η^-1ℒ_M(u) )_m
      ( 𝕃_η^-1ℒ_M(v) )_m' 
      (ψ_m, ψ_m');                                                                                                   =
      ∑_m=1^M  
     ( 𝕃_η^-1ℒ_M(u) )_m
      ( 𝕃_η^-1ℒ_M(v) )_m
    =
    ∑_p,p'=1^M
     ℓ_p^o(u)  ℓ_p'^o(v)  𝕎_p,p',;                                                                                                   ]
 
which proves the second statement.

We now prove (<ref>) (cf. statement 3).
First, we find that

    u ^2
    ≤ 2 (
      u - ℐ_M(u) ^2
    +
    ℐ_M(u) ^2
    )
    =
    2 u^2,

where the first inequality follows from 
the triangular inequality and
(a+b)^2 ≤ 2 a^2 + 2b^2, while the second identity follows 
 directly 
from the definition of ((·, ·)).
Second, exploiting  the fact that the interpolation operator is continuous in 𝒳, we find

    u^2
    =
      u - ℐ_M(u) ^2
    +
    ℐ_M(u) ^2
    ≤
    2  u ^2
    +
    3 ℐ_M(u)   ^2
    ≤(
    2 + 3 ℐ_M _ℒ(𝒳)^2
    )    u ^2.

Combining the latter two estimates, we obtain (<ref>).

We prove the fourth statement.
Since ϕ_m ∈𝒰_M, ϕ_m - ℐ_M(ϕ_m) ≡ 0. Therefore,

    ((  ϕ_m, v)) = ∑_p,p'=1^M  ℓ_p^o(ϕ_m)  ℓ_p'^o(v) 𝕎_p,p'.

By definition of the interpolation operator, 
we have 
ℓ_p'^o(v) = ℓ_p'^o(  ℐ_M(v)  ) for p'=1,…,M, and 
ℐ_M(v) = ∑_k'=1^M (
𝕃_η^-1ℒ_M(v)
)_k'ψ_k'. Then, we find

    [                                                                           (( ϕ_m, v   )) =         ∑_k,k'=1^M ∑_p,p'=1^M
     ( 𝕃_η)_m,kℓ_p^o(ψ_k) ℓ_p'^o(ψ_k')
    𝕎_p,p'( 𝕃_η^-1ℒ_M(v) )_k';                                                                                            =
    ∑_k,q=1^M (𝕃_η)_m,k(
    𝕃_η^-1)_k,qℓ_q^o(v)
    =
    ∑_q=1^M
     (
    𝕃_η𝕃_η^-1)_m,qℓ_q^o(v)
    =
    ℓ_m^o(v), ]

for all  v ∈𝒳.
Thesis follows.

In order to prove the fifth statement,
we observe that, 
for 𝒰_M =  span{ R_𝒳ℓ_m^o }_m=1^M,
ℐ_M(u) = Π_𝒰_M u.
Then, exploiting the properties 
of ((·, ·)) shown so far
and the projection theorem, we find

    u^2  =
     u -  ℐ_M(u)^2 + ℐ_M(u)^2
    =
    Π_𝒰_M^⊥(u)^2
    +  Π_𝒰_M (u)^2
    =
     u ^2.


The proof the sixth statement follows directly from the fourth property:

    (( Π_𝒰_M^· u, v)) =  (( u, v))  ∀  v ∈𝒰_M  ⇔ ℓ_m^o(  Π_𝒰_M^· u    ) = ℓ_m^o(v),  m=1,…,M.

The latter completes the proof.





The norm ℐ_M _ℒ(𝒳) is the Lebesgue constant in the 𝒳 norm, and it is the inverse of the inf-sup constant (<cit.>):

    γ_M:=
    inf_w ∈𝒰_Msup_v ∈𝒲_M (w,v)/ w   v ,
      𝒲_M
    =
     span{ 
    R_𝒳ℓ_m^o
     }_m=1^M.

The Lebesgue constant depends on the triple 𝒪= ({ℓ_m^o }_m=1^M,  𝒰_M, (𝒳, ·)   ).
For certain choices of the triple 𝒪, it   is possible to determine the asymptotic behavior of ℐ_M _ℒ(𝒳):
to provide a concrete example, 
we  refer  to <cit.> for an important result concerning the behavior of ℐ_M _ℒ(𝒳) for
one-dimensional polynomial interpolation.




Lemma <ref> can be exploited to 
recover an infinite-dimensional variational interpretation of the PBDW statement with user-defined update.
The proof is a straightforward consequence of <cit.>,  and is here omitted.



Let 𝒰_M be an M-dimensional space satisfying (<ref>); then, the PBDW solution to (<ref>) 
for 𝒰= 𝒰_M
solves the following problem:

    (z_ξ^⋆, η_ξ^⋆) :=
     arginf_  (z,η)  ∈𝒵_N ×𝒳 ξη^2
    +
    1/Mℒ_M  (z   +  η) -  𝐲_2^2,

where · is the norm induced by the inner product ((·, ·)) of 𝒳 defined in 
(<ref>).

If the inf-sup constant (<cit.>)

    β_N,M
    = inf_z ∈𝒵_Nsup_q ∈𝒰_M ((z,q))/zq

is strictly positive, then the solution 
(z_ξ^⋆, η_ξ^⋆)
to  (<ref>) is unique, and solves the following saddle-point problem:

    {[ 2ξ ((η_ξ^⋆, q)) + 
    2/M ∑_m=1^M
     (   ℓ_m^o(  z_ξ^⋆  +  η_ξ^⋆   )  - y_m   )  ℓ_m^o(q)                                                                                   =0                                                                           ∀  q ∈𝒰_M;;                                                                         ((η_ξ^⋆, p))                                                                                   =0                                                                          ∀  p  ∈𝒵_N.;                                                                                      ].
  
Furthermore,  the estimate u_ξ^⋆ belongs to the space:

    u_ξ^⋆∈𝒵_N ⊕𝒵_N^⊥, ·∩𝒰_M
    =
    {
    z + η: 
    z ∈𝒵_N,
      η∈𝒰_M,
      
    ((η, q)) = 0    ∀ q ∈𝒵_N
    };

where 
𝒵_N^⊥, · denotes the orthogonal complement  of
𝒵_N with respect to the · norm.



In the variational approach, we first choose an inner product (·, ·) for 𝒳, and then we appeal to an high-fidelity solver to generate the update space;
on the other hand, in the user-defined  approach, we first choose the space 𝒰_M, and then we use it to define the inner product ((·, ·)) (<ref>).
We remark that the second approach is used in the kernel methods' literature (see, e.g., <cit.>) for pointwise measurements. Given a positive definite kernel is possible to characterize the Sobolev regularity of the resulting ambient (native) space: the characterization of the properties of the ambient space can be then exploited 
to prove a priori error bounds for the state estimation error, and estimate the asymptotic convergence rate as M →∞.
On the other hand, the  construction presented in this section does not allow us to characterize the smoothness of the space induced  by the norm ·
in the limit M →∞, unless  ℐ_M _ℒ(𝒳) is bounded for M →∞.





 §.§ Algebraic formulation



Given the fields
z ∈𝒵_N and η∈𝒰_M, we introduce the vectors 𝐳∈ℝ^N and η∈ℝ^M such that
z = ∑_n=1^N z_n ζ_n, and η = ∑_m=1^M η_m ψ_m,
where ζ_1,…,ζ_N and
ψ_1,…,ψ_M are orthonormal bases
(with respect to ·)
 of 𝒵_N and 𝒰_M, respectively.
 We remark upfront that the derivation is independent of the particular construction of the update space.
Then, we can rewrite (<ref>) as a discrete optimization problem for the coefficients:

    (𝐳_ξ^⋆,  η_ξ^⋆)
    =
     argmin_  (𝐳,  η) ∈ℝ^N ×ℝ^M   ξη_2^2
     
    +
     1/M𝕃_ηη + 
     𝕃_z𝐳 - 𝐲_2^2;

where 𝕃_z ∈ℝ^M,N, and 
𝕃_η∈ℝ^M,M are given by
(𝕃_z)_m,n = ℓ_m^o(ζ_n),
(𝕃_η)_m,m' = ℓ_m^o(ψ_m').
We observe that   (<ref>) implies that 𝕃_η is invertible.
By differentiating the objective function with respect to 𝐳 and η, we obtain

    {[ (ξ M  𝕀 +   𝕃_η^T 𝕃_η     ) η_ξ^⋆
    +  𝕃_η^T 𝕃_zz_ξ^⋆
     =   𝕃_η^T    y;;                      𝕃_z^T 𝕃_ηη_ξ^⋆
    +  𝕃_z^T 𝕃_zz_ξ^⋆
     = 𝕃_z^T    y.;                                                                      ].

If we multiply (<ref>)_1 by 𝕃_z^T  𝕃_η^-T, and we exploit (<ref>)_2 we obtain

    ξ M    𝕃_η^T  𝕃_η^-Tη_ξ^⋆ + 
    𝕃_z ^T 𝕃_ηη_ξ^⋆ +
    𝕃_z^T 𝕃_zz_ξ^⋆
     =
    𝕃_z^T    y 
       =
    𝕃_z^T 𝕃_ηη_ξ^⋆ +  𝕃_z^T 𝕃_zz_ξ^⋆⇒𝕃_z^T 𝕃_η^-Tη_ξ^⋆  = 0

Finally, if we introduce η̃_ξ^⋆ =  
𝕃_η^-Tη_ξ^⋆ and we multiply (<ref>)_1 by 𝕃_η^-T, we obtain

    {[ ( ξ M  𝕀 + 𝕃_η𝕃_η^T) η̃_ξ^⋆ 
    + 𝕃_zz_ξ^⋆ = y,;                          𝕃_z^T  η̃_ξ^⋆  = 0;;                                              ].
        η_ξ^⋆ =  𝕃_η^Tη̃_ξ^⋆.

System (<ref>) can be used to efficiently compute the PBDW solution.


(spectrum of the PBDW system)
Following the argument in 
<cit.>, we observe that
the saddle-point system (<ref>)  is congruent to the block-diagonal matrix 
[
[                  𝕃_η𝕃_η^T  + ξ M 𝕀                                  0;                                  0 - 𝕃_z^T (𝕃_η𝕃_η^T   + ξ M 𝕀)^-1𝕃_z;                                    ]]. 
If the update satisfies (<ref>), 𝕃_η is invertible, and thus 
𝕃_η𝕃_η^T  + ξ M 𝕀 is positive definite for all ξ≥ 0.

On the other hand, the second block 
- 𝕃_z^T (𝕃_η𝕃_η^T   + ξ M 𝕀)^-1𝕃_z is negative definite if and only if 
the rank of 𝕃_z is equal to N: this condition is equivalent to the positivity of β_N,M in Proposition <ref>.


In view of the analysis of the method, and of the definition of the Greedy procedure for the selection of the transducers' location,  next Proposition  provides a computable expression for the inf-sup constant β_N,M defined in (<ref>), and for the norm ℐ_M _ℒ(𝒳) in (<ref>).




Suppose that ζ_1,…,ζ_N and ψ_1,…,ψ_M are orthonormal in ·.
The inf-sup constant β_N,M is the square root of the minimum eigenvalue of the following eigenproblem:

    𝕃_z^T 𝕎𝕃_z  𝐳_n
    =
    ν_n (
    𝕀 + 
    2 𝕃_z^T  𝕎𝕃_z
    - 2   sym( ℂ^T 𝕃_η^-1𝕃_z )
    )
    𝐳_n,
      
    n=1,…,N;

while the norm ℐ_M _ℒ(𝒳) is the inverse of the square root of the minimum eigenvalue of the following eigenproblem:

    𝕃_η^T 𝕂^-1𝕃_ηη_m =
    λ_m  η_m,   
    m=1,…,M.

Here,
𝕎  = 𝕃_η^-T𝕃_η^-1,
sym(A)=1/2(A+A^T), while the matrices 
ℂ∈ℝ^M,N and  
𝕂∈ℝ^M,M 
are given by
ℂ_m,n = (ψ_m, ζ_n), and 
𝕂_m,m' = (R_𝒳ℓ_m^o, R_𝒳ℓ_m'^o).



Given z = ∑_n=1^N z_n ζ_n, η = ∑_m=1^M η_m ψ_m,
we can write the interpolation operator as
ℐ_M(z) = ∑_k=1^M  (   𝕃_η^-1𝕃_z 𝐳)_k
ψ_k. Then, we obtain

    ((z, η   )) = 𝐳^T 𝕃_z^T 𝕎𝕃_ηη,
      η^2 = η^T η;

and

    z^2 = 
          z  ^2 +  ℐ_M(z)  ^2 - 2 (z, ℐ_M(z))
    + 
    𝐳^T 𝕃_z^T 𝕎𝕃_z  z
    =
    𝐳^T  𝔹 𝐳.

where 
𝔹 = 
𝕀 + 
2 𝕃_z^T  𝕎𝕃_z
-2   sym( ℂ^T 𝕃_η^-1𝕃_z ).
Recalling the definition of 𝕎,
we find

    β_N,M^2
    =
    inf_𝐳∈ℝ^Nsup_η∈ℝ^M   (   
    𝐳^T 𝕃_z^T 𝕎𝕃_ηη)^2
    / ( η^T  η )
     𝐳^T 𝔹𝐳
    =
    inf_𝐳∈ℝ^N 𝐳^T 𝕃_z^T 𝕎𝕃_z𝐳/𝐳^T 𝔹𝐳.

Introducing the Lagrangian multiplier ν, we can write the optimality conditions as

    {[ 𝕃_z^T 𝕎𝕃_z𝐳 
    =
    ν𝔹𝐳;;       𝐳^T  𝔹𝐳  = 1.;                     ].

Thesis follows.

The proof of (<ref>) exploits Remark <ref>, and the same  argument used to prove (<ref>). 
We omit the details.






 §.§  algorithm for the selection of the observation centers



Algorithm <ref> summarizes the Greedy procedure for the selection of the transducers' locations
 x_1^ obs, …, x_M^ obs.
 During the first (stability) stage, we maximize the constant β_N,M in a Greedy manner.
During the second (approximation) stage, we minimize  the fill distance 
h_M := sup_x ∈Ω min_m=1,…, M x - x_m^ obs_2 in a Greedy manner.
We switch from the first to the second stage when the inf-sup constant β_N,M is larger than an user-defined constant tol:
if tol=0, all the centers are selected through the approximation loop,
if tol=1, all centers are selected through the stability loop.
Representative values for tol used in the numerical simulations are 
tol ∈  [0.2, 0.6].

The Greedy procedure was first presented in
<cit.> for the variational update space.
We observe that the stability loop — for variational update — corresponds[
In the SGreedy procedure listed in <cit.>, steps 6 and 7 of Algorithm <ref> are replaced by 
ℓ_m+1^o=  argmin_ℓ∈ℒ | ℓ( z_ min,m - ℐ_M( z_ min,m  )  )|, and
𝒰_m+1 = 𝒰_m ∪ span{ R_𝒳ℓ_m+1^o  },
where ℒ = {ℓ(·, x, r_ w):   x ∈Ω} is the dictionary of available functionals. 
The approach presented here is easier to implement, and less computationally expensive.
Furthermore, the two procedures are asymptotically equivalent for  r_ w→ 0^+, and return similar results if r_ w is small enough compared to the characteristic length-scale of the elements in 𝒵_N.]
 to the SGreedy algorithm proposed in <cit.>; on the other hand,  the  strategy for the approximation step is strongly related to the so-called farthest-first traversal approach to the minimax facility location problem (see, e.g, <cit.>), first proposed by Rosenkrantz et al. in <cit.>.
For the variational update, since β_N,M is a non-decreasing function of M for a fixed value of N (see <cit.>), the stability constant remains above the threshold during the approximation stage.
On the other hand, for the user-defined update, there is in general no guarantee that the inf-sup constant will be above the threshold at the end of the Greedy procedure
since the inner product ((·, ·)) and the induced norm vary with m.
In the numerical experiments
we investigate the behavior  of the inf-sup constant for large values of M.







 §.§ Choice of ξ



The choice of ξ is performed using holdout validation.
We consider two mutually exclusive datasets 𝒟_M^ train = { (ℓ_m^o, y_m)  }_m=1^M and
𝒟_I^ val = { (ℓ_i^o, y_i)  }_i=1^I
;
given the finite-dimensional search space Ξ⊂ℝ_+, we then select ξ^⋆ such that

    ξ^⋆
    =
     argmin_ξ∈Ξ  MSE(I)
    :=
    1/I∑_i=1^I  ( y_i - ℓ_i^o(u_ξ^⋆   ) )^2,

where u_ξ̅^⋆ is the PBDW solution based on the training dataset 𝒟_M^ train, for ξ = ξ̅.

In this work, we consider  validation measurements of the form 
{ y_i =  ℓ(u, x_i^ obs, r_ w )  + ϵ_i  }_i=1^I, 
where ℓ is introduced in (<ref>), and
x_1^ obs, …, x_I^ obs are independent realizations of an uniformly-distributed random variable over Ω, X ∼ Uniform(Ω).
 If r_ w is small, for this choice of the observation centers, and assuming that 
ϵ_1, …,ϵ_I are independent realizations of ϵ∼ (0,σ^2),  it is possible to show that  (see <cit.>)

    𝔼[
     MSE(I)
    ]
    = 
    1/|Ω|  u^ true - u_ξ^⋆_L^2(Ω)^2 + σ^2
    +
    C(r_ w, u^ true - u_ξ^⋆),

where C(r_ w, u^ true - u_ξ^⋆) → 0 as r_ w→ 0^+ if 
∇ (u^ true - u_ξ^⋆) ∈ L^q(Ω), q>d.
Therefore, for sufficiently large values of I, the validation procedure approximately minimizes the L^2 error.


The choice of the number of validation measurements is a trade-off between experimental cost (given by the number of transducers dedicated to validation) and reliability of the validation procedure.
κ-fold cross-validation 
(see, e.g., <cit.> and  <cit.>)
reduces the experimental cost by partitioning the training  dataset 𝒟_M into κ equal-sized subsamples (folds) {𝒟_M^(k)}_k=1^κ. Of the κ folds,  a single fold is retained for testing and the remaining κ-1 folds are used for training. The procedure is then repeated κ times with each of the κ folds used once as the validation dataset.
After having selected the hyper-parameter ξ, the entire dataset is finally used for training.
We emphasize that cross-validation
relies on the assumption that the pairs (x_m^ obs, y_m) 
are independently sampled from the joint distribution of inputs and outputs
(random design).
On the other hand, in our framework the training centers { x_m^ obs}_m=1^M are chosen deterministically 
(fixed design).
For this reason, in this work, we simply consider holdout validation with I=M/2; we refer to a future work for more advanced cross-validation strategies.




 §.§ Extension to vector-valued problems



We can trivially extend the Greedy procedure for the selection of the observation centers to vector-valued fields
u^ true = [u_1^ true, …, u_J^ true] ∈𝒳=𝒳(Ω; ℝ^J).
Assuming that each transducer is able to measure all J components of the true field in all d directions, at each iteration of the Greedy algorithm, we select x_m^ obs such that
(compare with Algorithm <ref>, step 6)

    x_m+1^ obs =  argmax_x ∈Ω  z_ min(x) - ℐ_m (z_ min) (x)  _2.

Then, we add J modes ϕ_1(·, x_m+1^ obs), …, ϕ_J(·, x_m+1^ obs) ∈𝒳 in the update; note that  the subspace 𝒰_m is m × J dimensional).

As regards the choice of the J generators, we might consider the generator


    ϕ_i(·, x) = R_𝒳ℓ_i(·, x, R_ w),
       where ℓ_i(u,x, R_ w) = ℓ(u_i, x,R_ w),
      
    i=1,…,J;

and R_ w≥ r_ w. Alternatively, we might consider the update 

    ϕ_i(·, x) =
    Φ(· - x _2) 𝐞_i,
      
    i=1,…,J;

where Φ is a properly-defined kernel, and 
𝐞_1,…,𝐞_J are the canonical vectors in ℝ^J.


Note that, in the case of incompressible flows considered in section  <ref>,
for a proper choice of the ambient space 𝒳, the  update space (<ref>) is divergence-free, 
while the user-defined update space (<ref>) is not divergence-free.
For this reason, we might resort to matrix-valued divergence-free kernels (see, e.g., <cit.>).
We refer to a future work for the use of matrix-valued divergence-free kernels in the PBDW framework.










§ ERROR ANALYSIS


We present an error analysis for general linear functionals in the presence of noise, based on the identification of two sources:
the first source is related to the finite number of measurements, while the second source is related to the presence of noise.
For simplicity of notation, we consider scalar problems; however, the analysis can be trivially extended to vector-valued fields.
Given the measurements { y_m = ℓ_m^o(u^ true) + ϵ_m  }_m=1^M, 
we define the vector 
𝐮^ opt = [η̃^ opt, 𝐳^ opt] ∈ℝ^M+N corresponding to the solution 
u_ξ=0^ opt to  (<ref>)
fed with perfect observations { y_m^ true = ℓ_m^o(u^ true) }_m=1^M (see  (<ref>)).
On the other hand, we define  𝐮_ξ^⋆ = [η̃_ξ^⋆, 𝐳_ξ^⋆] corresponding to the solution
u_ξ^⋆  to the stabilized formulation (<ref>), fed with  imperfect observations. 


We first present the result for perfect measurements.



Let y_m = ℓ_m^o(u^ true), m=1,…,M, 
let β_N,M>0, and let 𝒰_M satisfy 
(<ref>).
We denote by u_ξ=0^ opt the solution to the PBDW formulation (<ref>).
Then, the following estimates hold:


    u^ true - u_ξ=0^ opt≤1/β_N,M inf_z ∈𝒵_N inf_  q ∈𝒰_M ∩𝒵_N^⊥, ·   u^ true  - z - q   ;


    u^ true - u_ξ=0^ opt≤√(4 + 6 ℐ_M _ℒ(𝒳)^2   )/β_N,M inf_z ∈𝒵_N inf_  q ∈𝒰_M ∩𝒵_N^⊥, ·   u^ true  - z - q .





Thanks to the  variational interpretation of the user-defined update
(see Proposition <ref>), 
proof of (<ref>) follows  
directly from <cit.>.
On the other hand, exploiting 
(<ref>) and
(<ref>), we find

    [                                                                       u^ true - u_ξ=0^ opt ≤√(2) u^ true - u_ξ=0^ opt≤√(2)/β_N,M inf_z ∈𝒵_N inf_  q ∈𝒰_M ∩𝒵_N^⊥, ·   u^ true  - z - q;                                                                                               ≤√(4 + 6 ℐ_M _ℒ(𝒳)^2   )/β_N,M inf_z ∈𝒵_N inf_  q ∈𝒰_M ∩𝒵_N^⊥, ·   u^ true  - z - q   ,;                                                                                            ]

which proves  (<ref>).


If we introduce the matrices and vectors

    𝔸(ξ) = 
    [
    [ ( ξ M  𝕀 + 𝕃_η𝕃_η^T)                  𝕃_z;                𝕃_z^T                    0;                      ]];
      𝐞  =
    [
    [ ϵ; 0;   ]];
      ϵ = [ϵ_1,…,ϵ_M]^T;

we obtain the following identity that links 𝐮_ξ^⋆ to 𝐮^ opt

    𝔸(ξ) 𝐮_ξ^⋆
    =
    𝔸(0) 𝐮^ opt
    +
     𝐞
     =
     𝔸(ξ) 𝐮^ opt
    -
     ξ M 
     [  
    [ η̃^ opt;       0 ]]
     +
     𝐞.

Identity (<ref>) can be used to estimate the bias of our estimate, and the expected error.
For simplicity, we present below estimates for the errors in the coefficients, 
𝐮_ξ^⋆ - 𝐮^ opt;
the estimates can also be extended to compute the state estimation error  in an integral norm of interest. We omit the details.



Suppose that ϵ_1,…,ϵ_M are independent realizations of
ϵ∼ (0,σ^2). Then, the following hold:

    𝔼[  𝐮_ξ^⋆   ]    -    𝐮^ opt_2
    ≤ξ M/s_ min ( 𝔸(ξ)  )   η̃^ opt_2;

and

    𝔼[
    𝐮_ξ^⋆ -  𝐮^ opt_2^2
    ]
    ≤(
    ξ M/s_ min ( 𝔸(ξ)  )   η̃^ opt_2
    )^2
    +
    σ^2
     trace(
    𝔸(ξ)^-1Σ𝔸(ξ)^-T);

where
s_ min( 𝔸(ξ) ) denotes the minimum singular value of 
𝔸(ξ), 
 Σ = 
[
[ 𝕀 0; 0 0;   ]].


Exploiting (<ref>), we find

    𝔸(ξ) 𝐮_ξ^⋆
     =
     𝔸(ξ) 𝐮^ opt
    -
     ξ M 
     [  
    [ η̃^ opt;       0 ]]
     +
     𝐞 ⇒ 𝔼[  𝐮_ξ^⋆   ]    -    𝐮^ opt_2
    ≤ξ M/s_ min ( 𝔸(ξ)  )   η̃^ opt_2;

which is (<ref>).

In order to show (<ref>), we first observe that
(cf. <cit.>):

    𝔼[
    𝔼[𝐮_ξ^⋆  ]  -
    𝐮_ξ^⋆_2^2
    ]
    =
    σ^2
     trace(
     𝔸(ξ)^-1Σ𝔸(ξ)^-T).

Then, 
we find

    [                                                 𝔼𝐮_ξ^⋆  - 𝐮^ opt_2^2 
    =                𝔼[𝐮_ξ^⋆  ]  - 𝐮^ opt_2^2
     +
    𝔼[
    𝔼[𝐮_ξ^⋆  ]  - 𝐮_ξ^⋆_2^2
    ];                                                                         ≤(
    ξ M/s_ min ( 𝔸(ξ)  )   η̃^ opt_2
    )^2
    +
    σ^2
     trace(
    𝔸(ξ)^-1Σ𝔸(ξ)^-1),;                                                                         ]

which is  (<ref>).


We observe that the bound for the mean squared error
(<ref>) is the sum of two contributions:
the former 
— 
ξ M/s_ min ( 𝔸(ξ)  )   η̃^ opt_2
—
involves the accuracy of the background space,
measured by η̃^ opt_2, 
and is monotonically increasing with ξ;
the latter 
—
σ^2
 trace( 𝔸(ξ)^-1Σ𝔸(ξ)^-T)
—
involves the accuracy of the measurements, and is monotonically decreasing with ξ according to our numerical experience.
Therefore, the optimal value of ξ depends on the ratio between “measurement inaccuracy"
and “model  inaccuracy"
σ/η̃^ opt_2:
 if u^ true∈𝒵_N, estimate (<ref>) suggests to pick   ξ→∞; 
 if σ=0, estimate (<ref>) suggests to pick  
 ξ→ 0^+. 
 Since in practice estimates of the ratio σ/η̃^ opt_2
 are rarely available, our analysis suggests that the value of ξ should be chosen adaptively using (cross-)validation techniques such as the one discussed in section <ref>.
These observations are in agreement with <cit.>.


(extension  to  complex-valued problems)
Proposition <ref> can  be trivially extended to the complex-valued case. Assuming that measurements are of the form

    y_m = ℓ_m^o(u^ true) + σ_ reϵ_m^ re  +  iσ_ imϵ_m^ im,
        ϵ_m^ re, ϵ_m^ im i.i.d.∼ (0,1), 
      σ_ re, σ_ im>0,

we find the same estimate for the bias (<ref>), and we find

    [                                                 𝔼[ 𝐮_ξ^⋆ - 𝐮^ opt_2^2    ]
    ≤ (
    ξ M/s_ min( 𝔸̃(ξ) )  η^ opt_2
    )^2
    +
    σ_ re^2 trace( 
    𝔸̃(ξ)^-1Σ^ re𝔸̃(ξ)^-T);                                                                                                                 +
    σ_ im^2 trace( 
    𝔸̃(ξ)^-1Σ^ im𝔸̃(ξ)^-T); ]

where 
𝔸̃(ξ) = [ 
[   Re[𝔸(ξ)] - Im[𝔸(ξ)];   Im[𝔸(ξ)]   Re[𝔸(ξ)];            ]],
 Σ^ re = [ 
[ Σ 0; 0 0;   ]],
 Σ^ im = [ 
[ 0 0; 0 Σ;   ]].










 
  

§ APPLICATION TO ACOUSTICS

 
 
 We illustrate the behavior of the PBDW formulation presented in this paper through the vehicle of two acoustic Helmholtz problems.
Since the objective of this work is to propose user-defined  update spaces that improve convergence with respect to the number of observations (M-convergence), we do not study the effect of the primary approximation provided by the background 𝒵_N: we refer to the PBDW literature for further results concerning the effect of the primary approximation.



 §.§ A two-dimensional model problem

The model problem is the same considered in 
<cit.>, and <cit.>. 



  §.§.§ Problem definition

Given the domain 
Ω = (0,1)^2, we define the acoustic model problem: 

    {[ -(1 +  iϵμ)  Δ u_g(μ)   -  μ^2 u_g(μ) =  μ(  2 x_1^2 + e^x_2)   +  μ g                                                                     Ω,;                                                       ∂_n   u_g(μ) = 0                                                                    ∂Ω,;                                                                        ].

where 
i is the imaginary unit, 
μ>0 is the wave number, ϵ=10^-2 is a fixed dissipation,  and g ∈ L^2(Ω) is a bias term that will be specified later.
Here, the parameter μ>0 constitutes the anticipated, parametric uncertainty in the system, which might model our uncertainty in the speed of sound, while the function g constitutes the unanticipated (non-parametric) uncertainty in the system.


To assess the performance of the PBDW formulation for various configurations,
we define the true field u^ true as
the solution to (<ref>) for some μ^ true∈𝒫^ bk and for the following   choice  of the bias g



    g:=  0.5 ( e^-x_1  +   1.3 cos ( 1.3 π x_2  )  ).

On the other hand, we define the bk manifold as

    ℳ^ bk :=  {  u_g=0(μ):  μ∈𝒫^ bk}.




To measure  performance, we introduce  the  relative  L^2  and H^1 errors  averaged over |𝒫_ test^ bk| = n_ test  fields associated with different choices of the parameter μ:

    E_ avg^ rel( n_ test )   := 1/  n_ test ∑_μ∈𝒫_ test^ bk  u^ true(μ) - u_ξ^⋆(μ)  _⋆/ u^ true(μ)_⋆;
        ·_⋆ = 
    ·_L^2(Ω)  or ·_H^1(Ω).

In all our tests, we consider n_ test=10 equispaced parameters in 𝒫^ bk.

We model the (synthetic) observations by a Gaussian convolution 
with standard deviation r_ w:



    ℓ_m^o(v) 
    =  Gauss (v, x_m^ obs; r_ w)
    =
    C(x_m^ obs)  ∫_Ω  e^-1/2 r_ w^2  x - x_m^ obs_2^2    v(x )   dx,
      
    m=1,…,M;

where C(x_m^ obs) is a normalization constant such that ℓ_m^o(1)  = 1, and
 x_m^ obs denotes the transducer location.
 In order to simulate noisy measurements, we corrupt the synthetic measurements with homoscedastic random noise:
 
    y_ℓ =  ℓ_m^o(u^ true) + 
    σ^ reϵ_ℓ^ re +   iσ^ imϵ_ℓ^ im,
        ϵ_ℓ^ re , ϵ_ℓ^ im  iid∼ 𝒩(0,1);

 where 
 
    σ^ re  =  SNR× std(  { Re[    ℓ_m^o(u^ true)   ]   }_m=1^M   ),
      σ^ im  = 
      SNR× std(  { Im[    ℓ_m^o(u^ true)   ]   }_m=1^M   ),

  and SNR denotes the signal-to-noise ratio in the measurements.





  §.§.§ PBDW spaces


We introduce the ambient space 𝒳= H^1(Ω) endowed with the inner product:

    (u,v) = ∫_Ω  u v̅ + ∇ u ·∇v̅  dx.

Note that (̅·̅)̅ denotes the complex conjugate of (·). The background space 
𝒵_N is built using the Weak-Greedy algorithm (see, e.g., <cit.>):
we refer to <cit.> for further details.

As regards the update space, we consider the variational update associated with the inner product (<ref>), and the user-defined updates 
𝒰_M =  span{ϕ (   λ· - x_m^ obs_2   )   }_m=1^M for 

    {[         ϕ(r) = 1/(1+r^2)^2   inverse   multiquadrics;; ϕ(r) =  (1-r)_+^4 (4r + 1)                     csRBF.;                            ].

The observation centers {  x_m^ obs}_m=1^M are chosen according to  Algorithm <ref>.
As regards the choice of the kernel scale, we here set  λ=1 for inverse multiquadrics, and λ=2 for csRBF. We remark that the optimal value for the kernel scale parameter λ strongly depends on the number of measurements, and also on the characteristic length-scale of the field u^ true - z_ξ^⋆: therefore, adaptation of the parameter λ might improve performance, particularly for large values of M. 
 




  §.§.§ Numerical results


In Figures <ref>(a)-(b)-(c),
we show the behavior of the inf-sup constant β_N,M with M 
 for the three different choices of the update space (and thus for the three choices of the norm ·),
r_ w=0.01 and N=6, and for 
 centers { x_m^ obs}_m=1^M selected by the SGreedy algorithm with tol=0.6;
 to measure performance, we compare the results with the ones obtained using  randomly-generated centers. For the second choice, we average over 35 different random choices of the M centers. 
 We observe that the SGreedy selection of the sensor locations leads to a more stable formulation compared to randomly-generated points.
 In Figures <ref>(d)-(e)-(f),
 we  show the centers { x_m^ obs}_m=1^M selected by the SGreedy algorithm for the three different choices of the update space: 
for all the three choices of the update space the Greedy procedure selects most points along the principal diagonal (0,0) → (1,1).  

 


In Figure <ref>, we show the behavior of β_N=6,M fo three
choices of the update space considered, for M=6,…,150.
We observe that, for the user-defined update, β_N,M is not monotonic increasing with M.
Nevertheless, we do not observe pathologic behaviors of the inf-sup constant as M increases.

 



In Figure <ref>, we show the behavior of the Lebesgue constant ℐ_M _ℒ(𝒳) defined in (<ref>), for the three choices of the update. As expected for the variational update space, the Lebesgue constant is equal to one.
We also observe that ℐ_M _ℒ(𝒳) is significantly larger for inverse multiquadrics than for csRBF: recalling estimate (<ref>), the use of inverse multiquadrics (and more in general of smooth kernels) is appropriate only for smooth fields.


 

Figure <ref> shows the behavior of the relative error E_ avg^ rel
(<ref>) with M for perfect observations, N=6, and two values of r_ w in (<ref>). We observe that the user-defined update based on inverse multiquadrics leads to more accurate state estimates compared to the standard H^1 PBDW, particularly for r_ w=0.01. 
We further observe that  inverse multiquadrics outperform csRBF in all tests considered:
since the true state is smooth in Ω,
the  inverse multiquadric kernel —  which is  C^∞ — exhibits superior approximation properties compared to the C^2  csRBF kernel (see <cit.>).

 

Figure <ref> shows results for the behavior of the L^2 relative error E_ avg^ rel with M for imperfect observations for the user-defined (inverse multiquadrics) update and for the variational H^1 update.
The error E_ avg^ rel is averaged over 25 realizations of the homoscedastic random noise.
More in detail, we here consider three different noise levels SNR; we set N=4, and we consider r_ w=0.01 in (<ref>). Furthermore, in order to select the value of ξ, we resort to the holdout validation procedure outlined in section 
<ref> based on additional I=M/2 measurements, associated with uniformly-generated points in Ω.
The rate of convergence with M seems to weakly depend on the choice of the update;
nevertheless, also for noisy measurements, the use of inverse multiquadrics strongly improves performance.

 


In Figure <ref>  we discuss the interpretation of the regularization parameter ξ.
Towards this end, we show the behavior of the mean squared error
MSE(I) introduced in section <ref>
for two different choices of the true field 
(u^ true = u_g(μ=5.8) and 
u^ true = u_g=0(μ=5.8)), 
  two different noise levels, and
N=5, M=100, I=50.
Note that for g=0 u^ true belongs to the bk manifold.
We observe that the optimal value of ξ increases as the noise increases, and decreases as the best-fit error increases:
this is in good agreement with the discussion in section <ref>, and also with the results in 
<cit.>
for pointwise measurements.
We further observe that,
as  anticipated in  section <ref>,
the error estimator is in good qualitative agreement with the true error
u^ true - u_ξ^⋆_L^2(Ω)^2: therefore, it can be used to select the optimal value of ξ.

 












 


 §.§ A three-dimensional acoustic problem
 
 

  §.§.§ Problem definition

We  consider the three-dimensional  model problem:

    {[ - (1+ ϵ i)  Δ  u_g(μ)
        - (2πμ)^2  u_g (μ) = g                                               Ω;;                                 ∂_n  u_g (μ)
    = 0                                              ∂Ω;;                                                  ].

where ϵ=10^-2,   Ω  = (-1.5,1.5) × (0,3) × (0,3) ∖Ω^ cut,
Ω^ cut = (-0.5,0.5) × (0.25,0.5) × (0,1). 
Figure <ref> shows the geometry. 
In this example, we consider the bk manifold ℳ^ bk = 
{ u^ bk(μ)
=u_g^ bk(μ): μ∈𝒫^ bk = [0.1,0.5] }, and we define the true field as the solution
to (<ref>) for some 
μ^ true∈𝒫^ bk
and g=g^ true, where 

    g^ bk(x) = 10   e^ -  x - p^ bk_2^2;
      
    g^ true(x) = 10   e^ -  x - p^ true_2^2;

and p^ bk = [0,2,1], p^ true = [-0.1,2,1].
Parameterized uncertainty in the system models uncertainty in the input frequency μ;
non-parametric or unanticipated uncertainty is here associated with the incorrect location of the acoustic source (p^ bk≠ p^ true).
Computations are based on a P2 FE discretization with roughly 𝒩 = 16000 degrees of freedom in Ω. 
Figure <ref> shows the bk and true solutions  for two values of μ.



 


As in the previous example, we model the synthetic observations by a Gaussian convolution with standard deviation r_ w, see (<ref>).
For simplicity, in the tests below, we only consider perfect measurements.
Furthermore, we measure performance by computing the relative
L^2  and H^1 errors  E_ avg^ rel
(<ref>) for n_ test=10  different choices of the parameter μ in 𝒫^ bk.




  §.§.§ PBDW spaces


We consider the ambient space 𝒳=H^1(Ω) endowed with the inner product (·, ·) (<ref>). 
Furthermore, the background space 𝒵_N is built using the Weak-Greedy algorithm based on the residual.


As regards the update space, we consider the variational  update associated with the inner product (<ref>), and the user-defined updates 
𝒰_M =  span{ϕ (   λ· - x_m^ obs_2   )   }_m=1^M for ϕ(r) = 1/1+r^2 (inverse multiquadrics) and λ=1.
The observation centers {  x_m^ obs}_m=1^M are chosen according to  Algorithm <ref>.





  §.§.§ Numerical results



In Figure <ref>, we perform the same test of
Figure <ref>.
We compare the behavior of the inf-sup constant β_N=30,M with M using the SGreedy procedure, and using randomly-generated centers. For the second choice, we average over 35 different random choices of the M centers. 
As in the previous example, the Greedy procedure   improves the stability of the PBDW formulation.

 

Figure <ref> shows the behavior of the relative error E_ avg^ rel
(<ref>) with M for perfect observations, N=30, and r_ w=0.02 in (<ref>). As in the previous case, the user-defined update based on inverse multiquadrics leads to more accurate state estimates compared to the standard H^1 PBDW. 

 



 
 

§ APPLICATION TO FLUID MECHANICS

 
  

 §.§ Problem statement


We consider the following model problem[
The configuration is similar to the one considered in <cit.>;
however, while here we consider a 2D domain, in <cit.> the authors consider a more realistic axi-symmetric model.
]:

    {[ - 1/ Re ∇·(
    ∇ u_g + 
    ∇ u_g^T
    )
     
    +
     (u_g ·∇) u_g + ∇ p_g = 0                                                         in Ω;                                                   ∇· u_g = 0                                                         in Ω;                      p_g 𝐧 - 1/ Re (∇ u_g + ∇ u_g ^T) 𝐧  = 0                                         on Γ_ out,1∪Γ_ out,2;                                                  u_g = g 𝐞_1                                                     on Γ_ in;                                                      u_g = 0                              on ∂Γ_ hom=
    Ω∖(   Γ_ in∪Γ_ out);                                                              ].

where the domain Ω is depicted in Figure <ref>(a).
 We then define the bk and true manifolds  as follows:
  
    ℳ^ bk
      =
      {
      u_g( Re):
       Re∈ [50, 350],
        
      g(x_2) = 4 (1-x_2)   x_2
      };

 and 
   
    ℳ^ true
      =
      {
      u_g( Re):
       Re∈ [50, 350],
        
      g(x_2) = 4 (1-x_2) x_2
      (   
      1 + 0.1 sin(2π x_2)
      ) 
      }.

Here,   uncertainty in μ =  Re constitutes the anticipated (parametric) uncertainty in the system, while 
  uncertainty in the inflow condition is the unanticipated (non-parametric) uncertainty.
  
  As in the previous examples, we model synthetic observations by  a Gaussian convolution with standard deviation r_ w=0.01. Furthermore, we resort to a conforming Taylor-Hood P3-P2 Finite Element discretization with 
  𝒩_ u=24038 degrees of freedom.
  

  





 §.§ PBDW spaces


We consider the ambient space 𝒳 = { v ∈ [H_0,Γ_ hom(Ω)]^2:  ∇· v =  0  } endowed with the inner product:

    (u, v)
    = ∫_Ω ∇ u : ∇ v   +   u · v   dx.

The background space is built using Proper Orthogonal Decomposition based on the (·, ·) inner product
(see, e.g., <cit.>).
Figure <ref>(b) shows the behavior of the POD eigenvalues.
On the other hand, we consider the update generators:

    ϕ_i(·, x)
    =
    R_𝒳 ℓ_i(·, x, R_ w),
      ℓ_i(u,x,R_ w)
    =
    ℓ(u_i, x, R_ w),
        
    i=1,2,

for several values of R_ w≥ r_ w.
We observe that the update space is divergence-free; furthermore, it satisfies the no-slip boundary conditions on Γ_ hom.




 §.§ Numerical results


In Figure <ref>, we show the 
performance of the SGreedy procedure for N=5 and R_ w=0.05 (tol=0.3).
Figure <ref>(a) shows the behavior of the inf-sup constant β_N,M with M, 
for the user-defined update with R_ w=0.05:
we observe that the SGreedy procedure outperforms on average the random uniform selector.
Figure <ref>(b)  shows the location of the points selected by the SGreedy (stabilization stage) Algorithm. We observe that  SGreedy selects   points in the boundary layer that develops near the junction:
this can be explained by recalling that the Reynolds number — the parameter associated with the solution manifold —
deeply affects the thickness of the boundary layer.

 

Figure <ref> shows the behavior of the relative error E_ avg^ rel 
(<ref>) with M for perfect measurements, N=5, and three values of R_ width in
(<ref>). We observe that for R_ width=0.01 the user-defined update corresponds to the variational update.
As for the previous test cases, the use of an user-defined update strongly improves performance,  particularly for moderate-to-large values of M.
 
     
 


 

§ CONCLUSIONS

 
 
In this paper, we presented a number of extensions  to the PBDW formulation for state estimation.
First, we proposed a Tikhonov regularization of the original PBDW statement for general linear functionals, which relies on holdout validation, to systematically deal with noisy measurements.
Second, we proposed user-defined update spaces, which guarantee rapid convergence with respect to the number of measurements M and also might not require the solution to M Riesz problems. 
Third, we presented an a priori error analysis that provides insights into the role of the regularization hyper-parameter ξ associated with the penalized formulation.

We identify a number of future research directions to improve the PBDW formulation  and extend its range of applications.
First, we wish to extend the formulation to different models for the experimental noise, and also to probabilistic background.
Towards this end, we wish to rely on the connection between PBDW and PSM established in <cit.>, and recapped in section <ref>.
Second, we wish to extend the formulation to time-dependent problems.
This might be accomplished by exploiting a space-time variational formulation (<cit.>)
to incorporate   the space-time structure of the evolution equations in the empirical
expansion.
 









unsrt










