Randomized residual-based error estimators for parametrized equations
Kathrin Smetana Olivier Zahm Anthony T Patera


1/f Noise Reduction using In-Pixel Chopping in CMOS Image Sensor
    Kapil Jainwal and Mukul Sarkar, Member IEEE
Manuscript received June X, XXXX.
Authors are with the Electrical Engineering Department, Indian Institute of Technology Delhi, 110016 New Delhi, India. (e-mail: KJ: kapiljainwal@gmail.com; MS: msarkar@ee.iitd.ac.in) 


     December 30, 2023
===========================================================================================================================================================================================================================================================================




We propose a randomized a posteriori error estimator for reduced order approximations of parametrized (partial) differential equations. The error estimator has several important properties: the effectivity is close to unity with prescribed lower and upper bounds at specified high probability; the estimator does not require the calculation of stability (coercivity, or inf-sup) constants; the online cost to evaluate the a posteriori error estimator is commensurate with the cost to find the reduced order approximation; the probabilistic bounds extend to many queries with only modest increase  in cost. To build this estimator, we first estimate the norm of the error with a Monte-Carlo estimator using Gaussian random vectors whose covariance is chosen according to the desired error measure, e.g. user-defined norms or quantity of interest. Then, we introduce a dual problem with random right-hand side the solution of which allows us to rewrite the error estimator in terms of the residual of the original equation. In order to have a fast-to-evaluate estimator, model order reduction methods can be used to approximate the random dual solutions. Here, we propose a greedy algorithm that is guided by a scalar quantity of interest depending on the error estimator. Numerical experiments on a multi-parametric Helmholtz problem demonstrate that this strategy yields rather low-dimensional reduced dual spaces. 




  A posteriori error estimation, parametrized equations, projection-based model order reduction, Monte-Carlo estimator, concentration phenomenon, goal-oriented error estimation.




  65N15, 65C05, 65N30, 68Q25, 62G15 






§ INTRODUCTION


Many models for engineering applications, life sciences, environmental issues, or finance depend on parameters which account for variation in the material or geometry but also uncertainty in the data. Often the respective applications require low marginal (i.e. per parameter) computational costs. This is for instance the case in “many query” settings where we require the computation of the solution of the corresponding parametrized equation for many different parameter values. 
Examples for model order reduction techniques that aim at computationally feasible approximations of such parametrized models are tensor-based methods <cit.> and the reduced basis (RB) method <cit.>. In order to ensure say functional safety of a structure, certification of such approximation is of high importance. 
Moreover, bounding the approximation error to get a handle on the uncertainty induced by the approximation is crucial when using it in the context of uncertainty quantification. The subject of this paper is thus certification of approximations to parametrized equations via an a posteriori error estimator for a large number of parameter queries. Our method is also well-suited to real-time contexts. Employing the a posteriori error estimator say within a greedy algorithm to construct the reduced space requires some (minor) modifications, which we will touch on only very briefly in this paper.

One of the most commonly used error estimators for inf-sup stable problems is the product of the dual norm of the residual and the inverse of the inf-sup constant. While the former can usually be computed rapidly, accurate estimation of the inf-sup constant is in general rather costly. For instance, the Successive Constraint Method (SCM) <cit.> computes a parameter-dependent lower bound of the inf-sup constant by employing the successive solution to appropriate linear optimization problems. This procedure is usually computationally demanding and can lead to pessimistic error bounds <cit.>.

In this paper we introduce a random a posteriori error estimator which does not require the estimation of stability constants.
The error estimator features several other desirable properties. First, it is both reliable and efficient at given high probability and often has an effectivity close to one. 
Secondly, the effectivity can be bounded from below and above at high probability with constants selected by the user, balancing computational costs and desired sharpness of the estimator.
Moreover, the presented framework yields error estimators with respect to user-defined norms, for instance the L^2-norm or the H^1-norm; the approach also permits error estimation of linear quantities of interest (QoI).
Finally, depending on the desired effectivity the computation of the error estimator is in general only as costly as the computation of the reduced order approximation or even less expensive, which makes our error estimator strategy attractive from a computational viewpoint.

To derive this error estimator, we consider a Gaussian random vector whose covariance matrix is chosen depending on the respective norm or QoI we wish to estimate. Summing the squares of the inner products of K independent copies of that random vector with the approximation error yields an unbiased Monte Carlo estimator. 
Using concentration inequalities, we control the effectivity of the resulting random error estimator with high probability.
This type of random subspace embedding is typically encountered in compressed sensing <cit.>.
The motivation for using these techniques is to create a high-to-low dimensional map which, in high probability, nearly preserves distances and is thus well-suited for norm estimation.
By exploiting the error-residual relationship we recognize that these inner products equal the inner products of the residual and the dual solutions of K dual problems with random right-hand sides. Approximating the dual problems via projection-based model order reduction yields an a posteriori error estimator of low marginal computation cost. To construct the dual reduced space we introduce a greedy algorithm 
driven by a scalar QoI that assesses how good the fast-to-evaluate a posteriori error estimator approximates the original Monte Carlo estimator. This goal-oriented strategy outperforms standard dual-residual based greedy algorithms or the Proper Orthogonal Decomposition (POD). We emphasize that the dual reduced space so obtained does generally not contain the primal reduced space as a subspace; the intersection can even be empty. Furthermore, the dimension of the dual reduced space can be smaller than the dimension of the primal reduced space.


Our a posteriori error estimator is inspired by the probabilistic error estimator for the approximation error in the solution of a system of ordinary differential equations introduced in <cit.> by Cao and Petzold. To estimate the norm of the error, they employ the small statistical sample method from Kenney and Laub <cit.>, which estimates the norm of a vector by its inner product with a random vector drawn uniformly at random on the unit sphere.
Rewriting that inner product using the error-residual relationship results in an adjoint (or dual) problem with random final time, whose solution is then invoked to estimate the error <cit.>. This approach is extended to ordinary differential equations via a POD by Homescu et al in <cit.> and differential algebraic equations in <cit.>. Also, the effect of perturbations in the initial conditions or parameters on the quality of the approximation of the reduced model is investigated <cit.>. 

In our work we extend these concepts to address the general norms of interest within the PDE context, to explicitly address accurate error estimation for any given parameter value within a finite parameter domain, and to address the limit of many queries.


Randomized methods for error estimation are gaining interest in the reduced order modeling community. 

For instance in <cit.>, randomized techniques are used to speed-up the computation of the dual norm of the residual used as an error indicator. By exploiting the fact that the residual manifold is included in a low-dimensional subspace, the authors need appeal to only a few random samples when constructing the random subspace embedding. Instead, our approach targets the true error which, in contrast to the residual, is in general not exactly included in a low-dimensional subspace for the problems we have at hand. Therefore, in our approach, we use different techniques and we determine the number of random sample we need via the cardinality of the parameter set on which we wish to estimate the error.

In <cit.> a probabilistic a posteriori error bound for linear scalar-valued quantities of interest is proposed, with application in sensitivity analysis. Contrary to the method presented in our work, the right-hand side of the dual problem in <cit.> is the linear functional associated with the QoI and randomization is done by assuming that the parameter is a random variable on the parameter set.

Another application of randomized techniques, in particular randomized numerical linear algebra <cit.>, to (localized) model order reduction is considered in <cit.>: a reliable and efficient probabilistic a posteriori error estimator for the difference between a finite-dimensional linear operator and its orthogonal projection onto a reduced space is derived; the main idea is to apply the operator to standard Gaussian random vectors and consider the norm of the result. Also in <cit.>, an interpolation of the operator inverse is built via a Frobenius-norm projection and computed efficiently using randomized methods. An error estimator is obtained by measuring the norm of residual multiplied by the interpolation of the operator inverse, used here as a preconditioner.




We note that also the hierarchical error estimator for the RB method presented in <cit.> does not require the estimation of any stability constants, such as the inf-sup constant. In <cit.> the error is estimated by the distance between two reduced approximations of different accuracies and the computational costs depend highly on the dimension of the (primal) reduced space and are always higher than the costs for the computation of the RB approximation. In contrast, in our approach, the costs associated with the dual problems, and hence estimator evaluation, are commensurate with the cost associated with the (primal) RB approximation.
Finally, the reduced-order-model error surrogates (ROMES) method introduced in <cit.> and the closely related approaches <cit.> aim at constructing a statistical model for the approximation error. In <cit.> the statistical model is learned via stochastic-process data-fit methods from a small number of computed error indicators. 


The remainder of this article is organized as follows. In <ref> we derive a randomized a posteriori error estimator that estimates the error for a finite number of parameter values at given high probability. As this error estimator still depends on the high-dimensional solutions of dual problems, <ref> is devoted to the reduced order approximation of the dual problems and the analysis of the fast-to-evaluate a posteriori error estimator. In <ref> we demonstrate several theoretical aspects of the error estimator numerically and finally draw some conclusions in <ref>.






§ RANDOMIZED ERROR ESTIMATOR FOR PARAMETER-DEPENDENT EQUATION



 §.§ Parameter-dependent equations and error measurement


Consider a real-valued[Throughout the paper we consider real-valued equations: the extension of our method to the case of complex-valued problems is straightforward using the isomorphy ℂ=ℝ^2.] parameter-dependent equation 

    A(μ)u(μ)=f(μ),

where the parameter μ belongs to a parameter set 𝒫⊂ℝ^P. 
For every queried parameter μ∈𝒫, A(μ)∈ℝ^N× N is an invertible matrix and f(μ)∈ℝ^N. We assume we are given an approximation u(μ) of the solution u(μ). In this paper, the goal is to estimate the error

    u(μ)- u(μ)_Σ .

Here, ·_Σ is either a norm defined by means of a symmetric positive-definite (SPD) matrix Σ∈ℝ^N× N via v_Σ^2 = v^T Σ v for all v∈ℝ^N, or a semi-norm if Σ is only symmetric positive semi-definite. We highlight that the framework presented in this paper encompasses the estimation of the error in various different norms or the error in some QoI as will be discussed in the remainder of this subsection; see <ref> for a brief summary.

By choosing Σ = I_N, the identity matrix of size N, ·_Σ becomes the canonical norm ·_2 of ℝ^N. If problem <ref> stems from the discretization of a parameter-dependent linear partial differential equation, there is usually a natural norm ·_X associated with a Hilbert space of functions X⊂ H^1(D) for some spatial domain D⊂ℝ^d, d∈{1,2,3}. In such a case, there exists a discrete Riesz map R_X∈ℝ^N× N which is a SPD matrix such that √( (·)^T R_X (·) )= ·_X.
The choice Σ = R_X implies ·_Σ = ·_X, which means that the error is measured with respect to the natural norm of the problem. We may also consider for instance the error in the L^2-norm by choosing Σ=R_L^2(D), where the discrete Riesz map R_L^2(D) is chosen such that (·)^T R_L^2(D) (·) = ·_L^2(D)^2.


In some cases one is not interested in the solution u(μ) itself but rather in some QoI defined as a linear function of u(μ), say 

    s(μ)=L  u(μ)  ∈ℝ^m,

for some L∈ℝ^m × N.
In this situation one would like to estimate the error s(μ) - L u(μ) _W, where ·_W is a given natural norm on ℝ^m associated with a SPD matrix R_W so that w_W^2 = w^TR_W w for all w∈ℝ^m. With the choice Σ = L^T R_W L we can write

    u(μ) - u(μ) _Σ^2
        = (u(μ)-u(μ))^T ( L^T R_W L ) (u(μ)-u(μ)) =   s(μ) - Lu(μ) _W^2 ,

so that measuring the error with respect to the norm ·_Σ gives the error associated with the QoI. Notice that if m<N the matrix Σ is singular and ·_Σ is a semi-norm. Finally, consider the scalar-valued QoI given by s(μ)= l^T u(μ) where l∈ℝ^N. This corresponds to previous situation with m=1 and L=l^T. The choice Σ = l l^T yields u(μ) - u(μ) _Σ^2= |s(μ) - Lu(μ) |, where |·| denotes the absolute value. 









 §.§ Estimating norms using Gaussian maps




In this section we show how the (semi-)norm ·_Σ can be approximated by Φ·_2 for some random matrix Φ∈ℝ^K× N with K≪ N.

Let Z∼𝒩(0,Σ) be a zero mean Gaussian random vector in ℝ^N whose covariance matrix is chosen to be the matrix Σ∈ℝ^N× N which defines the (semi-)norm ·_Σ, cf <ref>.
Given a vector v∈ℝ^N, for example v=u(μ)-u(μ) for some (fixed) parameter μ∈𝒫, we can write

    v_Σ^2 
     = v^T Σ v 
     = v^T 𝔼( ZZ^T )v 
     = 𝔼( (Z^T v)^2 ),

where 𝔼(·) denotes the expected value.
This means that (Z^T v)^2 is an unbiased estimator of v_Σ^2. Let Z_1,,Z_K be K independent copies of Z and define the random matrix Φ∈ℝ^K× N whose i-th row is (1/√(K))Z_i^T. The matrix Φ is sometimes called a Gaussian map. Denoting by ·_2 the canonical norm of ℝ^K, we can write

    Φ v _2^2 = 1/K∑_i=1^K  ( Z_i^T v )^2     for any  v∈ℝ^N.

In other words, Φ v _2^2 is a K-sample Monte-Carlo estimator of 𝔼( (Z^T v )^2 ) = v_Σ^2. By the independence of the Z_i's, we have ( Φ v _2^2  ) = 1/K( Z^T v  ) so that Φ v _2^2 is a lower variance estimator of v_Σ^2 compared to (Z ^T v )^2. However, the variance is not always the most relevant criteria to assess the performance of an estimator. In the context of this paper, we rather want to quantify the probability that Φ v _2^2 deviates from v_Σ^2. This can be done by noting that, provided v_Σ≠0, the random variables (Z_i^T v ) / v_Σ for i=1,,K are independent standard normal random variables so that we have

    Φ v _2^2 = v_Σ^2/K∑_i=1^K ( Z_i^T v/v_Σ)^2 ∼v_Σ^2/K Q ,

where Q ∼χ^2(K) follows a chi-squared distribution with K degrees of freedom. Denoting by ℙ{A} the probability of an event A and by A the complementary event of A, the previous relation yields

    ℙ{ w^-1v_Σ≤Φ v _2 ≤ w v_Σ}   =1-ℙ{ Kw^-2≤  Q ≤ Kw^2 } ,

for any w≥1. Then for any given (fixed) vector v∈ℝ^N, the probability that a realization of Φ v _2 lies between w^-1v_Σ and wv_Σ is independent of v but also independent of the dimension N. The following proposition gives an upper bound for ℙ{ Kw^-2≤  Q ≤ Kw^2 } in terms of w and K. The proof, given in <ref>, relies on the fact that we have closed form expressions for the law of Q∼χ^2(K).



 Let Q∼χ^2(K) be a chi-squared random variable with K≥3 degrees of freedom. For any w>√(e) we have 
 
    ℙ{ Kw^-2≤  Q ≤ Kw^2 }≤( √(e)/w)^K .



<ref> shows that the probability ℙ{ Kw^-2≤  Q ≤ Kw^2 } decays at least exponentially with respect to K, provided w≥√(e) and K≥3. 
Then for any v∈ℝ^N, the relation

    w^-1v_Σ≤Φ v _2 ≤ w v_Σ ,

holds with a probability greater than 1-(√(e)/w)^K. As expected, a large value of w is beneficial to ensure the probability of failure (√(e)/w)^K to be small. For instance with w=4 and K=6, relation <ref> holds with a probability larger than 0.995. However, we observe in <ref> that this theoretical result is rather pessimistic since it overestimates the true probability by one order of magnitude for small values of w. Also, we conjecture on <ref> that there is an exponential decay even when w≤√(e) (see the blue curve with w=1.5), which is not predicted by <ref>.






In many situations we want to estimate the norm of several vectors rather than just one vector solely. This is for instance the case if one has to estimate the norm of the error v= u(μ) - u(μ) for many different parameter values μ∈𝒫. In that case, one would like to quantify the probability that relation <ref> holds simultaneously for any vector in a set ℳ⊂ℝ^N. Assuming ℳ is finite, a union bound argument
— for a detailed proof see <ref> — yields the following result:


Given a finite collection of vectors ℳ={v_1,v_2,, v_#ℳ}⊂ℝ^N and a failure probability 0<δ<1. Then, for any w > √(e) and 

    K≥min{log( #ℳ) + log(δ^-1)/log(w/√(e)) ,  3 }

we have

    ℙ{ w^-1v_Σ≤Φ v _2  ≤ w  v_Σ , ∀ v∈ℳ}≥ 1 - δ.





<ref> gives numerical values of K that satisfy <ref> depending on δ, w and #ℳ. For example with δ=10^-4 and w=10, estimating simultaneously the norm of 10^9 vectors requires only K=17 samples. Again, we emphasize that this result is independent on the dimension N of the vectors to be estimated.








 
 The Johnson-Lindenstrauss (JL) lemma states that for any 0<ε<1 and any finite set ℳ⊂ℝ^N, the condition K≥ 8ε^-2log(#ℳ) ensures the existence of a linear map Φ: ℝ^N→ℝ^K such that
 (1-ε) v-u_2^2 ≤Φ v -Φ u_2^2 ≤ (1+ε) v-u _2^2,
 holds for all u,v∈ℳ. 
 Replacing ℳ by ℳ∪{0} and letting u=0, one has that K≥ 8ε^-2log(#ℳ+1) is sufficient to ensure the existence of a Φ∈ℝ^K× N such that
 
    √(1-ε)v_2 ≤Φ v _2 ≤√(1+ε) v _2,   for all  v∈ℳ.

The above relation differs from <ref> in the sense that the deviation of Φ v _2 from v_2 is controlled in an additive manner via a parameter ε instead of a multiplicative way via w. We highlight also the different dependencies of K on ε and w. In contrast to the requirement in the JL lemma Condition <ref> permits reduction in the number of required copies K of the random vectors by considering an increased w. Note that the computational complexity of the a posteriori error estimator we propose in this paper crucially depends on K, see <ref>. Since the goal in this paper is to estimate the error we do in general not have to insist on a very accurate estimation of v_2. Instead, in many situtations it might be preferable to accept a higher effectivity w of the a posteriori error estimator in favour of a faster computational time. We emphasize that the user has the choice here. 
 
Notice also that with the choice w=1/√(1-ε), Equation <ref> implies <ref>. Then, the JL lemma ensures that <ref> holds true if K≥ 8(1-w^-2)^-2log(#ℳ+1). Even if we have the same logarithmic dependence on #ℳ, this is much larger than what we obtained in <ref>, already for moderate but especially for large values of w. For example with w=4, #ℳ=10^3 and δ=10^-2, JL lemma requires K≥ 63 whereas Condition <ref> requires only K≥ 13.

Finally, we highlight that a similar result to <ref> has been obtained in <cit.> for random vectors that are uniformly and randomly selected from the sphere. The multiplicative type of estimates in <cit.> motivated us to derive similar results for Gaussian vectors. 


 




In actual practice we can draw efficiently from Z∼𝒩(0,Σ) using a factorization of the covariance matrix of the form of Σ=U^TU, e.g. a (sparse) Cholesky decomposition.
It is then sufficient to draw a standard Gaussian vector Z and to compute the matrix-vector product Z=U^TZ.
As pointed-out in <cit.>, one can take advantage of a potential block structure of Σ to build a (non-square) factorization U with a negligible computational cost.






 §.§ Randomized a posteriori error estimator


We apply the methodology described in the previous subsection to derive a residual-based randomized a posteriori error estimator for the error u(μ)-u(μ)_Σ. Let Φ = K^-1/2[Z_1, Z_K]^T be a random matrix in ℝ^K× N where Z_1, Z_K are independent copies of Z∼𝒩(0,Σ), and consider the error estimator Δ(μ)=Φ( u(μ)-u(μ) ) _2, or equivalently

    Δ(μ)  
     = ( 1/K∑_k=1^K(Z_i^T(u(μ)-u(μ) ) )^2)^1/2.

If the parameter set 𝒫 is finite, <ref> with ℳ={ u(μ)-u(μ) ; μ∈𝒫} permits control of the quality of the estimate Δ(μ) uniformly over μ∈𝒫. But in actual practice the parameter set is often of infinite cardinality. Using more sophisticated techniques than just a simple union bound argument should provide results also when 𝒫 has infinite cardinality. In this paper, we are however only interested in the case of a finite set of parameter values, as restated in the following corollary.



 Let 0 < δ < 1 and w > √(e). Given a finite set of parameter values 𝒮⊂𝒫, the condition
 
    K≥min{log( #𝒮) + log(δ^-1)/log(w/√(e)) ,  3 },

 is sufficient to ensure
 
    ℙ{ w^-1Δ(μ) ≤u(μ)-u(μ) _Σ≤ w Δ(μ)  , ∀μ∈𝒮}≥ 1 - δ.



It is important to note that Condition <ref> depends only on the cardinality of 𝒮. This means that K can be determined only knowing the number of parameters for which we need to estimate the error.
However, computing Δ(μ) requires the solution u(μ) of problem <ref>, which is infeasible in practice. By introducing the residual

    r(μ)= f(μ) - A(μ)u(μ),

associated with Problem <ref> and, similar to <cit.>, exploiting the error residual relationship we may albeit rewrite the terms Z_i^T (u(μ)- u(μ)), 1≤ i ≤ K as follows:

    Z_i^T (u(μ) - u(μ)) = Z_i^T A(μ)^-1 r(μ) = (A(μ)^-TZ_i)^Tr(μ).

The terms Z_i^T (u(μ)- u(μ)) thus equal the inner products of the (primal) residual and the solutions Y_i(μ) ∈ℝ^N of the random dual problems 

    A(μ)^T Y_i(μ) = Z_i,    1≤ i ≤ K.



Because of the random right hand side in <ref>, the solutions Y_1(μ),,Y_K(μ) are random vectors. 
Thanks to the above the error estimator Δ(μ) <ref> can be rewritten as

    Δ(μ) 
     = ( 1/K∑_i=1^K ( Y_i(μ)^T r(μ) )^2 )^1/2 .

This shows that Δ(μ) can be computed by applying K linear forms to the residual r(μ).
In that sense, Δ(μ) can be considered as an a posteriori error estimator.
Notice that computing the solutions to <ref> is in general as expensive as solving the primal problem <ref>. In the next section we show how to approximate the dual solutions Y_1(μ),,Y_K(μ) in order to obtain a fast-to-evaluate a posteriori error estimator.


 When estimating the error in scalar-valued QoIs of the form of s(μ)=l^T u(μ), the covariance matrix is Σ=l l^T, see  <ref>. In that case the random vector Z∼𝒩(0,Σ) follows the same distribution as X l where X∼𝒩(0,1) is a standard normal random variable (scalar). The random dual problem <ref> then becomes
 A(μ)^TY_i(μ) = X_i  l
 and the solution is Y_i(μ) = X_i   q(μ) where q(μ) is the solution of the deterministic dual problem A(μ)^T q(μ) = l. Dual problems of this form are commonly encountered for estimating linear quantities of interest, see <cit.> for a general presentation and <cit.> for the application in reduced order modeling.
 




Say that we use the a posteriori error estimator Δ(μ) to select a new parameter and use the associated solution to enrich the reduced space. Then, we wish to use Δ(μ) again for the enriched reduced space. However, now the problem occurs that the error between u(μ) and the reduced solution that uses the newly selected snapshot depends on the error estimator and thus Z_1,,Z_K; we lose independence. One solution would be to redraw the samples in each iteration, which is in general however computationally infeasible. Alternatively, as suggested in <cit.>, we can adapt the number of samples K in order to take into account (using union bound arguments) all possible outcomes of the greedy algorithm; for further details we refer to <cit.>.





§ A FAST-TO-EVALUATE RANDOMIZED A POSTERIORI ERROR ESTIMATOR


In order to obtain a fast-to-evaluate a posteriori error estimator whose computational complexity is independent of N, we employ projection-based model order reduction (MOR) techniques to compute approximations of the solutions Y_1(μ),,Y_K(μ) of the dual problems <ref>. To that end, let us assume that we are given a fixed realization of the K random vectors Z_1,,Z_K and that we have a reduced space 𝒴⊂ℝ^N at our disposal.
Different ways to construct 𝒴 will be discussed in <ref> and compared numerically in <ref>. In any case, 𝒴 will be built from dual solutions Y_i(μ), i=1,,K of (<ref>) for random right-hand sides
Z_1,,Z_K, the latter being fixed before constructing the dual reduced space. 𝒴 should thus be considered as a random subspace. 
Then, we define Y_i(μ) as the Galerkin projection of Y_i(μ) on 𝒴, meaning

    Y_i(μ) ∈𝒴 :   ⟨  A(μ)^TY_i(μ) ,v ⟩ = ⟨  Z_i ,v ⟩ ,   ∀ v∈𝒴.

Here, ⟨ v,w⟩:=v^Tw for all v,w ∈ℝ^N. 
We emphasize that we employ the same reduced space 𝒴 for the approximation of the K  dual solutions Y_1(μ),,Y_K(μ). Needless to say that a segregated strategy, where we construct and use K different dual reduced spaces 𝒴_i for the K different right-hand sides Z_i and associated dual reduced solutions Y_i(μ), i=1,,K, can also be considered. The advantage of a segregated strategy is that one can easily parallelize the computations, if needed. However, in this paper we focus exclusively on the monolithic approach <ref>, employing one single dual reduced space. 


By replacing Y_i(μ) in <ref> by the fast-to-evaluate approximation Y_i(μ), we define a fast-to-evaluate a posteriori error estimator as

    Δ(μ) := ( 1/K∑_i=1^K (Y_i(μ)^T r(μ))^2)^1/2.

We highlight that, in constrast to for instance the “standard” a posteriori error estimator being defined as the product of the reciprocal of a stability constant and the dual norm of the primal residual, Δ(μ) does not contain any constants that require estimation. Moreover, unlike hierarchical error estimators <cit.> the quality of the approximation used for the error estimator does not depend on the quality of the primal approximation; the dual reduced space does not in general contain the primal reduced space as a subspace and can even be of smaller dimension than the latter. For a more elaborate comparison we refer to <ref>.

Additionally, we shall show in <ref> that evaluating μ↦Δ(μ) requires only the solution of one linear system of size n_𝒴:=(𝒴), instead of K linear systems of size n_𝒴 as suggested by <ref>. 
However, before discussing the computational complexity of Δ(μ), we show in <ref> that under certain conditions Δ(μ) is both a reliable and efficient error estimator at high probability. Based on this analysis we propose in <ref> different greedy algorithms for constructing the reduced space 𝒴. 










 §.§ Analysis of the fast-to-evaluate a posteriori error estimator


First, we relate the relative error in the a posteriori error estimator to the error in the dual residual:





 Assume Σ is invertible. The fast-to-evaluate error estimator Δ(μ) defined by <ref> satisfies
 
    | Δ(μ) - Δ(μ) |/u(μ)-u(μ)_Σ   ≤max_1≤ i≤ K A^T(μ) Y_i(μ) - Z_i _Σ^-1     for all μ∈𝒫.

Here, ·_Σ^-1 denotes the norm on ℝ^N such that v_Σ^-1^2 = v^TΣ^-1v for all v∈ℝ^N.
 


The proof is given in <ref>. Notice that <ref> requires Σ to be invertible, which excludes the cases where one wants to estimate the error in a vector-valued QoI, see <ref>. <ref> allows us to control the error via Δ(μ), where the effectivity w is enlarged in an additive manner, as stated in the following corollary.



 Suppose we are given a finite set of parameter values 𝒮⊂𝒫 for which we want to estimate the error u(μ)-u(μ) _Σ.
 Let 0 < δ < 1, w > √(e) and assume
 
    K≥min{log( #𝒮) + log(δ^-1)/log(w/√(e)) ,  3 }.

 Furthermore, assume that Σ is invertible and that we have ε≤ w^-1, where
 
    ε = sup_μ∈𝒫{max_1≤ i≤ K A^T(μ) Y_i(μ) - Z_i _Σ^-1} .

 Then, we have
 
    ℙ{ 
      (w + ε)^-1Δ(μ) 
      ≤u(μ)-u(μ) _Σ≤w/1 - w  ε Δ(μ) ,
        ∀μ∈𝒮}≥ 1-δ .






The proof is given in <ref>.
<ref> gives a sufficient condition to control the quality of the estimator Δ(μ) over a finite set of parameter values 𝒮⊂𝒫 with high probability. It requires ε≤ w^-1, which is equivalent to A^T(μ) Y_i(μ) - Z_i _Σ^-1≤ w^-1 for all μ∈𝒫 and all 1≤ i ≤ K. 
To satisfy this condition, one has to design an algorithm which builds Y_i(μ) in a way that A^T(μ) Y_i(μ) is close to Z_i uniformly over 𝒮 and independently on the value taken by Z_i. Obtaining ε≤ w^-1 can however be challenging (from a computational perspective). To explain this, let us note that Z_i_Σ^-1 is, with high probability[To show this, note that Z_i_Σ^-1^2∼χ^2(N) so that, by <ref>, relation w'^-1√(N)≤Z_i_Σ^-1≤ w'√(N) holds with probability 1-(√(e)/w')^N for any w'≥√(e).], of the order of √(N).
Therefore ε≤ w^-1 means that the relative dual residual norm ought to be of the order of

    A^T(μ) Y_i(μ) - Z_i _Σ^-1/ Z_i _Σ^-1≃1/w√(N).

When N≫1, the condition ε≤ w^-1 means that we need a very accurate approximation of the dual variables. For instance with N=10^6, the dual residual norm A^T(μ) Y_i(μ) - Z_i _Σ^-1 /  Z_i _Σ^-1 has to be less that 10^-3 for all μ∈𝒫 and all 1≤ i≤ K, which can be too demanding in actual practice.

Next, we give an alternative way of controlling the quality of Δ(μ).
Contrarily to <ref>, which provides a additive type of control, the following proposition gives a control in an multiplicative manner. The proof in given in <ref> in the appendix.


 Suppose we are given a finite set of parameter values 𝒮⊂𝒫 over which we want to estimate the error u(μ)-u(μ) _Σ.
 Let 0 < δ < 1, w > √(e) and assume
 
    K≥min{log( #𝒮) + log(δ^-1)/log(w/√(e)) ,  3 }.

 Then the fast-to-evaluate estimator Δ(μ) satisfies
 
    ℙ{ 
      (α w)^-1Δ(μ) 
      ≤u(μ)-u(μ) _Σ≤ (α w)  Δ(μ),
        μ∈𝒮,
      }≥ 1-δ ,

 where
 
    α :=  max_μ∈𝒫( max{Δ(μ)/Δ(μ) , Δ(μ)/Δ(μ)}) ≥ 1.



<ref> shows that, with high probability, the error estimator Δ(μ) departs from the true error u(μ)-u(μ) _Σ at most by a multiplicative factor (α w)^-1 or (α w). 

Notice that α is a measure of the distance from μ↦Δ(μ) to μ↦Δ(μ): if it is close to 1 then Δ(μ) is close to Δ(μ) uniformly over the parameter set 𝒫. Unlike <ref>, <ref> does not require Σ to be invertible and, even more importantly, it does not put any restrictions on α. However, the computation of α can be expensive since it requires the exact error estimator Δ(μ) over the whole parameter set 𝒫. Therefore, we propose to use α as a stopping criterion when constructing the dual reduced space to ensure that (α w)^-1Δ(μ) 
  ≤u(μ)-u(μ) _Σ≤ (α w)  Δ(μ) holds true for a rich training set ⊂𝒫 as we will detail in <ref>.









 §.§ Greedy constructions of the dual reduced space 𝒴





  §.§.§ Vector point of view of the dual problems


A popular technique to build a reduced space is to take the span of snapshots of the solution. In order to handle the K distinct dual problems, the index “i” in <ref> in considered as an additional parameter. Thus, we define the augmented parameter set 𝒫_K = { 1,,K }×𝒫 and seek a n_𝒴-dimensional reduced space of the form of

    𝒴 = span{ Y_i_1(μ_1),,Y_i_n_𝒴(μ_n_𝒴) } ,

where the n_𝒴 elements (i_1,μ_1),,(i_n_𝒴,μ_n_𝒴) are to be chosen in 𝒫_K. The RB methodology (see for instance <cit.> for an introduction) consists in selecting (i_1,μ_1),,(i_n_𝒴,μ_n_𝒴) in a greedy fashion <cit.>. In detail, assuming that the j first parameters are given, the (j+1)-th parameter is defined as

    (i_j+1,μ_j+1) ∈ (i,μ)∈𝒫_K^trainargmax A(μ)^T Y_i(μ) - Z_i _* ,

where Y_i(μ) is the approximation of Y_i(μ) given by <ref> with 𝒴 defined as in <ref>. 
Here, 𝒫_K^train⊂𝒫_K is a sufficiently rich training set with finite cardinality and ·_* denotes an arbitrary norm of ℝ^N.

According to <ref>, it is natural to chose ·_* = ·_Σ^-1, provided Σ is invertible. After having computed the snapshot Y_i_j+1(μ_j+1), the reduced space 𝒴 is updated using <ref> with j← j+1. By selecting the parameter (i_j+1,μ_j+1) according to <ref>, the idea is to construct a reduced space that minimizes the dual residual norm A(μ)^T Y_i(μ) - Z_i _* uniformly over the training set (i,μ)∈𝒫_K^train. 


It remains to define a criterion to stop the greedy iterations.
Given a user-defined tolerance ≥0, the use of the stopping criterion

    max_(i,μ)∈𝒫_K^train A(μ)^T Y_i(μ) - Z_i _* ≤ ,

ensures that, at the end of the iteration procedure, the residual norm of the dual problem is below  everywhere on the training set 𝒫_K^train.
One can relax that criterion by replacing the max in <ref> by the quantile of order q ∈[0,1]:

    { A(μ)^T Y_i(μ) - Z_i _*  :   (i,μ)∈𝒫_K^train}≤ .

Here {A} denotes the ⌈ q#A ⌉-th largest entries of a (ordered and finite) set A.
With this stopping criterion, the iterations stop when at least a fraction of q points in 𝒫_K^train have a dual residual norm below . 
Notice that the  and the max coincides when q=1 so that <ref> generalizes <ref>.
The resulting greedy algorithm is summarized in <ref>.






  §.§.§ Matrix point of view of the dual problems


Next, we propose another greedy algorithm which relies on a matrix interpretation of the K dual problems <ref>. Let us denote by 

    𝐘(μ) =[ Y_1(μ),, Y_K(μ)] ∈ℝ^N× K ,

the matrix containing the dual solutions. Instead of constructing the reduced space 𝒴 as the span of vectors Y_i(μ), like in Equation <ref>, we now consider reduced spaces of the form of

    𝒴 = span{𝐘(μ_1) λ_1 , , 𝐘(μ_n_𝒴) λ_n_𝒴},

where λ_1,,λ_n_𝒴 are n_𝒴 vectors in ℝ^K and where μ_1,,μ_n_𝒴∈𝒫.
Notice that if the vectors λ_i are canonical vectors of ℝ^K we have that 𝒴 can be written as in <ref>.
In that sense, the approximation format <ref> is richer than <ref> and we can expect better performance. We also note that the greedy algorithm we propose here shares some similarities with the POD-greedy algorithm introduced in <cit.>.

We now propose a second greedy algorithm inspired by <ref>. Let 𝒫^train⊂𝒫 be again a finite training set and suppose that at step r in the greedy algorithm we have a reduced space 𝒴 as in <ref> at our disposal. The first step is to define the next evaluation point μ_j+1 as

    μ_j+1∈μ∈𝒫^trainargmax( max{Δ(μ)/Δ(μ) , Δ(μ)/Δ(μ)}) ,

where we recall that Δ(μ)=(1/K∑_i=1^K[Z_i^T(u(μ)-u(μ))]^2)^1/2. 
Finding μ_j+1 according to <ref> requires to compute the solution u(μ) over the training set μ∈𝒫^train. As this is in general not computationally efficient, we suggest replacing u(μ) by a reference solution u_ref(μ) such that u_ref(μ) - u(μ)_Σ≪u(μ) - u(μ)_Σ. We can choose as u_ref(μ) for instance a hierarchical approximation of u(μ), where we use a larger primal reduced space to determine u_ref(μ). Note that we only suggest using such a reference solution for the construction of the dual reduced space and not afterwards when certifying the reduced approximation in the online stage.
Then, we introduce the reference error estimator 

    Δ_ref(μ):=(1/K∑_i=1^K( Z_i^T(u_ref(μ)-u(μ)))^2)^1/2

and seek μ_j+1 as

    μ_r+1∈μ∈𝒫^trainargmax( max{Δ_ref(μ)/Δ(μ) , Δ(μ)/Δ_ref(μ)}).

Once the parameter μ_j+1 is found either with <ref> of with <ref>, we compute the dual solutions Y_1(μ_j+1),,Y_K(μ_j+1) and assemble 𝐘(μ_j+1). Here we need to solve K linear equations with the same operator A(μ_j+1)^T but with K different right-hand sides, see Equation <ref>. This can be done efficiently say be using a Cholesky or LU decomposition and reusing the factorization for the K problems.

The second step is to determine the vector λ_j+1. In order to maximize the improvement of the reduced space, we propose to define λ_j+1 as follows:

    λ_j+1∈λ∈ℝ^Kargmax 𝐘(μ_j+1) λ - 𝐘(μ_j+1) λ_2 /λ_2  ,

where 𝐘(μ_j+1) = [ Y_1(μ_j+1),, Y_K(μ_j+1)]. The rational behind <ref> is to align λ_j+1 with the direction where the matrix 𝐘(μ_j+1) differs the most from 𝐘(μ_j+1). One can easily show that λ_j+1 defined by <ref> is the first eigenvector of the K-by-K matrix 

    M(μ_j+1) = (  𝐘(μ_j+1) - 𝐘(μ_j+1)  )^T( 𝐘(μ_j+1) - 𝐘(μ_j+1) ) .

Once λ_j+1 is computed, we set j← j+1 and we update the reduced space 𝒴 using <ref>. We terminate the algorithm based on the following stopping criteria

    {max{Δ(μ)/Δ(μ) , Δ(μ)/Δ(μ)} : μ∈𝒫^train}≤.

The resulting greedy algorithm is summarized in <ref>.





Note that in the POD-greedy algorithm <cit.> one would consider the orthogonal projection on the reduced space 𝒴 instead of the actual reduced solutions in <ref>. However, for problems where the Galerkin projection deviates significantly from the orthogonal projection, we would expect that using the reduced solution gives superior results than the POD-greedy as the latter does not take into account the error due to the Galerkin projection which can be significant for instance close to resonances in a Helmholtz problem. We have performed numerical experiments for the same benchmark problem (parametrized Helmholtz equation) we consider in <ref> that confirm this conjecture.






 §.§ Computational aspects of the fast-to-evaluate error estimator



At a first glance the complexity for evaluating μ↦Δ(μ) is dominated by the solution of the K reduced problems <ref>, meaning K times the solution of a (dense) linear system of equations of size n_𝒴. The next proposition, inspired by Lemma 2.7 in <cit.>, shows that one can actually evaluate μ↦Δ(μ) by solving only one linear system of size n_𝒴, which reduces the previous complexity by a factor K; the proof is provided in <ref>. 
Note however that the complexity for evaluating μ↦Δ(μ) is not completely independent on K. Indeed, as we employ the same reduced space for the approximation of K dual problems, the dimension of n_𝒴 depends on K. The rate of the increase of n_𝒴 for growing K will be investigated in numerical experiments in <ref>.


 The error indicator Δ(μ) defined by (<ref>) can be written as
 
    Δ(μ) 
      = ( 1/K∑_i=1^K  (Z_i^T  e(μ) )^2 )^1/2 ,

 where e(μ) ∈𝒴 is the solution to
 
    e(μ)∈𝒴 ,   ⟨ A(μ) e(μ),v ⟩ 
      = ⟨ r(μ), v ⟩ ,   ∀ v∈𝒴.

 


Besides giving an alternative way of computing Δ(μ), <ref> also gives a new insight into the fast-to-evaluate error estimator. Reformulating Problem <ref> as

    e(μ)∈𝒴 ,   ⟨ A(μ) (u(μ) + e(μ) ),v ⟩ = ⟨ f(μ), v ⟩ ,   ∀ v∈𝒴 ,

demonstrates that e(μ)∈𝒴 may be interpreted as a correction of the primal approximation u(μ), so that u(μ) + e(μ) is an enriched solution of the original problem <ref> compared to u(μ). Since 𝒴 is not designed for improving the primal approximation u(μ), one cannot reasonably hope that the correction e(μ) improves significantly u(μ). 
However the norm of e(μ), estimated by the fast-to-evaluate error estimator <ref>, gives relevant information about the error u(μ)-u(μ)_Σ. Finally, we emphasize again that the primal reduced space is in general not a subspace of the dual reduced space and that the intersection of the primal and dual reduced space can even be empty. As a consequence, the right-hand side in (<ref>) is in general not zero for test functions from the dual reduced space.



 Assume A(μ) = ∑_q=1^Q_Aα_q(μ) A_q and f(μ) = ∑_q=1^Q_fζ_q(μ) f_q with A_q, f_q parameter-independent and consider u(μ) as the Galerkin projection onto some primal reduced order space 𝒳 of dimension n_𝒳.
Since all inner products involving high dimensional quantities can be preassembled, the marginal computational complexity of Δ(μ) is 𝒪(Q_An_𝒴^2 + Q_fn_𝒴 + Q_An_𝒴n_𝒳) for assembling <ref>, 𝒪(n_𝒴^3) for solving <ref> and 𝒪(K n_𝒴) for calculating <ref>. 
 For moderate Q_A the marginal computational complexity of Δ(μ) is thus dominated by 𝒪(n_𝒴^3), i.e. the costs for solving <ref>.




 An alternative strategy for estimating the error is to measure the distance between the approximation u(μ) ∈𝒳 and a reference solution u_ref(μ), which is an improved approximation of u(μ) compared to u(μ). 
 When using projection based model order reduction u_ref(μ) can be defined as a Galerkin projection onto an enriched reduced space of the form of 𝒳+𝒴, as proposed in <cit.>. Unlike our approach, the space 𝒴 ought to be adapted for capturing the error u(μ) - u(μ). The complexity for evaluating such a hierarchical error estimator is dominated by the solution of a dense system of equations of size dim(𝒳+𝒴). 
 In contrast, our approach requires the solution of a system of equations whose size is independent on the dimension of the primal reduced space 𝒳, see the above <ref>.








§ NUMERICAL EXPERIMENTS



We numerically demonstrate various theoretical aspects of the proposed error estimator.
Our benchmark is a parameterized Helmholtz equation for which a reduced order solution is obtained by the RB method. 
Estimating the error in this reduced order model is challenging because, around the resonances, we lose the coercivity of the operator which makes a posteriori error estimation quite difficult with standard methods.

Let us mention here that all the training sets 𝒫^train (or 𝒫^train_K) and all the online sets 𝒮 are comprised of snapshots selected independently and uniformly at random in 𝒫 (or in 𝒫_K).
Those (random) sets are redrawn at each new simulation, unless mentioned otherwise.






 §.§ Benchmark: Multi-parametric Helmholtz equation





Consider the parameterized Helmholtz equation

- ∂_x_1x_1- μ_1 ∂_x_2x_2 - μ_2     =         in  D:=(0,1)×(0,1),

   =0         on  (0,1)×{0},

         

∂_x_2    = cos(πx_1)             on  (0,1)×{1}, 

∂_x_1    = 0         on  {0,1}×(0,1).

The solution = (μ) is parameterized by μ=(μ_1,μ_2) ∈𝒫:=[0.2,1.2]× [10,50], where μ_1 accounts for anisotropy and μ_2 is the wavenumber squared.
The source term  is defined by (x_1,x_2)= _1(x_1)_2(x_2) for any (x_1,x_2)∈ D, where

_1(x_1) :=  
5       if  0 ≤x_1 ≤0.1,

-5     if  0.2 ≤x_1 ≤0.3,

10     if 0.45 ≤x_1 ≤0.55,

-5     if  0.7 ≤x_1 ≤8,

5     if  0.9 ≤x_1 ≤1,

0     else,
    and  
_2(x_2) := 
1       if  0.5 ≤x_2 ≤1, 

0     else.


A similar test case with a smaller parameter set has been considered in <cit.>.
The resonances can be determined analytically and are depicted by the black lines in <ref>. Because of the multi-parameter setting, we have resonance surfaces which are more difficult to deal with than a union of isolated resonance frequencies in the single-parameter setting; see <cit.>. Moreover, we observe that in the region [0.2,0.4] × [30,50] ⊂𝒫 there are quite a few resonance surfaces that are also relatively close together, making this an even more challenging situation both for the construction of suitable reduced models and even more for a posteriori error estimation.  

We employ the Finite Element (FE) method to discretize the weak solution of <ref>.
To that end, we define a FE space X^h⊂ X:={𝔳∈ H^1(D)  :  𝔳(x_1,0) = 0} by means of a regular mesh with square elements of edge length h=0.01 and FE basis functions that are piecewise linear in x_1 and x_2 direction, resulting in a FE space of N = (X^h)=10100. The FE approximation ^h(μ) is defined as the Galerkin projection of (μ) on X^h, and we denote by u(μ)∈ℝ^N the vector containing the coefficients of ^h(μ) when expressing it in the FE basis. Moreover, we denote by R_X∈ℝ^N× N the discrete Riesz map associated with the H^1-norm, which is such that u(μ) R_X u(μ) = ^h(μ) _H^1(D)^2 for any μ∈𝒫. By default the covariance matrix Σ is always chosen to be Σ=R_X, unless mentioned otherwise.


We may also consider a QoI defined as the trace of the FE solution on the boundary Γ = {0}× (0,1) ⊂∂ D, meaning ^h_|Γ(μ).
We denote by s(μ) ∈ℝ^100 the vector containing those entries of u(μ) ∈ℝ^N that are associated with the grid points on Γ. Then, we can write s(μ)=Lu(μ) where L∈ℝ^100× N is an extraction matrix. To measure the error associated with the QoI, we use the norm ·_W defined as the discretization of the L^2(Γ) norm, which is such that s(μ)_W = _|Γ^h(μ) _L^2(Γ) for any μ∈𝒫.

The primal RB approximation u(μ) is defined as the Galerkin projection of u(μ) onto the space of snapshots, meaning
u(μ) ∈𝒳 := span{ u(μ_1),u(μ_2),} ,
where the parameters μ_1,μ_2, are selected in a greedy way based on the dual norm of the residual associated with <ref>. Each time we run <ref>, we use a reference solution u_ref(μ) defined as an RB approximation of u(μ) using n_𝒳 + 10 basis functions, where n_𝒳:=(𝒳).
Note that this reference solution appears only in the offline stage.



 §.§ Randomized a posteriori error estimation with exact dual





We demonstrate here the statistical properties of the error estimator Δ(μ) defined by <ref>.

<ref> shows histograms of the effectivity indices {Δ(μ)/u(μ) - u(μ) _Σ , μ∈𝒮} for five different realizations of the vectors Z_1,,Z_K. Here, the same online set 𝒮 with #𝒮=10^4 is used.

We observe that for each of the five realizations, the effectivity indices Δ(μ)/u(μ) - u(μ) _Σ lie in the interval [1/w , w] for any μ∈𝒮, as predicted by <ref>.
This theoretical bound looks however pessimistic, as the effectivities for K=5 (resp. K=10) lie in the interval [1/w,w] that corresponds to K=10 (resp. K=20).
This might be due to the rather crude union bound argument.

The solid lines on <ref> represent the probability density function (pdf) of √(Q/K) where Q∼χ^2(K).
This is the pdf of Δ(μ)/u(μ) - u(μ) _Σ for any fixed μ.

Even though the histograms depicted on <ref> are not representing that pdf (instead they represent the distribution of the effectivity index among the set 𝒮), we observe good accordance with the black line. In particular we observe a concentration phenomenon of the histograms around 1 when K increases.








 §.§ Approximation of the dual problems



  §.§.§ Construction of the dual space


In <ref> we compare the maximum, the minimum, the 95% quantile and the 99% quantile of {Δ(μ)/u(μ)-u(μ)_Σ : μ∈𝒮} where the dual reduced space is constructed either by <ref> (with ·_*=·_R_X^-1), by <ref> or by a POD.

We observe that by using <ref> we need many fewer dual basis functions than for <ref> and for the POD. 

In detail, we see for instance in <ref>, <ref> and <ref> that for K=5 employing <ref> requires about n_𝒴=20 dual basis functions to have 99% of the samples in the interval [1/3,3], while when using the <ref> or the POD we need about 35 or 30 basis functions, respectively. We emphasize that for K=20 the difference is even larger. 
Moreover, when considering the QoI (last row of <ref>), the difference between <ref> and the <ref> and POD is less pronounced but still considerable. This significant disparity can be explained by the fact that while both the POD and the <ref> try to approximate the K dual solutions Y_1(μ),,Y_K(μ), <ref> is driven by the approximation of the error estimator Δ(μ) and thus a scalar quantity; compare the selection criteria <ref> and <ref>. 
This also explains why the discrepancy increases significantly for growing K: While POD and <ref> have to approximate a more complex object (the K dual solutions), we only obtain an additional summand in Δ(μ) for each additional random right-hand side. 
Let us also highlight the significant difference between the maximum value and the 99% quantile over the parameter set and the somewhat erratic behavior of the maximum, which both seem to be due to the resonance surfaces. As indicated above this motivates considering for instance the 99% quantile as a stopping criterion in both <ref> and <ref>. 







  §.§.§ Dimension of the dual space



<ref> shows statistics of the dimension of the dual reduced space 𝒴 obtained by <ref> with different stopping criterion. We consider =0.5 and relax this tolerance by multiplication with a varying relaxation parameter ρ taking the values 1,10 and 100. We observe that except for ·ρ=50 and moderate K the dimension of the dual space is in general quite large. Comparing with <ref>, we observe that choosing ·ρ=50 is albeit sufficient to obtain an effectivity close to 1.

Notice however that the use of <ref> requires ·ρ≈ε≤ 1/w ≤ 1, which excludes ·ρ=5 and ·ρ=50. <ref> shows the evolution of the stopping criteria during the first 80 iterations of <ref>. 

We observe a significant impact of K on the convergence profiles: with K=20 the curves do not attain the tolerance =0.5, which explains the results we observed in <ref>.





















In comparison, <ref> yields much smaller dual reduced spaces; compare <ref> and <ref>. We see in <ref> that, except for =1.5 the dimension of the dual RB space 𝒴 is smaller than the dimension of the primal RB space 𝒳 when using the 95%, 97.5%, 99%-quantiles for the stopping criterion. Moreover, for instance for =3 we see that for n_𝒳=20,30 we can use (significantly) less dual than primal basis functions. However, we also see that tight tolerances for  will lead in general to dual reduced spaces that have a larger dimension than the primal RB space. 
As larger tolerances ≥ 5 may lead to an significant underestimation of the error (see <ref>), tolerances for  between 1 and 4 seem to be preferable. <ref> shows the evolution of the stopping criteria during the first 80 iterations of <ref>. Note that for higher tolerances for  it may happen for a realization that <ref> terminates in a valley between two peaks.

Furthermore, we observe in <ref> a very large standard deviation of about 10 if we consider the maximum over the offline training set, while for the 95%,97.5%,99% quantiles we have often a standard deviation of about 2. Additionally, the dimension of the dual reduced spaces for the maximum is much larger than for the considered quantiles, but among the considered quantiles we observe only very moderate changes. Again, it seems that this behavior is due to the resonance surfaces. Moreover, as we obtain a very satisfactory effectivity of Δ(μ) when we use for instance the 99% quantile (see <ref>), we conclude that using quantiles between 97.5% and 99% as a stopping criterion in <ref> seems advisable.

Finally, we observe both a very moderate dependency of the dimension of the dual reduced space constructed by <ref> on K and a rather mild dependency on n_𝒳. Therefore, we conjecture that the proposed error estimator might also be applied rather complex problems. 






  §.§.§ Performance of Δ(μ) on an online parameter set



On <ref> we plot the histograms of {Δ(μ)/u(μ)-u(μ)_Σ : μ∈𝒮} for 5 realizations of the random vectors Z_1,,Z_K, where the dual reduced space 𝒴 is built via <ref>. 
We observe a similar behavior as for the error estimator Δ(μ) with the exact dual, see <ref>. In particular for all μ∈𝒮 the effectivity index Δ(μ)/u(μ)-u(μ)_Σ lies  between (α w)^-1 and (α w), see <ref>, where α is estimated by =2.
Finally, we highlight that <ref> demonstrates that with near certainty we obtain effectivities near unity with a dual space dimension on the same order as (or less than) the primal space dimension. Hence the costs for the a posteriori error estimator are about the same as those for constructing the primal approximation.




In order to understand the average performance of the online-efficient error indicator, we plot in <ref> the histograms of the concatenation of 100 realizations of the effectivity indices {Δ(μ)/u(μ)-u(μ)_Σ : μ∈𝒮}. Here, for each new realization, we redraw the K vectors Z_1,,Z_K, the training set 𝒫^train, then run <ref> to construct the dual reduced space 𝒴, and finally redraw the online set 𝒮.

We observe that for a larger tolerance  the histograms are shifted to the left, which seems to be a bit stronger for larger K (corresponding to smaller w).
This is due to the fact that <ref> is stopped earlier and the dimension of 𝒴 is not sufficiently large to approximate well the error estimator Δ(μ). 
Nevertheless, we observe that the effectivity indices Δ(μ)/u(μ)-u(μ)_Σ are always in the interval [(α w)^-1 , (α w)], where α≈, as expected thanks to <ref>.
This shows that, even with a rather crude approximation of the dual solutions, it is safe to use the fast-to-evaluate error estimator Δ(μ), as the grey area is taking into account the approximation error in the error estimator.


To guarantee that the effectivity indices lie in a user-defined interval of the form of [c^-1,c], it is sufficient to choose α and w such that α w=c, see <ref>. As a consequence there is a degree of freedom in the choice of α and w, meaning in the choice of ≈α and K=K(w) via relation <ref>. To avoid a too large shift of the histogram to the left as for instance observed for w=2.1 and =3 it seems advisable to choose α at least as small as w. Additionally, the plots corresponding to w =3.2,2.1 and α=2 highlight the importance of choosing α small enough compared to w if one is interested in rather tight estimates. However, decreasing α has, as anticipated, a much stronger effect on the dimension of the dual reduced space (see <ref>). Therefore, it seems that for the considered test case choosing α/w ∈ (1/3,1) seems to be a good compromise between computational costs and effectivity of the error estimator. We also see that for instance w=6.5 and α=3 or α=2 yield already very good results in this direction.











§ CONCLUSIONS

In this paper we introduced a randomized a posteriori error estimator for low-rank approximations, which is constant-free and is both reliable and efficient at given high probability. Here, the upper and lower bound of the effectivity is chosen by the user. To derive the error estimator we exploit the concentration phenomenon of Gaussian maps. 
Exploiting the error residual relationship and approximating the associated random dual problems via projection-based model order reduction yields a fast-to-evaluate a posteriori error estimator. 
We highlight that we had to put some effort in proving the concentration inequalities but regarding the parametrized problem we only relied on its well-posedness and the definition of the adjoint operator. Therefore, there is some chance that the presented framework might be extended quite easily to more complex problems.

To construct the dual reduced space we employed a greedy algorithm guided by a quantity of interest that assesses the quality of the fast-to-evaluate error estimator. The numerical experiments for a multi-parametric Helmholtz problem show that we obtain much smaller dual reduced spaces than with a standard greedy driven by the dual norm of the residual or with the POD. Moreover, the numerical experiments demonstrate that for moderate upper bounds for the effectivities of about 20 the dimension of the dual reduced space needs only to be a bit more than half of the dimension of the primal reduced space. If a very tight effectivity bound of about 2 or 3 is desired the dual reduced spaces have to be about twice as large as the primal approximation spaces. We emphasize however that even for larger bounds of the effectivity thanks to the concentration of measure  the effectivity is still very often close to one. Furthermore, we observed only a very moderate dependence of the dimension of the dual reduced space on the number of random vectors K, which controls the variance of the estimator and a very mild dependence on the dimension on the (primal) reduced space. This might indicate that the error estimator will also perform well for challenging problems. Finally, we showed that to compute the fast-to-evaluate a posteriori error estimator we need to solve one dense linear system of equations of the size of the dimension of the dual reduced space. 

Due to the above the proposed a posteriori error estimator features a very favorable computational complexity and its computational costs are often about the same as the costs for the low-rank approximation or even smaller for moderate effectivity bounds. The presented error estimator can thus be more advantageous from a computational viewpoint than error estimators based on the dual norm of the residual and a (costly to estimate) stability constant or hierarchical type error estimators.









§ PROOFS



 §.§ Proof of <ref>


First we give a bound for ℙ{ Q≤ Kw^-2}. This quantity corresponds to the cumulative distribution function of the χ^2(K) distribution evaluated at Kw^-2. We have ℙ{ Q≤ Kw^-2} = 1/Γ(K/2)γ( K/2 , K/2w^2), where Γ(·) is the gamma function such that Γ(a)=∫_0^∞ t^a-1 e^-tdt and γ(·,·) the lower incomplete gamma function defined by γ(a,x)=∫_0^x t^a-1 e^-tdt. Following the lines of <cit.>, we can write γ(a,x)≤∫_0^x t^a-1dt = 1/ax^a and 

    Γ(a) = ∫_0^a t^a-1 e^-tdt + ∫_a^∞ t^a-1 e^-tdt 
     ≥ e^-a∫_0^a t^a-1dt + a^a-1∫_a^∞  e^-tdt = 2 a^a-1 e^-a ,

whenever a≥1. Then, if K≥2 we have

    ℙ{ Q≤ Kw^-2} 
     = 1/Γ(K/2) γ( K/2 , K/2w^2) 
     ≤(K/2)e^K/2/2(K/2)^K/2·(K/(2w^2))^K/2/K/2
     = 1/2(√(e)/w)^K .

Now we give a bound for ℙ{  Q ≥ Kw^2 }.
Using a Markov inequality, for any 0≤ t < 1/2 we can write

    ℙ{  Q ≥ Kw^2 }   = ℙ{ e^ t Q ≥ e^t Kw^2}≤𝔼( e^ t Q   ) /e^ t Kw^2 
     =  (1-2t)^-K/2/e^t Kw^2  ,

where for the last equality we used the expression for the moment-generating function of χ^2(K). The minimum of the above quantity is attained for t=(w^2-1)/(2w^2) so we can write

    ℙ{  Q ≥ Kw^2 }   ≤ (w^2 e^ 1-w^2  )^K/2 
     = ( √(e)/w)^K(w^2 e^ -w^2/2  )^K≤( √(e)/w)^K2^K/e^K≤1/2( √(e)/w)^K,

for any K≥ 3. Together with <ref>, the previous inequalities allows writing

    ℙ{ Kw^-2≤  Q ≤ Kw^2 } 
     = ℙ{ Q≤ Kw^-2} + ℙ{  Q ≥ Kw^2 }≤( √(e)/w)^K ,

for any K≥3, which concludes the proof.





 §.§ Proof of <ref>

A union bound allows writing

    ℙ{ w^-1v_Σ≤Φ v _2  ≤ w  v_Σ , ∀ v∈ℳ}
         ≥ 1-∑_v∈ℳℙ{w^-1v_Σ≤Φ v _2 ≤ w v_Σ}
         = 1-(#ℳ)  ℙ{ Kw^-2≤  Q ≤ Kw^2 }≥ 1-(#ℳ)  ( √(e)/w)^K,

where, for the last inequality, we used Proposition <ref> (assuming w > √(e) and K≥3 hold).
Given 0<δ<1, condition 
K≥log( #ℳ) + log(δ^-1)/log(w/√(e)) ,
is equivalent to 1-(#ℳ)( √(e)/w )^K≥ 1-δ and ensures that (<ref>) holds for all v∈ℳ with probability larger than 1-δ. 



 §.§ Proof of <ref>


Let Ψ(μ)=K^-1/2[Y_1(μ),,Y_K(μ)]^T and Ψ(μ)=K^-1/2[Y_1(μ),,Y_K(μ)]^T so that, from Equations <ref> and <ref>, we can write Δ(μ) = Ψ(μ) r(μ)_2 and Δ(μ) = Ψ(μ) r(μ)_2.  Using a triangle inequality we can write

    | Δ(μ) - Δ(μ) | 
      = | Ψ(μ) r(μ) _2 - Ψ(μ) r(μ) _2 | 
      ≤Ψ(μ) r(μ) - Ψ(μ) r(μ) _2

Dividing by u(μ)-u(μ)_Σ we can write

    | Δ(μ) - Δ(μ) | /u(μ)-u(μ)_Σ   ≤ ( Ψ(μ) - Ψ(μ) ) r(μ) _2 /u(μ)-u(μ)_Σ =  ( Ψ(μ)- Ψ(μ)) A(μ) ( u(μ)-u(μ) ) _2 /u(μ)-u(μ)_Σ
       ≤sup_v∈ℝ^N \{0} ( Ψ(μ)- Ψ(μ)) A(μ) v _2 /v_Σ
       = sup_v_Σ=1√(1/K∑_i=1^K ( ( A(μ)^T Y_i(μ) - Z_i )^T v )^2 )
       ≤sup_v_Σ=1 max_1≤ i ≤ K | ( A(μ)^T Y_i(μ) - Z_i )^T v |  
       = max_1≤ i ≤ K ( A(μ)^T Y_i(μ) - Z_i )^T v _Σ^-1 ,

which yields <ref> and concludes the proof.





 §.§ Proof of <ref>


By <ref> we have | Δ(μ) - Δ(μ) | ≤εu(μ)-u(μ)_Σ, which is equivalent to

    Δ(μ) - εu(μ)-u(μ)_Σ≤Δ(μ) ≤Δ(μ) + εu(μ)-u(μ)_Σ.

By <ref>, it holds with probability larger than 1-δ that w^-1Δ(μ) ≤u(μ)-u(μ) _Σ≤ w Δ(μ) for all μ∈𝒮. Then with the same probability we have

    (w^-1 - ε ) u(μ)-u(μ)_Σ≤Δ(μ) ≤ (w + ε )u(μ)-u(μ)_Σ,

for all μ∈𝒮, which yields <ref> and concludes the proof.






 §.§ Proof of <ref>

By <ref>, it holds with probability larger than 1-δ that w^-1Δ(μ) ≤u(μ)-u(μ) _Σ≤ w Δ(μ) for all μ∈𝒮. Then with the same probability we have

    u(μ)-u(μ) ≤  w Δ(μ) 
     ≤  w ( sup_μ'∈𝒮Δ(μ')/Δ(μ')) Δ(μ) 
     <ref>≤ (α w) Δ(μ),

and

    u(μ)-u(μ) ≥ w^-1Δ(μ) 
     ≥ w^-1( inf_μ'∈𝒮Δ(μ')/Δ(μ')) Δ(μ) 
     <ref>≥ (α w)^-1Δ(μ),

for any μ∈𝒮, which yields <ref> and concludes the proof.




 §.§ Proof of <ref>


By construction, both Y_i(μ) and e(μ) belong to 𝒴. Then for all i=1,,K we can write

    Y_i(μ)^T r(μ)
        =⟨ r(μ), Y_i(μ) ⟩<ref>=⟨ A(μ)e(μ),Y_i(μ) ⟩
        = ⟨e(μ), A(μ)^T Y_i(μ) ⟩<ref>=⟨e(μ), Z_i ⟩  
     = Z_i^T e(μ).

Then, by definition <ref> we can write

    Δ(μ) 
     = ( 1/K∑_k=1^K(Y_i(μ)^T r(μ)  )^2)^1/2 
     = ( 1/K∑_k=1^K(Z_i^T e(μ)  )^2)^1/2

which gives the result.





siam


