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

HAL
open science

# Projection based model order reduction methods for the estimation of vector-valued variables of interest

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy

► To cite this version:

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy. Projection based model order reduction methods for the estimation of vector-valued variables of interest. SIAM Journal on Scientific Computing, 2017, 39 (4), pp.A1647-A1674. 10.1137/16M106385X. hal-01387236

HAL Id: hal-01387236

https://hal.science/hal-01387236

Submitted on 5 Jun 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_1.png" /></figure>

HAL
open science

# Projection based model order reduction methods for the estimation of vector-valued variables of interest

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy

► To cite this version:

Olivier Zahm, Marie Billaud-Friess, Anthony Nouy. Projection based model order reduction methods for the estimation of vector-valued variables of interest. SIAM Journal on Scientific Computing, 2017, 39 (4), pp.A1647-A1674. 10.1137/16M106385X. hal-01387236

HAL Id: hal-01387236

https://hal.science/hal-01387236v1

Submitted on 5 Jun 2024

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
# Projection based model order reduction methods for the estimation of vector-valued variables of interest*

Olivier ZAHM§, Marie BILLAUD-FRIESS†, Anthony NOUY†‡

November 1, 2016

## Abstract

We propose and compare goal-oriented projection based model order reduction methods for the estimation of vector-valued functionals of the solution of parameter-dependent equations. The first projection method is a generalization of the classical primal-dual method to the case of vector-valued variables of interest. We highlight the role played by three reduced spaces: the approximation space and the test space associated to the primal variable, and the approximation space associated to the dual variable. Then we propose a Petrov-Galerkin projection method based on a saddle point problem involving an approximation space for the primal variable and an approximation space for an auxiliary variable. A goal-oriented choice of the latter space, defined as the sum of two spaces, allows us to improve the approximation of the variable of interest compared to a primal-dual method using the same reduced spaces. Then, for both approaches, we derive computable error estimates for the approximations of the variable of interest and we propose greedy algorithms for the goal-oriented construction of reduced spaces. The performance of the algorithms are illustrated on numerical examples and compared to standard (non goal-oriented) algorithms.

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author (anthony.nouy@ec-nantes.fr).</footnote>

<footnote>*This work was supported by the French National Research Agency (Grant ANR CHORUS MONU-0005)</footnote>

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author anthony.nouy@ec-nantes.fr.</footnote>

<footnote>§Department of Aeronautics & Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.</footnote>
# Projection based model order reduction methods for the estimation of vector-valued variables of interest*

Olivier ZAHM§, Marie BILLAUD-FRIESS†, Anthony NOUY†‡

November 1, 2016

## Abstract

We propose and compare goal-oriented projection based model order reduction methods for the estimation of vector-valued functionals of the solution of parameter-dependent equations. The first projection method is a generalization of the classical primal-dual method to the case of vector-valued variables of interest. We highlight the role played by three reduced spaces: the approximation space and the test space associated to the primal variable, and the approximation space associated to the dual variable. Then we propose a Petrov-Galerkin projection method based on a saddle point problem involving an approximation space for the primal variable and an approximation space for an auxiliary variable. A goal-oriented choice of the latter space, defined as the sum of two spaces, allows us to improve the approximation of the variable of interest compared to a primal-dual method using the same reduced spaces. Then, for both approaches, we derive computable error estimates for the approximations of the variable of interest and we propose greedy algorithms for the goal-oriented construction of reduced spaces. The performance of the algorithms are illustrated on numerical examples and compared to standard (non goal-oriented) algorithms.

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author (anthony.nouy@ec-nantes.fr).</footnote>

<footnote>*This work was supported by the French National Research Agency (Grant ANR CHORUS MONU-0005)</footnote>

<footnote>†Ecole Centrale de Nantes, GeM, UMR CNRS 6183, France.</footnote>

<footnote>‡Corresponding author anthony.nouy@ec-nantes.fr.</footnote>

<footnote>§Department of Aeronautics & Astronautics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA.</footnote>
# 1 Introduction

This paper is concerned with the numerical solution of linear equations of the form

$$A(\xi)u(\xi) = b(\xi), \qquad (1)$$

where the operator $A(\xi)$ and right-hand side $b(\xi)$ depend on a parameter $\xi$ which takes values in some parameter set $\Xi$. Such equations arise in many contexts such as uncertainty quantification, optimization or control, where the solution of (1) have to be evaluated with many instances of the parameters (multi-query context). For large systems of equations (e.g. arising from a fine discretization of a parameter-dependent partial differential equation), solving (1) for one instance of the parameter can be very expensive, which leads to intractable computations in a multi-query context. Model order reduction methods aim at constructing an approximation of the solution map $u : \Xi \to V$ whose evaluation for a certain value of $\xi$ is cheaper than solving (1). Standard approaches rely on Galerkin-type projections of $u(\xi)$ on a low-dimensional subspace $V_r$ of the solution space $V$, a so-called reduced space. The reduced space can be generated from evaluations (snapshots) of the solution $u(\xi)$ at some selected (or randomly chosen) values of the parameter $\xi$, see [5, 14, 17, 19]. The Proper Orthogonal Decomposition method aims at constructing an optimal subspace for the approximation of the set of solutions $M = \{u(\xi) : \xi \in \Xi\}$ in a mean-square sense (see [14]). Reduced Basis (RB) methods (see [11] for a survey) aim at controlling the approximation uniformly over the parameter set. In this context, reduced spaces are usually constructed using greedy algorithms.

In many applications one is not interested in the solution $u(\xi)$ itself, but only in a variable of interest $s(\xi)$ which is a functional of $u(\xi)$. Here we assume that $s(\xi)$ depends linearly on $u(\xi)$. Efficient goal-oriented methods have been proposed for the estimation of a scalar-valued variable of interest $s(\xi)$. A standard method consists in computing an approximation of the solution of the so-called dual problem associated to (1) which is used to correct the estimation of $s(\xi)$. We refer to [16] for a general survey on primal-dual methods and to [6, 10, 11, 17] for the application in the context of RB methods.

In this paper, we propose projection based model order reduction methods for the estimation of a variable of interest $s(\xi)$ taking values in a vector space of finite or infinite dimension. We consider the case where

$$s(\xi) = L(\xi)u(\xi),$$

with $L(\xi)$ a parameter-dependent linear operator. For example, for boundary value problems, $L(\xi)$ can be defined as the trace operator providing the restriction of the solution
# 1 Introduction

This paper is concerned with the numerical solution of linear equations of the form

$$A(\xi)u(\xi) = b(\xi), \quad (1)$$

where the operator $A(\xi)$ and right-hand side $b(\xi)$ depend on a parameter $\xi$ which takes values in some parameter set $\Xi$. Such equations arise in many contexts such as uncertainty quantification, optimization or control, where the solution of (1) have to be evaluated with many instances of the parameters (multi-query context). For large systems of equations (e.g. arising from a fine discretization of a parameter-dependent partial differential equation), solving (1) for one instance of the parameter can be very expensive, which leads to intractable computations in a multi-query context. Model order reduction methods aim at constructing an approximation of the solution map $u: \Xi \to V$ whose evaluation for a certain value of $\xi$ is cheaper than solving (1). Standard approaches rely on Galerkin-type projections of $u(\xi)$ on a low-dimensional subspace $V_r$ of the solution space $V$, a so-called reduced space. The reduced space can be generated from evaluations (snapshots) of the solution $u(\xi)$ at some selected (or randomly chosen) values of the parameter $\xi$, see [5, 14, 17, 19]. The Proper Orthogonal Decomposition method aims at constructing an optimal subspace for the approximation of the set of solutions $M = \{u(\xi) : \xi \in \Xi\}$ in a mean-square sense (see [14]). Reduced Basis (RB) methods (see [11] for a survey) aim at controlling the approximation uniformly over the parameter set. In this context, reduced spaces are usually constructed using greedy algorithms.

In many applications one is not interested in the solution $u(\xi)$ itself, but only in a variable of interest $s(\xi)$ which is a functional of $u(\xi)$. Here we assume that $s(\xi)$ depends linearly on $u(\xi)$. Efficient goal-oriented methods have been proposed for the estimation of a scalar-valued variable of interest $s(\xi)$. A standard method consists in computing an approximation of the solution of the so-called dual problem associated to (1) which is used to correct the estimation of $s(\xi)$. We refer to [16] for a general survey on primal-dual methods and to [6, 10, 11, 17] for the application in the context of RB methods.

In this paper, we propose projection based model order reduction methods for the estimation of a variable of interest $s(\xi)$ taking values in a vector space of finite or infinite dimension. We consider the case where

$$s(\xi) = L(\xi)u(\xi),$$

with $L(\xi)$ a parameter-dependent linear operator. For example, for boundary value problems, $L(\xi)$ can be defined as the trace operator providing the restriction of the solution
to the boundary of the domain. In this case the variable of interest belongs to an infinite dimensional space or, after discretization, to a finite but possibly high dimensional space. The standard approach, which consists in treating $s(\xi)$ as a collection of scalar-valued variables of interest and in building one reduced dual space for each of them, has a complexity which grows proportionally to the dimension of $s(\xi)$. Our approach circumvents this issue by constructing a single reduced dual space, thus allowing to handle variables of interest with high and potentially infinite dimension. A similar approach can be found for parametric dynamical systems, see the monograph [3] for a general introduction. In this framework, projection-based model order reduction methods are used for the approximation of $s(\xi)$ which is an output of the dynamical system. Petrov-Galerkin methods have been proposed with different ways of constructing the reduced basis for the test and trial space, such as the balanced truncation methods, (balanced) Proper Orthogonal Decomposition method, moment matching methods, etc. We refer to [4] for a recent review on these methods. In the present paper, we aim at exploring other possibilities than the Petrov-Galerkin projection.

In a first part, we introduce and analyze different methods for computing projections of the solution and approximations of the variable of interest. We first present a non goal-oriented Petrov-Galerkin approach to compute an approximation of $u(\xi)$ from which an estimation of $s(\xi)$ is deduced. Then, we introduce a generalization of the standard primal-dual method to the case of a vector-valued variable of interest, which relies on the approximation of the primal variable $u(\xi)$ and of the solution $Q(\xi)$ of the dual problem

$$A(\xi)^*Q(\xi) = L(\xi)^*,$$

where $A(\xi)^*$ and $L(\xi)^*$ are the adjoints of operators $A(\xi)$ and $L(\xi)$ respectively. We show that the error on the variable of interest depends on three reduced spaces: the approximation space $V_r$ for the primal variable $u(\xi)$, the test space $W_r$ which is used for the Petrov-Galerkin projection of $u(\xi)$, and an approximation space $W_k^Q$ for the dual variable $Q(\xi)$ which is projected on the space of $W_k^Q$-valued linear operators. Finally, we present a Petrov-Galerkin method where the projection is obtained by solving a saddle point problem which involves an approximation space $V_r$ for $u(\xi)$ and an approximation space $T_p$ for an auxiliary variable. We show that if $T_p$ is defined by $T_p = W_r + W_k^Q$, then error bounds for both the projection of the primal variable on $V_r$ and the approximation of the variable of interest can be improved compared to error bounds of a primal-dual approach using the same spaces $V_r, W_r$ and $W_k^Q$. The proposed approach is a goal-oriented extension of the method proposed in [8].

In a second part, we derive (for both approaches) computable error estimates for the approximation of the variable of interest. Then, we propose greedy algorithms based on these
to the boundary of the domain. In this case the variable of interest belongs to an infinite dimensional space or, after discretization, to a finite but possibly high dimensional space. The standard approach, which consists in treating $s(\xi)$ as a collection of scalar-valued variables of interest and in building one reduced dual space for each of them, has a complexity which grows proportionally to the dimension of $s(\xi)$. Our approach circumvents this issue by constructing a single reduced dual space, thus allowing to handle variables of interest with high and potentially infinite dimension. A similar approach can be found for parametric dynamical systems, see the monograph [3] for a general introduction. In this framework, projection-based model order reduction methods are used for the approximation of $s(\xi)$ which is an output of the dynamical system. Petrov-Galerkin methods have been proposed with different ways of constructing the reduced basis for the test and trial space, such as the balanced truncation methods, (balanced) Proper Orthogonal Decomposition method, moment matching methods, etc. We refer to [4] for a recent review on these methods. In the present paper, we aim at exploring other possibilities than the Petrov-Galerkin projection.

In a first part, we introduce and analyze different methods for computing projections of the solution and approximations of the variable of interest. We first present a non goal-oriented Petrov-Galerkin approach to compute an approximation of $u(\xi)$ from which an estimation of $s(\xi)$ is deduced. Then, we introduce a generalization of the standard primal-dual method to the case of a vector-valued variable of interest, which relies on the approximation of the primal variable $u(\xi)$ and of the solution $Q(\xi)$ of the dual problem

$$A(\xi)^*Q(\xi) = L(\xi)^*,$$

where $A(\xi)^*$ and $L(\xi)^*$ are the adjoints of operators $A(\xi)$ and $L(\xi)$ respectively. We show that the error on the variable of interest depends on three reduced spaces: the approximation space $V_r$ for the primal variable $u(\xi)$, the test space $W_r$ which is used for the Petrov-Galerkin projection of $u(\xi)$, and an approximation space $W_k^Q$ for the dual variable $Q(\xi)$ which is projected on the space of $W_k^Q$-valued linear operators. Finally, we present a Petrov-Galerkin method where the projection is obtained by solving a saddle point problem which involves an approximation space $V_r$ for $u(\xi)$ and an approximation space $T_p$ for an auxiliary variable. We show that if $T_p$ is defined by $T_p = W_r + W_k^Q$, then error bounds for both the projection of the primal variable on $V_r$ and the approximation of the variable of interest can be improved compared to error bounds of a primal-dual approach using the same spaces $V_r, W_r$ and $W_k^Q$. The proposed approach is a goal-oriented extension of the method proposed in [8].

In a second part, we derive (for both approaches) computable error estimates for the approximation of the variable of interest. Then, we propose greedy algorithms based on these
error estimates for the construction of the reduced spaces $V_r$ and $W_k^Q$. We discuss different choices for the reduced space $W_r$. In particular, we introduce a parameter-dependent space depending on a preconditioner obtained by means of an interpolation of the inverse of the operator $A(\xi)$ proposed in [22].

This paper is organized as follows. In Section 2, we introduce and analyze the different projection methods for the estimation of vector-valued variables of interest for general linear equations of the form (1) formulated in a Hilbert setting. Then, in Section 3, we derive error estimates for the approximation of the variable of interest and we propose practical greedy algorithms for the construction of reduced spaces. Finally, in Section 4, numerical experiments illustrate the properties of the projection methods and of the greedy algorithms. In particular, we provide a simplified complexity analysis for the so-called offline phase (i.e. the construction of the reduced spaces) and for the online phase (i.e. the evaluation of $s(\xi)$ for a particular instance of $\xi$).

## 2 Projection methods for the estimation of a variable of interest

Let $V, W$ and $Z$ be three Hilbert spaces. For a Hilbert space $H$ equipped with a norm $\|\cdot\|_H$, we denote by $H'$ the topological dual space of $H$. We consider the linear equation

$$Au = b \qquad (2)$$

with $A \in \mathcal{L}(V, W')$ and $b \in W'$, and a variable of interest

$$s = Lu,$$

where $L \in \mathcal{L}(V, Z)$. We assume that $A$ is a norm-isomorphism<fnref n="1" /> such that for all $u \in V$,

$$\alpha \|u\|_V \le \|Au\|_{W'} \le \beta \|u\|_V,$$

where

$$\inf_{0 \neq v \in V} \sup_{0 \neq w \in W} \frac{\langle Av, w \rangle}{\|v\|_V \|w\|_W} := \alpha > 0 \qquad (3a)$$

$$\sup_{0 \neq v \in V} \sup_{0 \neq w \in W} \frac{\langle Av, w \rangle}{\|v\|_V \|w\|_W} := \beta < \infty, \qquad (3b)$$

<footnote n="1">¹A is a norm-isomorphism if it is a continuous and weakly coercive operator satisfying the assumptions of the Nečas' theorem [9, Chapter 2].</footnote>
error estimates for the construction of the reduced spaces $V_r$ and $W_k^Q$. We discuss different choices for the reduced space $W_r$. In particular, we introduce a parameter-dependent space depending on a preconditioner obtained by means of an interpolation of the inverse of the operator $A(\xi)$ proposed in [22].

This paper is organized as follows. In Section 2, we introduce and analyze the different projection methods for the estimation of vector-valued variables of interest for general linear equations of the form (1) formulated in a Hilbert setting. Then, in Section 3, we derive error estimates for the approximation of the variable of interest and we propose practical greedy algorithms for the construction of reduced spaces. Finally, in Section 4, numerical experiments illustrate the properties of the projection methods and of the greedy algorithms. In particular, we provide a simplified complexity analysis for the so-called offline phase (i.e. the construction of the reduced spaces) and for the online phase (i.e. the evaluation of $s(\xi)$ for a particular instance of $\xi$).

## 2 Projection methods for the estimation of a variable of interest

Let $V$, $W$ and $Z$ be three Hilbert spaces. For a Hilbert space $H$ equipped with a norm $\|\cdot\|_H$, we denote by $H'$ the topological dual space of $H$. We consider the linear equation

$$Au = b \qquad (2)$$

with $A \in \mathcal{L}(V, W')$ and $b \in W'$, and a variable of interest

$$s = Lu,$$

where $L \in \mathcal{L}(V, Z)$. We assume that $A$ is a norm-isomorphism<fnref n="1" /> such that for all $u \in V$,

$$\alpha \|u\|_V \le \|Au\|_{W'} \le \beta \|u\|_V,$$

where

$$\inf_{0 \neq v \in V} \sup_{0 \neq w \in W} \frac{\langle Av, w \rangle}{\|v\|_V \|w\|_W} := \alpha > 0 \qquad (3a)$$

$$\sup_{0 \neq v \in V} \sup_{0 \neq w \in W} \frac{\langle Av, w \rangle}{\|v\|_V \|w\|_W} := \beta < \infty, \qquad (3b)$$

<footnote n="1">¹A is a norm-isomorphism if it is a continuous and weakly coercive operator satisfying the assumptions of the Nečas' theorem [9, Chapter 2].</footnote>
which ensures the well-posedness of (2). In this section, we present different methods for constructing an approximation $\tilde{s}$ of $s$. First, in Section 2.1, we present a standard approach which consists in estimating the variable of interest from a Petrov-Galerkin projection of $u$. In Section 2.2, we present an extension of the primal-dual approach to the case of vector-valued variables of interest, where the variable of interest is estimated from a standard Petrov-Galerkin projection of the primal variable and a projection of the solution of a dual problem. Finally, in Section 2.3, we introduce a goal-oriented projection method based on a saddle-point formulation.

Before going further, let us introduce some additional notations. For a Hilbert space $H$, we denote by $R_H \in \mathcal{L}(H, H')$ the Riesz map such that $\|v\|_H^2 = \langle R_H v, v \rangle$, where $\langle \cdot, \cdot \rangle$ denotes the duality pairing. The dual norm $\| \cdot \|_{H'}$ on $H'$ is such that $R_{H'} = R_H^{-1}$. Then $\|v\|_H = \|R_H v\|_{H'}$ and $|\langle v, w \rangle| \le \|v\|_H \|w\|_{H'}$ hold for any $v \in H$ and $w \in H'$. For any operator $C \in \mathcal{L}(H_1, H_2')$, with $H_1$ and $H_2$ two Hilbert spaces, $C^* \in \mathcal{L}(H_2, H_1')$ denotes the adjoint of $C$, such that $\langle C v_1, v_2 \rangle = \langle v_1, C^* v_2 \rangle$ for any $v_1 \in H_1$ and $v_2 \in H_2$.

## 2.1 Petrov-Galerkin projection

Suppose that we are given a subspace $V_r \subset V$ of finite dimension $r$ in which we seek an approximation of $u$. The orthogonal projection $u_r^\perp$ of $u$ on $V_r$, given by $\|u - u_r^\perp\|_V = \min_{v \in V_r} \|u - v\|_V$, is characterized by

$$ \langle u - u_r^\perp, R_V v \rangle = 0, \quad \forall v \in V_r. \qquad (4) $$

In practice, an approximation $u_r \in V_r$ can be defined as a Petrov-Galerkin projection of $u$ characterized by

$$ \langle Au_r - b, y \rangle = 0, \quad \forall y \in W_r, \qquad (5) $$

where $W_r \subset W$ is a test space of dimension $r$. Under the assumption that

$$ \alpha_{V_r, W_r} = \inf_{0 \neq v \in V_r} \sup_{0 \neq y \in W_r} \frac{\langle Av, y \rangle}{\|v\|_V \|y\|_W} > 0, \qquad (6) $$

the next proposition provides a quasi-optimality result for $u_r$ and gives an error bound for the approximation of the variable of interest. In what follows, notation min (resp. max) is used in place of inf (resp. sup) when the minimum (resp. the maximum) is reached.

**Proposition 2.1.** Under assumption (6), the solution $u_r$ of equation (5) satisfies

$$ \|u - u_r\|_V \le \frac{1}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_V. \qquad (7) $$
which ensures the well-posedness of (2). In this section, we present different methods for constructing an approximation $\tilde{s}$ of $s$. First, in Section 2.1, we present a standard approach which consists in estimating the variable of interest from a Petrov-Galerkin projection of $u$. In Section 2.2, we present an extension of the primal-dual approach to the case of vector-valued variables of interest, where the variable of interest is estimated from a standard Petrov-Galerkin projection of the primal variable and a projection of the solution of a dual problem. Finally, in Section 2.3, we introduce a goal-oriented projection method based on a saddle-point formulation.

Before going further, let us introduce some additional notations. For a Hilbert space $H$, we denote by $R_H \in \mathcal{L}(H, H')$ the Riesz map such that $\|v\|_H^2 = \langle R_H v, v \rangle$, where $\langle \cdot, \cdot \rangle$ denotes the duality pairing. The dual norm $\| \cdot \|_{H'}$ on $H'$ is such that $R_{H'} = R_H^{-1}$. Then $\|v\|_H = \|R_H v\|_{H'}$ and $|\langle v, w \rangle| \le \|v\|_H \|w\|_{H'}$ hold for any $v \in H$ and $w \in H'$. For any operator $C \in \mathcal{L}(H_1, H_2')$, with $H_1$ and $H_2$ two Hilbert spaces, $C^* \in \mathcal{L}(H_2, H_1')$ denotes the adjoint of $C$, such that $\langle C v_1, v_2 \rangle = \langle v_1, C^* v_2 \rangle$ for any $v_1 \in H_1$ and $v_2 \in H_2$.

## 2.1 Petrov-Galerkin projection

Suppose that we are given a subspace $V_r \subset V$ of finite dimension $r$ in which we seek an approximation of $u$. The orthogonal projection $u_r^\perp$ of $u$ on $V_r$, given by $\|u - u_r^\perp\|_V = \min_{v \in V_r} \|u - v\|_V$, is characterized by

$$ \langle u - u_r^\perp, R_V v \rangle = 0, \quad \forall v \in V_r. \quad (4) $$

In practice, an approximation $u_r \in V_r$ can be defined as a Petrov-Galerkin projection of $u$ characterized by

$$ \langle Au_r - b, y \rangle = 0, \quad \forall y \in W_r, \quad (5) $$

where $W_r \subset W$ is a test space of dimension $r$. Under the assumption that

$$ \alpha_{V_r, W_r} = \inf_{0 \neq v \in V_r} \sup_{0 \neq y \in W_r} \frac{\langle Av, y \rangle}{\|v\|_V \|y\|_W} > 0, \quad (6) $$

the next proposition provides a quasi-optimality result for $u_r$ and gives an error bound for the approximation of the variable of interest. In what follows, notation min (resp. max) is used in place of inf (resp. sup) when the minimum (resp. the maximum) is reached.

**Proposition 2.1.** Under assumption (6), the solution $u_r$ of equation (5) satisfies

$$ \|u - u_r\|_V \le \frac{1}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_V. \quad (7) $$
where

$$ \delta_{V_r, W_r} = \max_{0 \neq v \in V_r} \min_{y \in W_r} \frac{\|v - R_V^{-1} A^* y\|_V}{\|v\|_V} < 1. \qquad (8) $$

Furthermore,

$$ \|s - Lu_r\|_Z \le \frac{\delta_{W_r}^L}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_V, \qquad (9) $$

with

$$ \delta_{W_r}^L = \sup_{0 \neq z' \in Z'} \min_{y \in W_r} \frac{\|L^*z' - A^*y\|_{V'}}{\|z'\|_{Z'}}. \qquad (10) $$

Proof. With $u_r^\perp$ the orthogonal projection of $u$ on $V_r$, for any $v \in V_r \setminus \{0\}$ and $y \in W_r$, we have

$$ \begin{aligned} \langle u_r^\perp - u_r, R_V v \rangle &\stackrel{(4)}{=} \langle u - u_r, R_V v \rangle = \langle b - Au_r, A^{-*}R_V v \rangle \\ &\stackrel{(5)}{=} \langle b - Au_r, A^{-*}R_V v - y \rangle = \langle u - u_r, R_V v - A^*y \rangle \\ &\le \|u - u_r\|_V \|R_V v - A^*y\|_{V'}. \end{aligned} $$

Taking the minimum over $y \in W_r$, dividing by $\|v\|_V$ and taking the maximum over $v \in V_r \setminus \{0\}$, we obtain $\|u_r^\perp - u_r\|_V \le \delta_{V_r, W_r} \|u - u_r\|_V$, where $\delta_{V_r, W_r}$ is defined by (8). Thanks to the orthogonality condition (4) we have $\|u - u_r\|_V^2 = \|u - u_r^\perp\|_V^2 + \|u_r^\perp - u_r\|_V^2$, from which we deduce that $(1 - \delta_{V_r, W_r}^2) \|u - u_r\|_V^2 \le \|u - u_r^\perp\|_V^2$. To prove (7), it remains to prove that $\delta_{V_r, W_r} < 1$. Noting that

$$ \min_{y \in W_r} \|v - R_V^{-1} A^* y\|_V^2 = \min_{0 \neq y \in W_r} \min_{\lambda \in \mathbb{R}} \|v - \lambda R_V^{-1} A^* y\|_V^2 = \min_{0 \neq y \in W_r} \|v\|_V^2 - \frac{\langle v, A^* y \rangle^2}{\|A^* y\|_{V'}^2}, $$

we obtain

$$ \delta_{V_r, W_r}^2 = 1 - \min_{0 \neq v \in V_r} \max_{0 \neq y \in W_r} \frac{\langle Av, y \rangle^2}{\|v\|_V^2 \|A^* y\|_{V'}^2}. \qquad (11) $$

Let introduce $\beta_{W_r} = \sup_{0 \neq y_r \in W_r} \|A^*y\|_{V'} / \|y\|_W$ which, from assumption (3b), satisfies $\beta_{W_r} \le \beta$. Then using assumption (6) we obtain

$$ \delta_{V_r, W_r}^2 \le 1 - \frac{\alpha_{V_r, W_r}^2}{\beta_{W_r}^2} \le 1 - \frac{\alpha_{V_r, W_r}^2}{\beta^2} < 1. \qquad (12) $$

Furthermore for any $z' \in Z' \setminus \{0\}$ and $y \in W_r$, we have

$$ \begin{aligned} \langle s - Lu_r, z' \rangle &= \langle b - Au_r, A^{-*}L^*z' \rangle \stackrel{(5)}{=} \langle b - Au_r, A^{-*}L^*z' - y \rangle \\ &\le \|u - u_r\|_V \|L^*z' - A^*y\|_{V'}. \end{aligned} $$
where

$$ \delta_{V_r, W_r} = \max_{0 \neq v \in V_r} \min_{y \in W_r} \frac{\|v - R_V^{-1} A^* y\|_V}{\|v\|_V} < 1. \qquad (8) $$

Furthermore,

$$ \|s - Lu_r\|_Z \le \frac{\delta_{W_r}^L}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_V, \qquad (9) $$

with

$$ \delta_{W_r}^L = \sup_{0 \neq z' \in Z'} \min_{y \in W_r} \frac{\|L^*z' - A^*y\|_{V'}}{\|z'\|_{Z'}}. \qquad (10) $$

Proof. With $u_r^\perp$ the orthogonal projection of $u$ on $V_r$, for any $v \in V_r \setminus \{0\}$ and $y \in W_r$, we have

$$ \begin{aligned} \langle u_r^\perp - u_r, R_V v \rangle &\stackrel{(4)}{=} \langle u - u_r, R_V v \rangle = \langle b - Au_r, A^{-*}R_V v \rangle \\ &\stackrel{(5)}{=} \langle b - Au_r, A^{-*}R_V v - y \rangle = \langle u - u_r, R_V v - A^*y \rangle \\ &\le \|u - u_r\|_V \|R_V v - A^*y\|_{V'}. \end{aligned} $$

Taking the minimum over $y \in W_r$, dividing by $\|v\|_V$ and taking the maximum over $v \in V_r \setminus \{0\}$, we obtain $\|u_r^\perp - u_r\|_V \le \delta_{V_r, W_r} \|u - u_r\|_V$, where $\delta_{V_r, W_r}$ is defined by (8). Thanks to the orthogonality condition (4) we have $\|u - u_r\|_V^2 = \|u - u_r^\perp\|_V^2 + \|u_r^\perp - u_r\|_V^2$, from which we deduce that $(1 - \delta_{V_r, W_r}^2) \|u - u_r\|_V^2 \le \|u - u_r^\perp\|_V^2$. To prove (7), it remains to prove that $\delta_{V_r, W_r} < 1$. Noting that

$$ \min_{y \in W_r} \|v - R_V^{-1} A^* y\|_V^2 = \min_{0 \neq y \in W_r} \min_{\lambda \in \mathbb{R}} \|v - \lambda R_V^{-1} A^* y\|_V^2 = \min_{0 \neq y \in W_r} \|v\|_V^2 - \frac{\langle v, A^* y \rangle^2}{\|A^* y\|_{V'}^2}, $$

we obtain

$$ \delta_{V_r, W_r}^2 = 1 - \min_{0 \neq v \in V_r} \max_{0 \neq y \in W_r} \frac{\langle Av, y \rangle^2}{\|v\|_V^2 \|A^* y\|_{V'}^2}. \qquad (11) $$

Let introduce $\beta_{W_r} = \sup_{0 \neq y_r \in W_r} \|A^*y\|_{V'} / \|y\|_W$ which, from assumption (3b), satisfies $\beta_{W_r} \le \beta$. Then using assumption (6) we obtain

$$ \delta_{V_r, W_r}^2 \le 1 - \frac{\alpha_{V_r, W_r}^2}{\beta_{W_r}^2} \le 1 - \frac{\alpha_{V_r, W_r}^2}{\beta^2} < 1. \qquad (12) $$

Furthermore for any $z' \in Z' \setminus \{0\}$ and $y \in W_r$, we have

$$ \begin{aligned} \langle s - Lu_r, z' \rangle &= \langle b - Au_r, A^{-*}L^*z' \rangle \stackrel{(5)}{=} \langle b - Au_r, A^{-*}L^*z' - y \rangle \\ &\le \|u - u_r\|_V \|L^*z' - A^*y\|_{V'}. \end{aligned} $$
Taking the infimum over $y \in W_r$, dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (9) thanks to (7). ■

The error bound (9) for the approximation of the variable of interest $s$ is the product of three terms:

(a) $\inf_{v \in V_r} \|u - v\|_V$, which suggests that the approximation space $V_r$ should be defined such that $u$ can be well approximated in $V_r$,

(b) $(1 - (\delta_{V_r, W_r})^2)^{-1/2}$, which suggests that the test space $W_r$ should be chosen such that any element of $V_r$ can be well approximated by an element of $R_V^{-1}A^*W_r$, and

(c) $\delta_{W_r}^L$, which suggests that any element of $\mathrm{range}(L^*)$ should be well approximated by an element of $A^*W_r$.

As already noticed in [19, Section 11.1], $W_r$ plays a double role: a test space for the definition of $u_r$ (point (b)) and an approximation space for the range of $A^{-*}L^*$ (point (c)).

**Remark 2.2.** The proposed Petrov-Galerkin projection method coincides with the interpo-
latory projection method used in the context of parametric dynamical systems (see [2, 4]).
Our analysis provides quasi-optimality results on $s(\xi)$ for any parameter value $\xi$. Also, the
condition $\delta_{V_r, W_r}(\xi) > 0$ ensures the invertibility of the reduced operator $A_r(\xi) : V_r \to W_r'$
defined by $\langle A_r(\xi)v, y \rangle = \langle A(\xi)u_r, y \rangle$ for all $v \in V_r$ and $y \in W_r$. In [2], the invertibility of
$A_r(\xi)$ is not discussed in the time-independent case.

**Remark 2.3** (Comparison with the Céa's Lemma). Under assumption (6), the classical
Céa's lemma states that

$$ \|u - u_r\|_V \le (1 + \frac{\beta}{\alpha_{V_r, W_r}}) \min_{v \in V_r} \|u - v\|_V. \quad (13) $$

It has been shown in [21] that this can be improved to

$$ \|u - u_r\|_V \le \frac{\beta}{\alpha_{V_r, W_r}} \min_{v \in V_r} \|u - v\|_V. \quad (14) $$

Noting that Equation (12) yields

$$ \frac{1}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \le \frac{\beta_{W_r}}{\alpha_{V_r, W_r}} \le \frac{\beta}{\alpha_{V_r, W_r}}, $$

we observe that (7) provides a sharper bound than in (14), where the constants differ by a factor $\beta_{W_r}/\beta$.
Taking the infimum over $y \in W_r$, dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (9) thanks to (7). ■

The error bound (9) for the approximation of the variable of interest $s$ is the product of three terms:

(a) $\inf_{v \in V_r} \|u - v\|_V$, which suggests that the approximation space $V_r$ should be defined such that $u$ can be well approximated in $V_r$,

(b) $(1 - (\delta_{V_r, W_r})^2)^{-1/2}$, which suggests that the test space $W_r$ should be chosen such that any element of $V_r$ can be well approximated by an element of $R_V^{-1}A^*W_r$, and

(c) $\delta_{W_r}^L$, which suggests that any element of $\text{range}(L^*)$ should be well approximated by an element of $A^*W_r$.

As already noticed in [19, Section 11.1], $W_r$ plays a double role: a test space for the definition of $u_r$ (point (b)) and an approximation space for the range of $A^{-*}L^*$ (point (c)).

**Remark 2.2.** The proposed Petrov-Galerkin projection method coincides with the interpo-
latory projection method used in the context of parametric dynamical systems (see [2, 4]).
Our analysis provides quasi-optimality results on $s(\xi)$ for any parameter value $\xi$. Also, the
condition $\delta_{V_r, W_r}(\xi) > 0$ ensures the invertibility of the reduced operator $A_r(\xi) : V_r \to W_r'$
defined by $\langle A_r(\xi)v, y \rangle = \langle A(\xi)u_r, y \rangle$ for all $v \in V_r$ and $y \in W_r$. In [2], the invertibility of
$A_r(\xi)$ is not discussed in the time-independent case.

**Remark 2.3** (Comparison with the Céa's Lemma). Under assumption (6), the classical
Céa's lemma states that

$$ \|u - u_r\|_V \leq \left(1 + \frac{\beta}{\alpha_{V_r, W_r}}\right) \min_{v \in V_r} \|u - v\|_V. \quad (13) $$

It has been shown in [21] that this can be improved to

$$ \|u - u_r\|_V \leq \frac{\beta}{\alpha_{V_r, W_r}} \min_{v \in V_r} \|u - v\|_V. \quad (14) $$

Noting that Equation (12) yields

$$ \frac{1}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \leq \frac{\beta_{W_r}}{\alpha_{V_r, W_r}} \leq \frac{\beta}{\alpha_{V_r, W_r}}, $$

we observe that (7) provides a sharper bound than in (14), where the constants differ by a factor $\beta_{W_r}/\beta$.
**Remark 2.4** (Symmetric coercive case and compliant case). We suppose that A is a symmetric coercive operator, with $V = W$ and $\| \cdot \|_V = \| \cdot \|_W$ the norm induced by the operator A such that $R_V = A$. Then $\delta_{V_r,W_r}$ defined by (8) admits the following simple expression

$$ \delta_{V_r,W_r} = \min_{0 \neq v \in V_r} \max_{y \in W_r} \frac{\|v - y\|_V}{\|v\|_V}. $$

If the test space $W_r$ is defined by $W_r = V_r$, we obtain $\delta_{V_r,W_r} = 0$, and from (7), we obtain $u_r = u_r^\perp$. In other words, the standard Galerkin projection coincides with the orthogonal projection.

In the case where the variable of interest $s$ is scalar-valued, we have $Z = \mathbb{R}$ and $\mathcal{L}(V, Z) = V'$. The so-called compliant case corresponds to $Lv = \langle b, v \rangle$ for any $v \in V$. Then, by definition (10), we have

$$ \delta_{W_r}^L = \min_{v \in V_r} \|b - Av\|_{V'} = \min_{v \in V_r} \|u - v\|_V = \|u - u_r\|_V, $$

and thanks to (9), we recover the so-called “squared effect”

$$ |s - Lu_r| = \|s - Lu_r\|_Z \le \|u - u_r\|_V^2. $$

## 2.2 Primal-dual approach

We now extend the classical primal-dual approach [16] for the estimation of a vector-valued variable of interest.

Let us introduce the dual variable $Q \in \mathcal{L}(Z', W)$ defined by $A^*Q = L^*$. The relation

$$ s = Lu = Q^*Au = Q^*b $$

shows that the variable of interest can be exactly determined if either the primal variable $u$ or the dual variable $Q$ is known.

Now, for given approximations $\tilde{u}$ of $u$ and $\tilde{Q}$ of $Q$, we define the approximation $\tilde{s}$ of $s$ by

$$ \tilde{s} = L\tilde{u} + \tilde{Q}^*(b - A\tilde{u}), \quad (15) $$

where $L\tilde{u}$ is the standard estimation of the variable of interest and where $\tilde{Q}^*(b - A\tilde{u})$ is a correction using the approximation of the dual variable. The following proposition provides an error bound on the variable of interest, which is a generalization of the classical error bound for scalar-valued variables of interest (see [16]) to vector-valued variables of interest.
**Remark 2.4** (Symmetric coercive case and compliant case). We suppose that A is a symmetric coercive operator, with $V = W$ and $\|\cdot\|_V = \|\cdot\|_W$ the norm induced by the operator A such that $R_V = A$. Then $\delta_{V_r,W_r}$ defined by (8) admits the following simple expression

$$ \delta_{V_r,W_r} = \min_{0 \neq v \in V_r} \max_{y \in W_r} \frac{\|v - y\|_V}{\|v\|_V}. $$

If the test space $W_r$ is defined by $W_r = V_r$, we obtain $\delta_{V_r,W_r} = 0$, and from (7), we obtain $u_r = u_r^\perp$. In other words, the standard Galerkin projection coincides with the orthogonal projection.

In the case where the variable of interest $s$ is scalar-valued, we have $Z = \mathbb{R}$ and $\mathcal{L}(V, Z) = V'$. The so-called compliant case corresponds to $Lv = \langle b, v \rangle$ for any $v \in V$. Then, by definition (10), we have

$$ \delta_{W_r}^L = \min_{v \in V_r} \|b - Av\|_{V'} = \min_{v \in V_r} \|u - v\|_V = \|u - u_r\|_V, $$

and thanks to (9), we recover the so-called “squared effect”

$$ |s - Lu_r| = \|s - Lu_r\|_Z \le \|u - u_r\|_V^2. $$

## 2.2 Primal-dual approach

We now extend the classical primal-dual approach [16] for the estimation of a vector-valued variable of interest.

Let us introduce the dual variable $Q \in \mathcal{L}(Z', W)$ defined by $A^*Q = L^*$. The relation

$$ s = Lu = Q^*Au = Q^*b $$

shows that the variable of interest can be exactly determined if either the primal variable $u$ or the dual variable $Q$ is known.

Now, for given approximations $\tilde{u}$ of $u$ and $\tilde{Q}$ of $Q$, we define the approximation $\tilde{s}$ of $s$ by

$$ \tilde{s} = L\tilde{u} + \tilde{Q}^*(b - A\tilde{u}), \quad (15) $$

where $L\tilde{u}$ is the standard estimation of the variable of interest and where $\tilde{Q}^*(b - A\tilde{u})$ is a correction using the approximation of the dual variable. The following proposition provides an error bound on the variable of interest, which is a generalization of the classical error bound for scalar-valued variables of interest (see [16]) to vector-valued variables of interest.
**Proposition 2.5.** The approximation $\tilde{s}$ of $s$ defined by (15) satisfies

$$ \|s - \tilde{s}\|_{Z} \le \|u - \tilde{u}\|_{V} \|L^{*} - A^{*}\tilde{Q}\|_{Z' \to V'}, \quad (16) $$

where

$$ \|L^* - A^*\tilde{Q}\|_{Z' \to V'} = \sup_{0 \neq z' \in Z'} \frac{\|(L^* - A^*\tilde{Q})z'\|_{V'}}{\|z'\|_{Z'}}. \quad (17) $$

*Proof.* For any $z' \in Z'$, we have

$$ \begin{aligned} \langle s - \tilde{s}, z' \rangle &= \langle Lu - L\tilde{u} - \tilde{Q}^*(b - A\tilde{u}), z' \rangle = \langle (L - \tilde{Q}^*A)(u - \tilde{u}), z' \rangle \\ &= \langle u - \tilde{u}, (L^* - A^*\tilde{Q})z' \rangle \le \|u - \tilde{u}\|_V \|(L^* - A^*\tilde{Q})z'\|_{V'}. \end{aligned} $$

Dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (16). ■

In practice, the approximation $\tilde{u}$ can be defined as the Petrov-Galerkin projection $u_r$ of $u$ on a given approximation space $V_r$ with a given test space $W_r$, see equation (5). For the approximation $\tilde{Q}$ of $Q \in \mathcal{L}(Z', W)$, the bound (16) suggests that $\|L^* - A^*\tilde{Q}\|_{Z' \to V'}$ should be small. We then propose to choose $\tilde{Q}$ as a solution of

$$ \min_{\tilde{Q} \in \mathcal{L}(Z', W_k^Q)} \|L^* - A^*\tilde{Q}\|_{Z' \to V'}, \quad (18) $$

where $W_k^Q \subset W$ is a given approximation space (different from $W_r$). The next proposition shows how to construct a solution of (18).

**Proposition 2.6.** The operator $Q_k : Z' \to W_k^Q$ defined for $z' \in Z'$ by

$$ Q_k z' = \arg \min_{y_k \in W_k^Q} \|L^* z' - A^* y_k\|_{V'} \quad (19) $$

is linear and is a solution of (18). Moreover $Q_k z' \in W_k^Q$ is characterized by

$$ \langle L^* z' - A^* Q_k z', R_V^{-1} A^* y_k \rangle = 0, \quad \forall y_k \in W_k^Q. \quad (20) $$

*Proof.* We easily prove that the optimization problem (19) admits a unique solution which depends linearly and continuously on $z'$, so that $Q_k$ defined by (19) is a linear operator in $\mathcal{L}(Z', W_k^Q)$. Equation (20) is the Euler equation associated to the minimization problem (19). Furthermore for any $\tilde{Q} \in \mathcal{L}(Z', W_k^Q)$ and $z' \in Z' \setminus \{0\}$, we have

$$ \frac{\|L^* z' - A^* Q_k z'\|_{V'}}{\|z'\|_{Z'}} \stackrel{(19)}{\le} \frac{\|L^* z' - A^* \tilde{Q} z'\|_{V'}}{\|z'\|_{Z'}} \le \|L^* - A^* \tilde{Q}\|_{Z' \to V'}. $$
**Proposition 2.5.** The approximation $\tilde{s}$ of $s$ defined by (15) satisfies

$$ \|s - \tilde{s}\|_{Z} \le \|u - \tilde{u}\|_{V} \|L^{*} - A^{*}\tilde{Q}\|_{Z' \to V'}, \quad (16) $$

where

$$ \|L^* - A^*\tilde{Q}\|_{Z' \to V'} = \sup_{0 \neq z' \in Z'} \frac{\|(L^* - A^*\tilde{Q})z'\|_{V'}}{\|z'\|_{Z'}}. \quad (17) $$

*Proof.* For any $z' \in Z'$, we have

$$ \begin{aligned} \langle s - \tilde{s}, z' \rangle &= \langle Lu - L\tilde{u} - \tilde{Q}^*(b - A\tilde{u}), z' \rangle = \langle (L - \tilde{Q}^*A)(u - \tilde{u}), z' \rangle \\ &= \langle u - \tilde{u}, (L^* - A^*\tilde{Q})z' \rangle \le \|u - \tilde{u}\|_V \|(L^* - A^*\tilde{Q})z'\|_{V'}. \end{aligned} $$

Dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (16). ■

In practice, the approximation $\tilde{u}$ can be defined as the Petrov-Galerkin projection $u_r$ of $u$ on a given approximation space $V_r$ with a given test space $W_r$, see equation (5). For the approximation $\tilde{Q}$ of $Q \in \mathcal{L}(Z', W)$, the bound (16) suggests that $\|L^* - A^*\tilde{Q}\|_{Z' \to V'}$ should be small. We then propose to choose $\tilde{Q}$ as a solution of

$$ \min_{\tilde{Q} \in \mathcal{L}(Z', W_k^Q)} \|L^* - A^*\tilde{Q}\|_{Z' \to V'}, \quad (18) $$

where $W_k^Q \subset W$ is a given approximation space (different from $W_r$). The next proposition shows how to construct a solution of (18).

**Proposition 2.6.** The operator $Q_k : Z' \to W_k^Q$ defined for $z' \in Z'$ by

$$ Q_k z' = \arg \min_{y_k \in W_k^Q} \|L^* z' - A^* y_k\|_{V'} \quad (19) $$

is linear and is a solution of (18). Moreover $Q_k z' \in W_k^Q$ is characterized by

$$ \langle L^* z' - A^* Q_k z', R_V^{-1} A^* y_k \rangle = 0, \quad \forall y_k \in W_k^Q. \quad (20) $$

*Proof.* We easily prove that the optimization problem (19) admits a unique solution which depends linearly and continuously on $z'$, so that $Q_k$ defined by (19) is a linear operator in $\mathcal{L}(Z', W_k^Q)$. Equation (20) is the Euler equation associated to the minimization problem (19). Furthermore for any $\tilde{Q} \in \mathcal{L}(Z', W_k^Q)$ and $z' \in Z' \setminus \{0\}$, we have

$$ \frac{\|L^* z' - A^* Q_k z'\|_{V'}}{\|z'\|_{Z'}} \stackrel{(19)}{\le} \frac{\|L^* z' - A^* \tilde{Q} z'\|_{V'}}{\|z'\|_{Z'}} \le \|L^* - A^* \tilde{Q}\|_{Z' \to V'}. $$
Taking the supremum over $z' \in Z' \setminus \{0\}$ and then the infimum over $\tilde{Q} \in \mathcal{L}(Z', W_k^Q)$, we obtain that $\|L^* - A^*Q_k\|_{Z' \to V'} \le \inf_{\tilde{Q} \in \mathcal{L}(Z', W_k^Q)} \|L^* - A^*\tilde{Q}\|_{Z' \to V'}$, which means that $Q_k \in \mathcal{L}(Z', W_k^Q)$ is a solution of (18). ■

In practice, for computing the approximation of the variable of interest (15) with $\tilde{Q} = Q_k$, we only need to compute $Q_k^*(b - Au_r)$. The following lemma shows how this can be performed without computing the operator $Q_k$.

**Lemma 2.7.** Let $Q_k$ be defined by (19). Then for $r = b - Au_r \in W'$,

$$Q_k^*r = LR_V^{-1}A^*y_k^{\perp}, \qquad (21)$$

where $y_k^{\perp} \in W_k^Q$ is defined by

$$\langle AR_V^{-1}A^*y_k^{\perp}, y_k \rangle = \langle r, y_k \rangle, \quad \forall y_k \in W_k^Q. \qquad (22)$$

*Proof.* For any $z' \in Z'$, since $Q_k z' \in W_k^Q$, we have

$$\langle Q_k z', AR_V^{-1}A^*y_k^{\perp} \rangle \stackrel{(22)}{=} \langle Q_k z', r \rangle. \qquad (23)$$

Furthermore, by definition of $Q_k$ we have

$$\langle Q_k z', AR_V^{-1}A^*y_k^{\perp} \rangle \stackrel{(20)}{=} \langle L^*z', R_V^{-1}A^*y_k^{\perp} \rangle. \qquad (24)$$

Combining (23) and (24), we obtain $\langle z', Q_k^*r \rangle = \langle z', LR_V^{-1}A^*y_k^{\perp} \rangle$ for all $z' \in Z'$, which concludes the proof. ■

We give now a new bound of the error on the variable of interest.

**Proposition 2.8.** Let $\tilde{u} = u_r$ be the Petrov-Galerkin projection defined by (5) and let $\tilde{Q} = Q_k$ be defined by (20). Then the approximation $\tilde{s}$ defined by (15) satisfies

$$\|s - \tilde{s}\|_Z \le \delta_{W_k^Q}^L \min_{y_k \in W_k^Q} \|u - u_r - R_V^{-1}A^*y_k\|_V, \qquad (25)$$

where

$$\delta_{W_k^Q}^L = \sup_{0 \neq z' \in Z'} \min_{y_k \in W_k^Q} \frac{\|L^*z' - A^*y\|_{V'}}{\|z'\|_{Z'}}. \qquad (26)$$
Taking the supremum over $z' \in Z' \setminus \{0\}$ and then the infimum over $\tilde{Q} \in \mathcal{L}(Z', W_k^Q)$, we obtain that $\|L^* - A^*Q_k\|_{Z' \to V'} \le \inf_{\tilde{Q} \in \mathcal{L}(Z', W_k^Q)} \|L^* - A^*\tilde{Q}\|_{Z' \to V'}$, which means that $Q_k \in \mathcal{L}(Z', W_k^Q)$ is a solution of (18). ■

In practice, for computing the approximation of the variable of interest (15) with $\tilde{Q} = Q_k$, we only need to compute $Q_k^*(b - Au_r)$. The following lemma shows how this can be performed without computing the operator $Q_k$.

**Lemma 2.7.** Let $Q_k$ be defined by (19). Then for $r = b - Au_r \in W'$,

$$Q_k^*r = LR_V^{-1}A^*y_k^{\perp}, \qquad (21)$$

where $y_k^{\perp} \in W_k^Q$ is defined by

$$\langle AR_V^{-1}A^*y_k^{\perp}, y_k \rangle = \langle r, y_k \rangle, \quad \forall y_k \in W_k^Q. \qquad (22)$$

*Proof.* For any $z' \in Z'$, since $Q_k z' \in W_k^Q$, we have

$$\langle Q_k z', AR_V^{-1}A^*y_k^{\perp} \rangle \stackrel{(22)}{=} \langle Q_k z', r \rangle. \qquad (23)$$

Furthermore, by definition of $Q_k$ we have

$$\langle Q_k z', AR_V^{-1}A^*y_k^{\perp} \rangle \stackrel{(20)}{=} \langle L^*z', R_V^{-1}A^*y_k^{\perp} \rangle. \qquad (24)$$

Combining (23) and (24), we obtain $\langle z', Q_k^*r \rangle = \langle z', LR_V^{-1}A^*y_k^{\perp} \rangle$ for all $z' \in Z'$, which concludes the proof. ■

We give now a new bound of the error on the variable of interest.

**Proposition 2.8.** Let $\tilde{u} = u_r$ be the Petrov-Galerkin projection defined by (5) and let $\tilde{Q} = Q_k$ be defined by (20). Then the approximation $\tilde{s}$ defined by (15) satisfies

$$\|s - \tilde{s}\|_{Z} \le \delta_{W_k^Q}^L \min_{y_k \in W_k^Q} \|u - u_r - R_V^{-1}A^*y_k\|_{V}, \qquad (25)$$

where

$$\delta_{W_k^Q}^L = \sup_{0 \neq z' \in Z'} \min_{y_k \in W_k^Q} \frac{\|L^*z' - A^*y\|_{V'}}{\|z'\|_{Z'}}. \qquad (26)$$
Moreover,

$$ \|s - \tilde{s}\|_{Z} \le \frac{\delta_{W_k}^L}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_{V}. \quad (27) $$

*Proof.* For any $z' \in Z'$, and for any $y_k \in W_k^Q$ we have

$$ \begin{aligned} \langle s - \tilde{s}, z' \rangle &\stackrel{(15)}{=} \langle u - u_r, (L^* - A^*Q_k)z' \rangle \\ &\stackrel{(20)}{=} \langle u - u_r - R_V^{-1}A^*y_k, (L^* - A^*Q_k)z' \rangle \\ &\le \|u - u_r - R_V^{-1}A^*y_k\|_V \|(L^* - A^*Q_k)z'\|_{V'}. \end{aligned} \quad (28) $$

From (19), we have

$$ \|(L^* - A^*Q_k)z'\|_{V'} = \min_{y_k \in W_k^Q} \|L^*z' - A^*y_k\|_{V'}. $$

Dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$ in (28), we obtain

$$ \|s - \tilde{s}\|_{Z} \le \delta_{W_k}^L \|u - u_r - R_V^{-1} A^* y_k \|_{V} $$

Then, taking the minimum over $y_k \in W_k^Q$, we obtain (25). Finally, taking $y_k = 0$ in (25), we obtain (27) from (7). ■

**Remark 2.9.** Observing that $\delta_{W_k}^L \le \|L^* - A^*Q_k\|_{V'}$ and $\min_{y_k \in W_k^Q} \|u - u_r - R_V^{-1}A^*y_k\|_V \le \|u - u_r\|_V$, Proposition 2.8 provides a sharper bound of the error on the variable of interest by taking advantage of the orthogonality property (20).

## 2.3 Projection based on a saddle point problem

In this section we extend the method proposed in [8] for the approximation of (vector-valued) variables of interest. The idea is to define the projection of $u$ on the reduced space $V_r$ by means of a saddle point problem. We first define and analyze this saddle point problem. Then we use the solution of this problem for the estimation of the variable of interest. ” Let us equip $W$ with a norm $\|\cdot\|_W$ such that the relation $\|y\|_W = \|A^*y\|_{V'}$ holds for any $y \in W$, which is equivalent to the following relation between the Riesz maps $R_W$ and $R_V$:

$$ R_W = A R_V^{-1} A^*. \quad (29) $$
Moreover,

$$ \|s - \tilde{s}\|_{Z} \leq \frac{\delta_{W_k}^L}{\sqrt{1 - (\delta_{V_r, W_r})^2}} \min_{v \in V_r} \|u - v\|_{V}. \qquad (27) $$

*Proof.* For any $z' \in Z'$, and for any $y_k \in W_k^Q$ we have

$$ \begin{aligned} \langle s - \tilde{s}, z' \rangle &\stackrel{(15)}{=} \langle u - u_r, (L^* - A^*Q_k)z' \rangle \\ &\stackrel{(20)}{=} \langle u - u_r - R_V^{-1}A^*y_k, (L^* - A^*Q_k)z' \rangle \\ &\leq \|u - u_r - R_V^{-1}A^*y_k\|_V \|(L^* - A^*Q_k)z'\|_{V'}. \end{aligned} \qquad (28) $$

From (19), we have

$$ \|(L^* - A^*Q_k)z'\|_{V'} = \min_{y_k \in W_k^Q} \|L^*z' - A^*y_k\|_{V'}. $$

Dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$ in (28), we obtain

$$ \|s - \tilde{s}\|_{Z} \leq \delta_{W_k}^L \|u - u_r - R_V^{-1} A^* y_k\|_{V} $$

Then, taking the minimum over $y_k \in W_k^Q$, we obtain (25). Finally, taking $y_k = 0$ in (25), we obtain (27) from (7). ■

**Remark 2.9.** Observing that $\delta_{W_k}^L \leq \|L^* - A^*Q_k\|_{V'}$ and $\min_{y_k \in W_k^Q} \|u - u_r - R_V^{-1}A^*y_k\|_V \leq \|u - u_r\|_V$, Proposition 2.8 provides a sharper bound of the error on the variable of interest by taking advantage of the orthogonality property (20).

## 2.3 Projection based on a saddle point problem

In this section we extend the method proposed in [8] for the approximation of (vector-valued) variables of interest. The idea is to define the projection of $u$ on the reduced space $V_r$ by means of a saddle point problem. We first define and analyze this saddle point problem. Then we use the solution of this problem for the estimation of the variable of interest. ” Let us equip $W$ with a norm $\|\cdot\|_W$ such that the relation $\|y\|_W = \|A^*y\|_{V'}$ holds for any $y \in W$, which is equivalent to the following relation between the Riesz maps $R_W$ and $R_V$:

$$ R_W = A R_V^{-1} A^*. \qquad (29) $$
The orthogonal projection $u_r^\perp$ of $u$ on $V_r$ satisfies

$$
\begin{align*}
\|u - u_r^\perp\|_V &= \min_{v \in V_r} \|u - v\|_V = \min_{v \in V_r} \sup_{0 \neq w \in V} \frac{\langle u - v, R_V w \rangle}{\|w\|_V} \\
&= \min_{v \in V_r} \max_{0 \neq w \in V} \frac{\langle Av - b, A^{-*} R_V w \rangle}{\|w\|_V} = \min_{v \in V_r} \max_{0 \neq y \in W} \frac{\langle Av - b, y \rangle}{\|R_V^{-1} A^* y\|_V} \\
&= \min_{v \in V_r} \max_{0 \neq y \in W} \frac{\langle Av - b, y \rangle}{\|y\|_W}.
\end{align*}
$$

Starting from this observation, we introduce a subspace $T_p \subset W$ of dimension $p$ and we define the projection $u_{r,p}$ in $V_r$ as the solution of the saddle point problem

$$
\min_{v \in V_r} \max_{w \in T_p} \frac{\langle Av - b, w \rangle}{\|w\|_W} \quad (30)
$$

In the following proposition, we prove the well-posedness of (30) under the condition (discrete inf-sup condition)

$$
\inf_{0 \neq v \in V_r} \sup_{0 \neq y \in T_p} \frac{\langle Av, y \rangle}{\|v\|_V \|y\|_W} =: \alpha_{V_r, T_p} > 0, \quad (31)
$$

and we provide a practical characterization of $u_{r,p}$.

**Proposition 2.10.** Under assumption (31), there exists a unique solution $(u_{r,p}, y_{r,p})$ in $V_r \times T_p$ to

$$
\langle R_W y_{r,p}, y \rangle + \langle A u_{r,p}, y \rangle = \langle b, y \rangle \quad \forall y \in T_p, \quad (32a)
$$

$$
\langle A^* y_{r,p}, v \rangle = 0 \quad \forall v \in V_r, \quad (32b)
$$

and $(u_{r,p}, \frac{y_{r,p}}{\|y_{r,p}\|_W})$ is the unique solution of (30).

*Proof.* Since the Riesz map $R_W$ defined by (29) is coercive and under the discrete inf-sup condition (31) on operator $A$, Theorem 2.34 of [9] gives that (32) is a well-posed problem whose solution $(u_{r,p}, y_{r,p})$ is the unique solution of the saddle-point problem

$$
\min_{v \in V_r} \max_{y \in T_p} -\frac{1}{2} \langle R_W y, y \rangle + \langle b, y \rangle - \langle Av, y \rangle.
$$

Denoting $y = \lambda w$ with $\|w\|_W = 1$, this saddle point problem is equivalent to

$$
\min_{v \in V_r} \max_{w \in T_p} \max_{\lambda \in \mathbb{R}} -\frac{1}{2} \lambda^2 + \lambda \langle b - Av, w \rangle = \min_{v \in V_r} \max_{w \in T_p} \frac{1}{2} \langle Av - b, w \rangle^2,
$$

which coincides with problem (30).
The orthogonal projection $u_r^\perp$ of $u$ on $V_r$ satisfies

$$
\begin{align*}
\|u - u_r^\perp\|_V &= \min_{v \in V_r} \|u - v\|_V = \min_{v \in V_r} \sup_{0 \neq w \in V} \frac{\langle u - v, R_V w \rangle}{\|w\|_V} \\
&= \min_{v \in V_r} \max_{0 \neq w \in V} \frac{\langle Av - b, A^{-*} R_V w \rangle}{\|w\|_V} = \min_{v \in V_r} \max_{0 \neq y \in W} \frac{\langle Av - b, y \rangle}{\|R_V^{-1} A^* y\|_V} \\
&= \min_{v \in V_r} \max_{0 \neq y \in W} \frac{\langle Av - b, y \rangle}{\|y\|_W}.
\end{align*}
$$

Starting from this observation, we introduce a subspace $T_p \subset W$ of dimension $p$ and we
define the projection $u_{r,p}$ in $V_r$ as the solution of the saddle point problem

$$
\min_{v \in V_r} \max_{w \in T_p} \frac{\langle Av - b, w \rangle}{\|w\|_W} \quad (30)
$$

In the following proposition, we prove the well-posedness of (30) under the condition (discrete inf-sup condition)

$$
\inf_{0 \neq v \in V_r} \sup_{0 \neq y \in T_p} \frac{\langle Av, y \rangle}{\|v\|_V \|y\|_W} =: \alpha_{V_r, T_p} > 0, \quad (31)
$$

and we provide a practical characterization of $u_{r,p}$.

**Proposition 2.10.** Under assumption (31), there exists a unique solution $(u_{r,p}, y_{r,p})$ in $V_r \times T_p$ to

$$
\langle R_W y_{r,p}, y \rangle + \langle A u_{r,p}, y \rangle = \langle b, y \rangle \quad \forall y \in T_p, \quad (32a)
$$

$$
\langle A^* y_{r,p}, v \rangle = 0 \quad \forall v \in V_r, \qquad (32b)
$$

and $(u_{r,p}, \frac{y_{r,p}}{\|y_{r,p}\|_W})$ is the unique solution of (30).

Proof. Since the Riesz map $R_W$ defined by (29) is coercive and under the discrete inf-sup condition (31) on operator $A$, Theorem 2.34 of [9] gives that (32) is a well-posed problem whose solution $(u_{r,p}, y_{r,p})$ is the unique solution of the saddle-point problem

$$
\min_{v \in V_r} \max_{y \in T_p} -\frac{1}{2} \langle R_W y, y \rangle + \langle b, y \rangle - \langle Av, y \rangle.
$$

Denoting $y = \lambda w$ with $\|w\|_W = 1$, this saddle point problem is equivalent to

$$
\min_{v \in V_r} \max_{w \in T_p} \max_{\lambda \in \mathbb{R}} -\frac{1}{2} \lambda^2 + \lambda \langle b - Av, w \rangle = \min_{v \in V_r} \max_{w \in T_p} \frac{1}{2} \langle Av - b, w \rangle^2, \quad \|w\|_W = 1
$$

which coincides with problem (30).
The following proposition provides a quasi-optimality result for the projection $u_{r,p} \in V_r$ of $u$ onto $V_r$.

**Proposition 2.11.** Under assumption (31), the solution $u_{r,p}$ of (32) is such that

$$ \|u - u_{r,p}\|_V \le \frac{1}{\sqrt{1 - \delta_{V_r, T_p}^2}} \min_{v \in V_r} \|u - v\|_V, \quad (33) $$

where

$$ \delta_{V_r, T_p} = \max_{0 \neq v \in V_r} \min_{y \in T_p} \frac{\|v - R_V^{-1} A^* y\|_V}{\|v\|_V} \quad (34) $$

is such that

$$ \delta_{V_r, T_p}^2 = 1 - \alpha_{V_r, T_p}^2 < 1. \quad (35) $$

*Proof.* Let $(u_{r,p}, y_{r,p})$ be the solution of (32). For any $v \in V_r$ and $y \in T_p$, we have

$$ \begin{align} \langle u_r^\perp - u_{r,p}, R_V v \rangle &\stackrel{(4)}{=} \langle u - u_{r,p}, R_V v \rangle = \langle b - A u_{r,p}, A^{-*} R_V v \rangle \nonumber \\ &\stackrel{(32a)}{=} \langle b - A u_{r,p}, A^{-*} R_V v - y \rangle + \langle R_W y_{r,p}, y \rangle \nonumber \\ &\stackrel{(32b)}{=} \langle b - A u_{r,p}, A^{-*} R_V v - y \rangle - \langle R_W y_{r,p}, A^{-*} R_V v - y \rangle \nonumber \\ &\stackrel{(29)}{=} \langle u - u_{r,p} - R_V^{-1} A^* y_{r,p}, R_V v - A^* y \rangle \nonumber \\ &\le \|u - u_{r,p} - R_V^{-1} A^* y_{r,p}\|_V \|R_V v - A^* y\|_{V'}. \tag{36} \end{align} $$

Equation (32a) implies that

$$ y_{r,p} = \arg \min_{y \in W} \|R_W^{-1}(b - Au_{r,p}) - y\|_W, $$

so that $\|R_W^{-1}(b - Au_{r,p}) - y_{r,p}\|_W \le \|R_W^{-1}(b - Au_{r,p})\|_W$. Using (29), it follows

$$ \|u - u_{r,p} - R_V^{-1} A^* y_{r,p}\|_V \le \|u - u_{r,p}\|_V, \quad (37) $$

Using (37) in (36), taking the infimum over $y \in T_p$, dividing by $\|v\|_V$ and taking the supremum over $v \in V_r \setminus \{0\}$, we obtain

$$ \|u_r^\perp - u_{r,p}\|_V \le \delta_{V_r, T_p} \|u - u_{r,p}\|. $$

From (11) (with $W_r$ replaced by $T_p$) and (29) (which implies $\|A^*y\|_{V'} = \|y\|_W$), we obtain (35). Then, from the definition of $u_r^\perp$, we have

$$ \|u - u_{r,p}\|^2 = \|u - u_r^*\|^2 + \|u_r^* - u_{r,p}\|^2 \le \|u - u_r^*\|^2 + \delta_{V_r, T_p}^2 \|u - u_{r,p}\|^2, $$
The following proposition provides a quasi-optimality result for the projection $u_{r,p} \in V_r$ of $u$ onto $V_r$.

**Proposition 2.11.** Under assumption (31), the solution $u_{r,p}$ of (32) is such that

$$ \|u - u_{r,p}\|_V \le \frac{1}{\sqrt{1 - \delta_{V_r, T_p}^2}} \min_{v \in V_r} \|u - v\|_V, \quad (33) $$

where

$$ \delta_{V_r, T_p} = \max_{0 \neq v \in V_r} \min_{y \in T_p} \frac{\|v - R_V^{-1} A^* y\|_V}{\|v\|_V} \quad (34) $$

is such that

$$ \delta_{V_r, T_p}^2 = 1 - \alpha_{V_r, T_p}^2 < 1. \quad (35) $$

*Proof.* Let $(u_{r,p}, y_{r,p})$ be the solution of (32). For any $v \in V_r$ and $y \in T_p$, we have

$$ \begin{align} \langle u_r^\perp - u_{r,p}, R_V v \rangle &\stackrel{(4)}{=} \langle u - u_{r,p}, R_V v \rangle = \langle b - A u_{r,p}, A^{-*} R_V v \rangle \nonumber \\ &\stackrel{(32a)}{=} \langle b - A u_{r,p}, A^{-*} R_V v - y \rangle + \langle R_W y_{r,p}, y \rangle \nonumber \\ &\stackrel{(32b)}{=} \langle b - A u_{r,p}, A^{-*} R_V v - y \rangle - \langle R_W y_{r,p}, A^{-*} R_V v - y \rangle \nonumber \\ &\stackrel{(29)}{=} \langle u - u_{r,p} - R_V^{-1} A^* y_{r,p}, R_V v - A^* y \rangle \nonumber \\ &\le \|u - u_{r,p} - R_V^{-1} A^* y_{r,p}\|_V \|R_V v - A^* y\|_{V'}. \tag{36} \end{align} $$

Equation (32a) implies that

$$ y_{r,p} = \arg \min_{y \in W} \|R_W^{-1}(b - Au_{r,p}) - y\|_W, $$

so that $\|R_W^{-1}(b - Au_{r,p}) - y_{r,p}\|_W \le \|R_W^{-1}(b - Au_{r,p})\|_W$. Using (29), it follows

$$ \|u - u_{r,p} - R_V^{-1} A^* y_{r,p}\|_V \le \|u - u_{r,p}\|_V, \quad (37) $$

Using (37) in (36), taking the infimum over $y \in T_p$, dividing by $\|v\|_V$ and taking the supremum over $v \in V_r \setminus \{0\}$, we obtain

$$ \|u_r^\perp - u_{r,p}\|_V \le \delta_{V_r, T_p} \|u - u_{r,p}\|. $$

From (11) (with $W_r$ replaced by $T_p$) and (29) (which implies $\|A^*y\|_{V'} = \|y\|_W$), we obtain (35). Then, from the definition of $u_r^\perp$, we have

$$ \|u - u_{r,p}\|^2 = \|u - u_r^*\|^2 + \|u_r^* - u_{r,p}\|^2 \le \|u - u_r^*\|^2 + \delta_{V_r, T_p}^2 \|u - u_{r,p}\|^2, $$
from which we deduce (33).

From the definition (34) of $\delta_{V_r, T_p}$, we easily deduce the following corollary.

**Corollary 2.12.** If $T_p$ is such that $R_W^{-1}AV_r \subset T_p$ (or equivalently $V_r \subset R_V^{-1}A^*T_p$) then $\delta_{V_r, T_p} = 0$ and $u_{r,p}$ coincides with the best approximation $u_r^\perp$ of $u$ in $V_r$.

**Remark 2.13.** Note that (33) and (35) give

$$ \|u - u_{r,p}\|_V \le \frac{1}{\alpha_{V_r, T_p}} \min_{v \in V_r} \|u - v\|_V, $$

which is sharper than the classical error bound obtained by the Cea's lemma

$$ \|u - u_{r,p}\|_V \le (1 + \frac{1}{\alpha_{V_r, T_p}}) \min_{v \in V_r} \|u - v\|_V. $$

Now, we consider the approximation $\tilde{s}$ of $s$ defined by

$$ \tilde{s} = Lu_{r,p} + LR_V^{-1}A^*y_{r,p}, \quad (38) $$

where $(u_{r,p}, y_{r,p}) \in V_r \times T_p$ is the solution of the saddle point problem (32). The following proposition provides an error bound for the approximation of the variable of interest.

**Proposition 2.14.** The approximation $\tilde{s}$ defined by (38) satisfies

$$ \|s - \tilde{s}\|_Z \le \delta_{T_p}^L \|u - u_{r,p} - R_V^{-1}A^*y_{r,p}\|_V, \quad (39) $$

with

$$ \delta_{T_p}^L = \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L^*z' - A^*y\|_{Z'}}{\|z'\|_{Z'}}, \quad (40) $$

and

$$ \|s - \tilde{s}\|_Z \le \frac{\delta_{T_p}^L}{\sqrt{1 - (\delta_{V_r, T_p})^2}} \min_{v \in V_r} \|u - v\|_V. \quad (41) $$

*Proof.* For any $z' \in Z'$ and $y \in T_p$, we have

$$ \begin{aligned} \langle s - \tilde{s}, z' \rangle &\stackrel{(38)}{=} \langle u - u_{r,p} - R_V^{-1}A^*y_{r,p}, L^*z' \rangle \\ &\stackrel{(32a)\&(29)}{=} \langle u - u_{r,p} - R_V^{-1}A^*y_{r,p}, L^*z' - A^*y \rangle \\ &\le \|u - u_{r,p} - R_V^{-1}A^*y_{r,p}\|_V \|L^*z' - A^*y\|_{V'}. \end{aligned} $$
from which we deduce (33).

From the definition (34) of $\delta_{V_r, T_p}$, we easily deduce the following corollary.

**Corollary 2.12.** If $T_p$ is such that $R_W^{-1}AV_r \subset T_p$ (or equivalently $V_r \subset R_V^{-1}A^*T_p$) then $\delta_{V_r, T_p} = 0$ and $u_{r,p}$ coincides with the best approximation $u_r^\perp$ of $u$ in $V_r$.

**Remark 2.13.** Note that (33) and (35) give

$$ \|u - u_{r,p}\|_V \le \frac{1}{\alpha_{V_r, T_p}} \min_{v \in V_r} \|u - v\|_V, $$

which is sharper than the classical error bound obtained by the Cea's lemma

$$ \|u - u_{r,p}\|_V \le (1 + \frac{1}{\alpha_{V_r, T_p}}) \min_{v \in V_r} \|u - v\|_V. $$

Now, we consider the approximation $\tilde{s}$ of $s$ defined by

$$ \tilde{s} = Lu_{r,p} + LR_V^{-1}A^*y_{r,p}, \quad (38) $$

where $(u_{r,p}, y_{r,p}) \in V_r \times T_p$ is the solution of the saddle point problem (32). The following proposition provides an error bound for the approximation of the variable of interest.

**Proposition 2.14.** The approximation $\tilde{s}$ defined by (38) satisfies

$$ \|s - \tilde{s}\|_Z \le \delta_{T_p}^L \|u - u_{r,p} - R_V^{-1}A^*y_{r,p}\|_V, \quad (39) $$

with

$$ \delta_{T_p}^L = \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L^*z' - A^*y\|_{Z'}}{\|z'\|_{Z'}}, \quad (40) $$

and

$$ \|s - \tilde{s}\|_Z \le \frac{\delta_{T_p}^L}{\sqrt{1 - (\delta_{V_r, T_p})^2}} \min_{v \in V_r} \|u - v\|_V. \quad (41) $$

*Proof.* For any $z' \in Z'$ and $y \in T_p$, we have

$$ \begin{align*} \langle s - \tilde{s}, z' \rangle &\stackrel{(38)}{=} \langle u - u_{r,p} - R_V^{-1}A^*y_{r,p}, L^*z' \rangle \\ &\stackrel{(32a)\&(29)}{=} \langle u - u_{r,p} - R_V^{-1}A^*y_{r,p}, L^*z' - A^*y \rangle \\ &\le \|u - u_{r,p} - R_V^{-1}A^*y_{r,p}\|_V \|L^*z' - A^*y\|_{V'}. \end{align*} $$
Taking the minimum over $y \in T_p$, dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (39). Finally, thanks to (39), (37) and (33), we obtain (41). ■

We observe that $T_p$ impacts both the quality of the projection of $u$ (via the constant $\delta_{V_r,T_p}$ in (33)) and the quality of the approximation of the variable of interest (via constants $\delta_{V_r,T_p}$ and $\delta_{T_p}^L$ in (41)). Then, we will consider for $T_p$ spaces of the form

$$T_p = W_r + W_k^Q, \qquad (42)$$

with $\dim(W_r) = r$. This implies

$$\delta_{T_p}^L \le \delta_{W_k^Q}^L \quad \text{and} \quad \delta_{V_r,T_p} \le \delta_{V_r,W_r},$$

so that the error bound (41) for the variable of interest is better than the error bound (27) of the primal-dual method with primal approximation space $V_r$, primal test space $W_r$ and dual approximation space $W_k^Q$. Therefore, we expect the approximation $u_{r,p}$ to be closer to the solution $u$ than the Petrov-Galerkin projection $u_r$. Also, the approximation of the quantity of interest is expected to be improved.

**Remark 2.15 (Symmetric coercive case).** Let us consider the case where A is symmetric and coercive, $R_V = R_W = A$ and $W_r = V_r$. The choice (42) implies that $V_r \subset T_p$, so that $T_p$ admits the orthogonal decomposition $T_p = V_r \oplus (T_p \cap V_r^\perp)$. Equation (32b) implies that $y_{r,p} \in T_p \cap V_r^\perp$. Let $t_{r,p} = y_{r,p} + u_{r,p} \in T_p$. Equation (32a) gives $\langle R_V t_{r,p}, y \rangle = \langle R_V u, y \rangle$ for all $y \in T_p$, which implies that $t_{r,p}$ is the orthogonal projection of $u$ on $T_p$, where $u_{r,p}$ and $y_{r,p}$ are the orthogonal projections of $u$ on $V_r$ and $T_p \cap V_r^\perp$ respectively. Furthermore, the approximation of the variable of interest (38) is given by $\tilde{s} = Lt_{r,p}$. We conclude that in this particular setting, the saddle point approach can be simply interpreted as an orthogonal projection of $u$ on the enriched space $T_p = V_r + W_k^Q$, followed by a standard estimation of the variable of interest.

## 3 Goal-oriented projections for parameter-dependent equations

We now consider a parameter-dependent equation $A(\xi)u(\xi) = b(\xi)$ where $\xi$ denotes a parameter taking values in a set $\Xi \subset \mathbb{R}^d$, $A(\xi) \in \mathcal{L}(V, W')$ and $b(\xi) \in W'$. The variable of interest is defined by $s(\xi) = L(\xi)u(\xi)$, with $L(\xi) \in \mathcal{L}(V, Z)$.
Taking the minimum over $y \in T_p$, dividing by $\|z'\|_{Z'}$ and taking the supremum over $z' \in Z' \setminus \{0\}$, we obtain (39). Finally, thanks to (39), (37) and (33), we obtain (41). ■

We observe that $T_p$ impacts both the quality of the projection of $u$ (via the constant $\delta_{V_r, T_p}$ in (33)) and the quality of the approximation of the variable of interest (via constants $\delta_{V_r, T_p}$ and $\delta_{T_p}^L$ in (41)). Then, we will consider for $T_p$ spaces of the form

$$T_p = W_r + W_k^Q, \quad (42)$$

with $\dim(W_r) = r$. This implies

$$\delta_{T_p}^L \le \delta_{W_k^Q}^L \quad \text{and} \quad \delta_{V_r, T_p} \le \delta_{V_r, W_r},$$

so that the error bound (41) for the variable of interest is better than the error bound (27) of the primal-dual method with primal approximation space $V_r$, primal test space $W_r$ and dual approximation space $W_k^Q$. Therefore, we expect the approximation $u_{r,p}$ to be closer to the solution $u$ than the Petrov-Galerkin projection $u_r$. Also, the approximation of the quantity of interest is expected to be improved.

**Remark 2.15 (Symmetric coercive case).** Let us consider the case where $A$ is symmetric and coercive, $R_V = R_W = A$ and $W_r = V_r$. The choice (42) implies that $V_r \subset T_p$, so that $T_p$ admits the orthogonal decomposition $T_p = V_r \oplus (T_p \cap V_r^\perp)$. Equation (32b) implies that $y_{r,p} \in T_p \cap V_r^\perp$. Let $t_{r,p} = y_{r,p} + u_{r,p} \in T_p$. Equation (32a) gives $\langle R_V t_{r,p}, y \rangle = \langle R_V u, y \rangle$ for all $y \in T_p$, which implies that $t_{r,p}$ is the orthogonal projection of $u$ on $T_p$, where $u_{r,p}$ and $y_{r,p}$ are the orthogonal projections of $u$ on $V_r$ and $T_p \cap V_r^\perp$ respectively. Furthermore, the approximation of the variable of interest (38) is given by $\tilde{s} = Lt_{r,p}$. We conclude that in this particular setting, the saddle point approach can be simply interpreted as an orthogonal projection of $u$ on the enriched space $T_p = V_r + W_k^Q$, followed by a standard estimation of the variable of interest.

# 3 Goal-oriented projections for parameter-dependent equations

We now consider a parameter-dependent equation $A(\xi)u(\xi) = b(\xi)$ where $\xi$ denotes a parameter taking values in a set $\Xi \subset \mathbb{R}^d$, $A(\xi) \in \mathcal{L}(V, W')$ and $b(\xi) \in W'$. The variable of interest is defined by $s(\xi) = L(\xi)u(\xi)$, with $L(\xi) \in \mathcal{L}(V, Z)$.
In Section 2, we have presented different projection methods for the estimation of the variable of interest which rely on the introduction of three spaces: the primal approximation space $V_r$, the primal test space $W_r$ and the dual approximation space $W_k^Q$. We recall that for the saddle point approach, we introduce the space $T_p = W_r + W_k^Q$. We adopt an offline/online strategy. Reduced (low-dimensional) spaces $V_r$, $W_r$ and $W_k^Q$ are constructed during the offline phase. Then, the projections on these reduced spaces and the evaluations of the variable of interest are rapidly computed for any parameter value $\xi \in \Xi$ during the online phase.

In Section 3.1, we will first consider the construction of the test space $W_r$. For scalar-valued variables of interest, reduced spaces $V_r$ and $W_k^Q$ are classically defined as the span of snapshots of the primal and dual solutions $u(\xi)$ and $Q(\xi)$. These snapshots can be selected at random, using samples drawn according a certain probability measure over $\Xi$, see e.g. [17]. Another popular method is to select the snapshots in a greedy way [7, 11, 19], with a uniform control of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ over $\Xi$. This method requires an estimation of the error on the variable of interest. In the same lines, we introduce error estimates for vector-valued variables of interest in Section 3.2, and we propose greedy algorithms for the construction of $V_r$ and $W_k^Q$ in Section 3.3.

## 3.1 Construction of the test space $W_r$

Assuming that the primal approximation space $V_r$ is given, we know from the previous section that $W_r$ should be chosen such that $\delta_{V_r, W_r}$ is as close to zero as possible (see Propositions 2.1, 2.8, 2.11 and 2.14). In the literature, $W_r = V_r$ is a common choice (standard Galerkin projection). When the operator $A(\xi)$ is symmetric and coercive, we can choose $W_r = V_r$ which is the optimal test space with respect to the norm induced by $A(\xi)$ (see Remark 2.4). However, this choice may lead to an inaccurate projection of the primal variable when the operator is ill-conditioned (i.e. $\frac{\beta}{\alpha} \gg 1$). In the case of non coercive operators, a parameter-dependent test space is generally defined by $W_r = W_r(\xi) = R_W^{-1}A(\xi)V_r$, where $R_W^{-1}A(\xi)$ is called the “supremizer operator” (see e.g. [20, 15]). This approach is no more than a minimal residual method since the resulting Petrov-Galerkin projection defined by (5) is $u_r(\xi) = \arg\min_{v_r \in V_r} \|A(\xi)v_r - b(\xi)\|_{W'}$. In Section 2.1, we have seen that the Petrov-Galerkin projection with an ideal test space

$$W_r(\xi) = A(\xi)^{-*} R_V V_r \quad (43)$$

coincides with the best approximation. Having a basis $v_1, \dots, v_r$ of $V_r$, the computation of this ideal parameter-dependent test space would require the computation of $A^{-*}(\xi)R_V v_i$
In Section 2, we have presented different projection methods for the estimation of the variable of interest which rely on the introduction of three spaces: the primal approximation space $V_r$, the primal test space $W_r$ and the dual approximation space $W_k^Q$. We recall that for the saddle point approach, we introduce the space $T_p = W_r + W_k^Q$. We adopt an offline/online strategy. Reduced (low-dimensional) spaces $V_r$, $W_r$ and $W_k^Q$ are constructed during the offline phase. Then, the projections on these reduced spaces and the evaluations of the variable of interest are rapidly computed for any parameter value $\xi \in \Xi$ during the online phase.

In Section 3.1, we will first consider the construction of the test space $W_r$. For scalar-valued variables of interest, reduced spaces $V_r$ and $W_k^Q$ are classically defined as the span of snapshots of the primal and dual solutions $u(\xi)$ and $Q(\xi)$. These snapshots can be selected at random, using samples drawn according a certain probability measure over $\Xi$, see e.g. [17]. Another popular method is to select the snapshots in a greedy way [7, 11, 19], with a uniform control of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ over $\Xi$. This method requires an estimation of the error on the variable of interest. In the same lines, we introduce error estimates for vector-valued variables of interest in Section 3.2, and we propose greedy algorithms for the construction of $V_r$ and $W_k^Q$ in Section 3.3.

## 3.1 Construction of the test space $W_r$

Assuming that the primal approximation space $V_r$ is given, we know from the previous section that $W_r$ should be chosen such that $\delta_{V_r, W_r}$ is as close to zero as possible (see Propositions 2.1, 2.8, 2.11 and 2.14). In the literature, $W_r = V_r$ is a common choice (standard Galerkin projection). When the operator $A(\xi)$ is symmetric and coercive, we can choose $W_r = V_r$ which is the optimal test space with respect to the norm induced by $A(\xi)$ (see Remark 2.4). However, this choice may lead to an inaccurate projection of the primal variable when the operator is ill-conditioned (i.e. $\frac{\beta}{\alpha} \gg 1$). In the case of non coercive operators, a parameter-dependent test space is generally defined by $W_r = W_r(\xi) = R_W^{-1}A(\xi)V_r$, where $R_W^{-1}A(\xi)$ is called the “supremizer operator” (see e.g. [20, 15]). This approach is no more than a minimal residual method since the resulting Petrov-Galerkin projection defined by (5) is $u_r(\xi) = \arg\min_{v_r \in V_r} \|A(\xi)v_r - b(\xi)\|_{W'}$. In Section 2.1, we have seen that the Petrov-Galerkin projection with an ideal test space

$$W_r(\xi) = A(\xi)^{-*} R_V V_r \quad (43)$$

coincides with the best approximation. Having a basis $v_1, \dots, v_r$ of $V_r$, the computation of this ideal parameter-dependent test space would require the computation of $A^{-*}(\xi)R_V v_i$
for all $1 \le i \le r$ for each parameter's value $\xi$, which is unfeasible in practice. Up to our knowledge, the only attempt to construct quasi-optimal test spaces for non symmetric and weakly coercive operators can be found in [8], where the authors proposed a greedy algorithm for the construction of a (parameter independent) test space which ensures the quasi-optimality constant to be uniformly bounded by an arbitrarily small constant. Here, we adopt an alternative approach where the (parameter-dependent) test space is defined by

$$W_r(\xi) = P_m(\xi)^* R_V V_r, \qquad (44)$$

where $P_m(\xi)$ is an interpolation of the inverse of $A(\xi)$ using $m$ interpolation points in the parameter set $\Xi$. In practice, when $A(\xi)$ is a matrix, algorithms developed in [22] can be used. This will be detailed later on. The underlying idea is to obtain a test space as close as possible to the ideal test space $A^{-*}(\xi)R_VV_r$ defined in (43). For $m = 0$, with $P_0(\xi) = R_V^{-1}$ by convention, we have $W_r = V_r$, which yields the standard Galerkin projection.

## 3.2 Error estimates for vector-valued variables of interest

In this section, we propose practical error estimates for the variable of interest, first for the primal-dual approach and then for the saddle point method.

### 3.2.1 Primal-dual approach

Given approximations $\tilde{u}$ and $\tilde{Q}$ of the primal solution $u$ and the dual solution $Q$ respectively, a standard approach is to start from the error bound

$$\|s(\xi) - \tilde{s}(\xi)\|_Z \le \|u(\xi) - \tilde{u}(\xi)\|_V \|L(\xi)^* - A(\xi)^* \tilde{Q}(\xi)\|_{Z' \to V'}$$

which is provided by Proposition 2.5. This suggests to measure the norm of the residuals associated to the primal and dual variables. In practice, we distinguish two cases.

In the case where the operator $A(\xi)$ is symmetric and coercive, it is natural to choose the parameter-dependent norm $\| \cdot \|_V$ as the one induced by the operator, i.e. $R_V = R_V(\xi) = A(\xi)$. However, neither the primal error $\|u(\xi) - \tilde{u}(\xi)\|_V = \langle b(\xi) - A(\xi)\tilde{u}(\xi), u(\xi) - \tilde{u}(\xi) \rangle$ nor the dual residual norm $\|L(\xi)^* - A(\xi)^* \tilde{Q}(\xi)\|_{Z' \to V'} = \sup_{\|z\|_Z=1} \|Q(\xi)z - \tilde{Q}(\xi)z\|_Z$ can be computed without computing the primal and dual solutions $u(\xi)$ and $Q(\xi)$. The classical way to circumvent this issue is to introduce a parameter-independent norm $\| \cdot \|_{V_0}$, which is in general the “natural” norm associated to the space $V$, and to measure residuals with the associated dual norm $\| \cdot \|_{V'_0}$. Here we assume that the operator $A(\xi)$ satisfies

$$\alpha(\xi) \|v\|_{V_0} \le \|A(\xi)v\|_{V'_0} \qquad (45)$$
for all $1 \le i \le r$ for each parameter's value $\xi$, which is unfeasible in practice. Up to our knowledge, the only attempt to construct quasi-optimal test spaces for non symmetric and weakly coercive operators can be found in [8], where the authors proposed a greedy algorithm for the construction of a (parameter independent) test space which ensures the quasi-optimality constant to be uniformly bounded by an arbitrarily small constant. Here, we adopt an alternative approach where the (parameter-dependent) test space is defined by

$$W_r(\xi) = P_m(\xi)^* R_V V_r, \qquad (44)$$

where $P_m(\xi)$ is an interpolation of the inverse of $A(\xi)$ using $m$ interpolation points in the parameter set $\Xi$. In practice, when $A(\xi)$ is a matrix, algorithms developed in [22] can be used. This will be detailed later on. The underlying idea is to obtain a test space as close as possible to the ideal test space $A^{-*}(\xi)R_VV_r$ defined in (43). For $m = 0$, with $P_0(\xi) = R_V^{-1}$ by convention, we have $W_r = V_r$, which yields the standard Galerkin projection.

## 3.2 Error estimates for vector-valued variables of interest

In this section, we propose practical error estimates for the variable of interest, first for the primal-dual approach and then for the saddle point method.

### 3.2.1 Primal-dual approach

Given approximations $\tilde{u}$ and $\tilde{Q}$ of the primal solution $u$ and the dual solution $Q$ respectively, a standard approach is to start from the error bound

$$\|s(\xi) - \tilde{s}(\xi)\|_Z \le \|u(\xi) - \tilde{u}(\xi)\|_V \|L(\xi)^* - A(\xi)^* \tilde{Q}(\xi)\|_{Z' \to V'}$$

which is provided by Proposition 2.5. This suggests to measure the norm of the residuals associated to the primal and dual variables. In practice, we distinguish two cases.

In the case where the operator $A(\xi)$ is symmetric and coercive, it is natural to choose the parameter-dependent norm $\| \cdot \|_V$ as the one induced by the operator, i.e. $R_V = R_V(\xi) = A(\xi)$. However, neither the primal error $\|u(\xi) - \tilde{u}(\xi)\|_V = \langle b(\xi) - A(\xi)\tilde{u}(\xi), u(\xi) - \tilde{u}(\xi) \rangle$ nor the dual residual norm $\|L(\xi)^* - A(\xi)^* \tilde{Q}(\xi)\|_{Z' \to V'} = \sup_{\|z\|_Z=1} \|Q(\xi)z - \tilde{Q}(\xi)z\|_Z$ can be computed without computing the primal and dual solutions $u(\xi)$ and $Q(\xi)$. The classical way to circumvent this issue is to introduce a parameter-independent norm $\| \cdot \|_{V_0}$, which is in general the “natural” norm associated to the space $V$, and to measure residuals with the associated dual norm $\| \cdot \|_{V'_0}$. Here we assume that the operator $A(\xi)$ satisfies

$$\alpha(\xi) \|v\|_{V_0} \le \|A(\xi)v\|_{V'_0} \qquad (45)$$
for all $v \in V$, where $\alpha(\xi) > 0$. By definition of the norm $\| \cdot \|_V$, we can write

$$ \|v\|_V^2 = \langle A(\xi)v, v \rangle \le \|A(\xi)v\|_{V_0'} \|v\|_{V_0} \le \alpha(\xi)^{-1} \|A(\xi)v\|_{V_0'}^2 \quad \forall v \in V. $$

Then we have $\|u(\xi) - \tilde{u}(\xi)\|_V \le \alpha(\xi)^{-1/2} \|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'}$. In the same way, we can prove that $\|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V'} \le \alpha(\xi)^{-1/2} \|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V_0'}$. Finally, we obtain

$$ \|s(\xi) - \tilde{s}(\xi)\|_Z \le \frac{\|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'} \|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V_0'}}{\alpha(\xi)} := \Delta(\xi), \quad (46) $$

where $\Delta(\xi)$ is a certified error bound for the variable of interest, which involves computable primal and dual residual norms.

In the general case, we consider for $\| \cdot \|_V$ the natural norm on $V$, i.e. $\| \cdot \|_V = \| \cdot \|_{V_0}$. As a consequence, the norm of the dual residual is computable, but the computation of the error $\|u(\xi) - \tilde{u}(\xi)\|_{V_0}$ requires the primal solution $u(\xi)$ which is not available in practice. Once again, we assume that the operator satisfies the property (45) so that we can write $\|u(\xi) - \tilde{u}(\xi)\|_{V_0} \le \alpha(\xi)^{-1} \|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'}$. Then we end up with the same error bound (46) for the variable of interest.

### 3.2.2 Saddle point method

We now derive new error bounds in the case where the approximation $\tilde{s}(\xi)$ is obtained by the saddle point method introduced in Section 2.3. Let us start from the error bound

$$ \|s(\xi) - \tilde{s}(\xi)\|_Z \le \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^*z' - A(\xi)^*y\|_{V'}}{\|z'\|_{Z'}} \|u(\xi) - u_{r,p}(\xi) - R_V^{-1}A(\xi)^*y_{r,p}(\xi)\|_V \quad (47) $$

provided by Proposition 2.14. Once again, we distinguish two cases.

For the case where the operator $A(\xi)$ is symmetric and coercive, we consider for $\| \cdot \|_V$ the norm induced by the operator, i.e. $R_V = A$. According to Remark 2.15, the quantity $t_{r,p}(\xi) = u_{r,p}(\xi) - R_V^{-1}(\xi)A(\xi)^*y_{r,p}(\xi) = u_{r,p}(\xi) + y_{r,p}(\xi)$ is nothing but the orthogonal projection of $u(\xi)$ onto $T_p = W_r + W_k^Q$, with $W_r = V_r$. Then for any $\tilde{t}_{r,p} \in T_p$ we have

$$ \|u(\xi) - t_{r,p}(\xi)\|_V^2 \le \|u(\xi) - \tilde{t}_{r,p}\|_V^2 \le \alpha(\xi)^{-1} \|b(\xi) - A(\xi)\tilde{t}_{r,p}\|_{V_0'}^2, $$

where the norm $\| \cdot \|_{V_0}$ is the natural norm on $V$ such that (45) holds. Then, taking the infimum over $\tilde{t}_{r,p} \in T_p$ we obtain

$$ \|u(\xi) - t_{r,p}(\xi)\|_V \le \alpha(\xi)^{-1/2} \inf_{\tilde{t}_{r,p} \in T_p} \|A(\xi)\tilde{t}_{r,p} - b(\xi)\|_{V_0'}. $$
for all $v \in V$, where $\alpha(\xi) > 0$. By definition of the norm $\| \cdot \|_V$, we can write

$$ \|v\|_V^2 = \langle A(\xi)v, v \rangle \le \|A(\xi)v\|_{V_0'} \|v\|_{V_0} \le \alpha(\xi)^{-1} \|A(\xi)v\|_{V_0'}^2 \quad \forall v \in V. $$

Then we have $\|u(\xi) - \tilde{u}(\xi)\|_V \le \alpha(\xi)^{-1/2} \|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'}$. In the same way, we can prove that $\|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V'} \le \alpha(\xi)^{-1/2} \|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V_0'}$. Finally, we obtain

$$ \|s(\xi) - \tilde{s}(\xi)\|_Z \le \frac{\|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'} \|L(\xi)^* - A(\xi)^*\tilde{Q}(\xi)\|_{Z' \to V_0'}}{\alpha(\xi)} := \Delta(\xi), \quad (46) $$

where $\Delta(\xi)$ is a certified error bound for the variable of interest, which involves computable primal and dual residual norms.

In the general case, we consider for $\| \cdot \|_V$ the natural norm on $V$, i.e. $\| \cdot \|_V = \| \cdot \|_{V_0}$. As a consequence, the norm of the dual residual is computable, but the computation of the error $\|u(\xi) - \tilde{u}(\xi)\|_{V_0}$ requires the primal solution $u(\xi)$ which is not available in practice. Once again, we assume that the operator satisfies the property (45) so that we can write $\|u(\xi) - \tilde{u}(\xi)\|_{V_0} \le \alpha(\xi)^{-1} \|A(\xi)\tilde{u}(\xi) - b(\xi)\|_{V_0'}$. Then we end up with the same error bound (46) for the variable of interest.

### 3.2.2 Saddle point method

We now derive new error bounds in the case where the approximation $\tilde{s}(\xi)$ is obtained by the saddle point method introduced in Section 2.3. Let us start from the error bound

$$ \|s(\xi) - \tilde{s}(\xi)\|_Z \le \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^*z' - A(\xi)^*y\|_{V'}}{\|z'\|_{Z'}} \|u(\xi) - u_{r,p}(\xi) - R_V^{-1}A(\xi)^*y_{r,p}(\xi)\|_V \quad (47) $$

provided by Proposition 2.14. Once again, we distinguish two cases.

For the case where the operator $A(\xi)$ is symmetric and coercive, we consider for $\| \cdot \|_V$ the norm induced by the operator, i.e. $R_V = A$. According to Remark 2.15, the quantity $t_{r,p}(\xi) = u_{r,p}(\xi) - R_V^{-1}(\xi)A(\xi)^*y_{r,p}(\xi) = u_{r,p}(\xi) + y_{r,p}(\xi)$ is nothing but the orthogonal projection of $u(\xi)$ onto $T_p = W_r + W_k^Q$, with $W_r = V_r$. Then for any $\tilde{t}_{r,p} \in T_p$ we have

$$ \|u(\xi) - t_{r,p}(\xi)\|_V^2 \le \|u(\xi) - \tilde{t}_{r,p}\|_V^2 \le \alpha(\xi)^{-1} \|b(\xi) - A(\xi)\tilde{t}_{r,p}\|_{V_0'}^2, $$

where the norm $\| \cdot \|_{V_0}$ is the natural norm on $V$ such that (45) holds. Then, taking the infimum over $\tilde{t}_{r,p} \in T_p$ we obtain

$$ \|u(\xi) - t_{r,p}(\xi)\|_V \le \alpha(\xi)^{-1/2} \inf_{\tilde{t}_{r,p} \in T_p} \|A(\xi)\tilde{t}_{r,p} - b(\xi)\|_{V_0'}. $$
Finally, we obtain that

$$ \begin{aligned} & \|s(\xi) - \tilde{s}(\xi)\|_{Z} \\ & \le \frac{1}{\alpha(\xi)} \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^{*}z' - A(\xi)^{*}y\|_{V_0'}}{\|z'\|_{Z'}} \min_{\tilde{t}_{r,p} \in T_p} \|A(\xi)\tilde{t}_{r,p} - b(\xi)\|_{V_0'} := \Delta(\xi). \end{aligned} \quad (48) $$

Note that the main difference between this error estimate and the previous one (46) is the minimization problem over $T_p$ in both primal and dual residuals. The solution of those minimization problems lead to additional computational costs, but sharper error bounds will be obtained, as illustrated by the numerical examples in the next section.

For the general case, we consider $\|\cdot\|_V = \|\cdot\|_{V_0}$. Starting from (47) and using the relation (45) to bound the primal error by the primal residual norm, we obtain the following error estimate

$$ \begin{aligned} & \|s(\xi) - \tilde{s}(\xi)\|_{Z} \\ & \le \frac{1}{\alpha(\xi)} \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^{*}z' - A(\xi)^{*}y\|_{V_0'}}{\|z'\|_{Z'}} \|A(\xi)t_{r,p}(\xi) - b(\xi)\|_{V_0'} := \Delta(\xi), \end{aligned} \quad (49) $$

where $t_{r,p}(\xi) = u_{r,p}(\xi) + R_{V_0}^{-1}A(\xi)^{*}y_{r,p}(\xi)$.

**Remark 3.1.** All the proposed error estimates rely on the knowledge of $\alpha(\xi)$. In the case where $\alpha(\xi)$ can not be easily computed, we can replace it by a lower bound $\alpha^{LB}(\xi) \le \alpha(\xi)$, e.g. provided by a SCM procedure [13]. This option will not be considered here. Another option is to remove $\alpha(\xi)$ from the definitions of $\Delta(\xi)$, therefore leading to error estimates which are no more certified error bounds.

## 3.3 Greedy construction of the reduced spaces

Here, we propose different greedy algorithms for the construction of the reduced spaces $V_r$ and $W_k^Q$. At each iteration, we search for a parameter value $\xi^* \in \Xi$ where the error estimate $\Delta(\xi)$ is maximum, i.e.

$$ \xi^* \in \arg\max_{\xi \in \Xi} \Delta(\xi). \quad (50) $$

A first strategy is to simultaneously enrich both the primal approximation space

$$ V_{r+1} = V_r + \text{span}(u(\xi^*)) \quad (51) $$

and the dual approximation space

$$ W_{k+l}^Q = W_k^Q + \text{range}(Q(\xi^*)) \quad (52) $$
Finally, we obtain that

$$ \begin{aligned} & \|s(\xi) - \tilde{s}(\xi)\|_Z \\ & \le \frac{1}{\alpha(\xi)} \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^* z' - A(\xi)^* y\|_{V_0'}}{\|z'\|_{Z'}} \min_{\tilde{t}_{r,p} \in T_p} \|A(\xi) \tilde{t}_{r,p} - b(\xi)\|_{V_0'} := \Delta(\xi). \end{aligned} \quad (48) $$

Note that the main difference between this error estimate and the previous one (46) is the minimization problem over $T_p$ in both primal and dual residuals. The solution of those minimization problems lead to additional computational costs, but sharper error bounds will be obtained, as illustrated by the numerical examples in the next section.

For the general case, we consider $\|\cdot\|_V = \|\cdot\|_{V_0}$. Starting from (47) and using the relation (45) to bound the primal error by the primal residual norm, we obtain the following error estimate

$$ \begin{aligned} & \|s(\xi) - \tilde{s}(\xi)\|_Z \\ & \le \frac{1}{\alpha(\xi)} \sup_{0 \ne z' \in Z'} \min_{y \in T_p} \frac{\|L(\xi)^* z' - A(\xi)^* y\|_{V_0'}}{\|z'\|_{Z'}} \|A(\xi) t_{r,p}(\xi) - b(\xi)\|_{V_0'} := \Delta(\xi), \end{aligned} \quad (49) $$

where $t_{r,p}(\xi) = u_{r,p}(\xi) + R_{V_0}^{-1} A(\xi)^* y_{r,p}(\xi)$.

**Remark 3.1.** All the proposed error estimates rely on the knowledge of $\alpha(\xi)$. In the case where $\alpha(\xi)$ can not be easily computed, we can replace it by a lower bound $\alpha^{LB}(\xi) \le \alpha(\xi)$, e.g. provided by a SCM procedure [13]. This option will not be considered here. Another option is to remove $\alpha(\xi)$ from the definitions of $\Delta(\xi)$, therefore leading to error estimates which are no more certified error bounds.

## 3.3 Greedy construction of the reduced spaces

Here, we propose different greedy algorithms for the construction of the reduced spaces $V_r$ and $W_k^Q$. At each iteration, we search for a parameter value $\xi^* \in \Xi$ where the error estimate $\Delta(\xi)$ is maximum, i.e.

$$ \xi^* \in \arg\max_{\xi \in \Xi} \Delta(\xi). \quad (50) $$

A first strategy is to simultaneously enrich both the primal approximation space

$$ V_{r+1} = V_r + \operatorname{span}(u(\xi^*)) \quad (51) $$

and the dual approximation space

$$ W_{k+l}^Q = W_k^Q + \operatorname{range}(Q(\xi^*)) \quad (52) $$
at each iteration. This strategy is referred as the simultaneous construction, as opposed to the alternate construction which consists in enriching $W_k^Q$ (resp. $V_r$) if $V_r$ (resp. $W_k^Q$) were enriched at the previous greedy iteration step.

**Remark 3.2.** In the literature, and for scalar-valued variables of interest, the classical approaches are either a separated construction of $V_r$ and $W_k^Q$ (using two independent greedy algorithms, see for e.g. [10, 19]), or a simultaneous construction (see e.g. [18]). The latter option can take advantage of a single factorization of the operator $A(\xi^*)$ to compute both the primal and dual variables. The alternate construction proposed here is not usual. This possibility is mentioned in remark 2.47 of the tutorial [11].

For vector-valued variables of interest (dim($Z$) > 1), the enrichment strategy (52) makes sense only if $dim(\text{range}(Q(\xi^*)) < \infty$, in which case $l = \text{dim}(W_{k+l}^Q) - \text{dim}(W_l^Q) < \infty$. However, if $\text{dim}(\text{range}(Q(\xi^*))$ is finite but very high, the enrichment strategy (52) may lead to a rapid increase of the dimension of the dual approximation space. Therefore, when $\text{dim}(\text{range}(Q(\xi^*))$ is infinite or very high, we propose to replace the enrichment strategy (52) by

$$W_{k+1}^Q = W_k^Q + \text{span}(Q(\xi^*)z'), \quad (53)$$

where the space $W_k^Q$ is enriched with a single vector $Q(\xi^*)z'$, with $z' \in Z'$ such that

$$z' \in \arg\max_{\tilde{z}' \in Z'} \frac{\|(L(\xi^*)^* - A(\xi^*)^* \tilde{Q}(\xi))\tilde{z}'\|_{V_0'}}{\|\tilde{z}'\|_{Z'}} \quad (\text{for primal-dual method}), \text{ or} \quad (54)$$

$$z' \in \arg\max_{\tilde{z}' \in Z'} \min_{y \in W_k^Q} \frac{\|L(\xi^*)^*\tilde{z}' - A(\xi^*)^*y\|_{V_0'}}{\|\tilde{z}'\|_{Z'}} \quad (\text{for saddle point method}). \quad (55)$$

Contrarily to the full enrichment (52), this partial enrichment does not necessarily lead to a zero error at the point $\xi^*$ for the next iterations. Then we expect that (53) will deteriorate the convergence properties of the algorithm, but for $\text{dim}(Z) \gg 1$, the space $W_{k+1}^Q$ defined by (53) will have a much lower dimension than the space $W_{k+l}^Q$ defined by (52). It is worth mentioning that in [8], the authors propose a similar partial enrichment strategy for the test space $T_p$ but not in a goal-oriented framework.

The definition (44) of the test space $W_r$ requires the definition of a preconditioner $P_m(\xi)$ which is here constructed by interpolation of the inverse of $A(\xi)$. Following the idea of [22], the interpolation points for the preconditioner are chosen as the points where solutions (primal and dual) have already been computed, i.e. the points given by (50). The resulting algorithms are summarized in Algorithm 1 and Algorithm 2 respectively for the simultaneous and the alternate constructions of $V_r$ and $W_k^Q$.
at each iteration. This strategy is referred as the simultaneous construction, as opposed to the alternate construction which consists in enriching $W_k^Q$ (resp. $V_r$) if $V_r$ (resp. $W_k^Q$) were enriched at the previous greedy iteration step.

**Remark 3.2.** In the literature, and for scalar-valued variables of interest, the classical approaches are either a separated construction of $V_r$ and $W_k^Q$ (using two independent greedy algorithms, see for e.g. [10, 19]), or a simultaneous construction (see e.g. [18]). The latter option can take advantage of a single factorization of the operator $A(\xi^*)$ to compute both the primal and dual variables. The alternate construction proposed here is not usual. This possibility is mentioned in remark 2.47 of the tutorial [11].

For vector-valued variables of interest (dim($Z$) > 1), the enrichment strategy (52) makes sense only if $dim(\text{range}(Q(\xi^*)) < \infty$, in which case $l = \text{dim}(W_{k+l}^Q) - \text{dim}(W_l^Q) < \infty$. However, if $\text{dim}(\text{range}(Q(\xi^*))$ is finite but very high, the enrichment strategy (52) may lead to a rapid increase of the dimension of the dual approximation space. Therefore, when $\text{dim}(\text{range}(Q(\xi^*))$ is infinite or very high, we propose to replace the enrichment strategy (52) by

$$W_{k+1}^Q = W_k^Q + \text{span}(Q(\xi^*)z'), \qquad (53)$$

where the space $W_k^Q$ is enriched with a single vector $Q(\xi^*)z'$, with $z' \in Z'$ such that

$$z' \in \arg \max_{\tilde{z}' \in Z'} \frac{\| (L(\xi^*)^* - A(\xi^*)^* \tilde{Q}(\xi)) \tilde{z}' \|_{V_0'}}{\| \tilde{z}' \|_{Z'}} \quad (\text{for primal-dual method}), \text{ or} \quad (54)$$

$$z' \in \arg \max_{\tilde{z}' \in Z'} \min_{y \in W_k^Q} \frac{\| L(\xi^*)^* \tilde{z}' - A(\xi^*)^* y \|_{V_0'}}{\| \tilde{z}' \|_{Z'}} \quad (\text{for saddle point method}). \quad (55)$$

Contrarily to the full enrichment (52), this partial enrichment does not necessarily lead to a zero error at the point $\xi^*$ for the next iterations. Then we expect that (53) will deteriorate the convergence properties of the algorithm, but for $dim(Z) \gg 1$, the space $W_{k+1}^Q$ defined by (53) will have a much lower dimension than the space $W_{k+l}^Q$ defined by (52). It is worth mentioning that in [8], the authors propose a similar partial enrichment strategy for the test space $T_p$ but not in a goal-oriented framework.

The definition (44) of the test space $W_r$ requires the definition of a preconditioner $P_m(\xi)$ which is here constructed by interpolation of the inverse of $A(\xi)$. Following the idea of [22], the interpolation points for the preconditioner are chosen as the points where solutions (primal and dual) have already been computed, i.e. the points given by (50). The resulting algorithms are summarized in Algorithm 1 and Algorithm 2 respectively for the simultaneous and the alternate constructions of $V_r$ and $W_k^Q$.
**Algorithm 1** Simultaneous construction of $V_r$ and $W_k^Q$

**Input:** Error estimator $\Delta(\cdot)$, a samples set $\Xi$, maximum iteration $I$

1: Initialize $r, k = 0$ and the spaces $V_r = \{0\}$ and $W_k^Q = \{0\}$

2: **for** $i = 1$ to $I$ **do**

3: Find $\xi_i \in \arg\max_{\xi \in \Xi} \Delta(\xi)$

4: Compute a factorization of $A(\xi_i)$ and update the preconditioner if needed

5: Solve $u(\xi_i) = A(\xi_i)^{-1}b(\xi_i)$

6: Update $V_{r+1} = V_r + \text{span}(u(\xi_i))$, and $r \leftarrow r + 1$

7: **if** Full dual enrichment **then**

8: Solve $Q(\xi_i) = A(\xi_i)^{-*}L(\xi_i)^*$

9: Update $W_{k+l}^Q = W_k^Q + \text{range}(Q(\xi_i))$, and $k \leftarrow k + l$

10: **else if** Partial dual enrichment **then**

11: Find $z'$ according to (55) or (54)

12: Solve $y(\xi_i) = A(\xi_i)^{-*}(L(\xi_i)^*z')$

13: Update $W_{k+1}^Q = W_k^Q + \text{span}(y(\xi_i))$, and $k \leftarrow k + 1$

14: **end if**

15: **end for**

# 4 Numerical results

In this section, we present numerical applications of the methods proposed in Sections 2 and 3. We first describe the applications in Section 4.1. Then we compare the projection methods for the estimation of a variable of interest in Section 4.2. Finally, we study the behavior of the proposed greedy algorithms for the construction of the reduced spaces in Section 4.3.

## 4.1 Applications

### 4.1.1 Application 1 : a symmetric problem

We consider a linear elasticity problem<fnref n="2" /> $\text{div}(K(\xi) : \varepsilon(\mathbf{u}(\xi))) = 0$ over a domain $\Omega$ (represented in Figure 1(a)), where $\mathbf{u}(\xi) : \Omega \to \mathbb{R}^3$ is the displacement field and $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) \in \mathbb{R}^{3 \times 3}$ is the strain tensor associated to the displacement field $\mathbf{u}$. The Hooke tensor $K(\xi)$ is such that

$$K(\xi) : \varepsilon(\mathbf{u}(\xi)) = \frac{E(\xi)}{1 + \nu} \left( \varepsilon(\mathbf{u}(\xi)) + \frac{\nu}{1 - 2\nu} \text{trace}(\varepsilon(\mathbf{u}(\xi)))I_3 \right),$$

<footnote n="2">²The authors thank Mathilde Chevreuil for having proposed this benchmark problem.</footnote>
**Algorithm 1** Simultaneous construction of $V_r$ and $W_k^Q$

**Input:** Error estimator $\Delta(\cdot)$, a samples set $\Xi$, maximum iteration $I$

1: Initialize $r, k = 0$ and the spaces $V_r = \{0\}$ and $W_k^Q = \{0\}$

2: **for** $i = 1$ to $I$ **do**

3: Find $\xi_i \in \arg\max_{\xi \in \Xi} \Delta(\xi)$

4: Compute a factorization of $A(\xi_i)$ and update the preconditioner if needed

5: Solve $u(\xi_i) = A(\xi_i)^{-1}b(\xi_i)$

6: Update $V_{r+1} = V_r + \text{span}(u(\xi_i))$, and $r \leftarrow r + 1$

7: **if** Full dual enrichment **then**

8: Solve $Q(\xi_i) = A(\xi_i)^{-*}L(\xi_i)^*$

9: Update $W_{k+l}^Q = W_k^Q + \text{range}(Q(\xi_i))$, and $k \leftarrow k + l$

10: **else if** Partial dual enrichment **then**

11: Find $z'$ according to (55) or (54)

12: Solve $y(\xi_i) = A(\xi_i)^{-*}(L(\xi_i)^*z')$

13: Update $W_{k+1}^Q = W_k^Q + \text{span}(y(\xi_i))$, and $k \leftarrow k + 1$

14: **end if**

15: **end for**

# 4 Numerical results

In this section, we present numerical applications of the methods proposed in Sections 2 and 3. We first describe the applications in Section 4.1. Then we compare the projection methods for the estimation of a variable of interest in Section 4.2. Finally, we study the behavior of the proposed greedy algorithms for the construction of the reduced spaces in Section 4.3.

## 4.1 Applications

### 4.1.1 Application 1 : a symmetric problem

We consider a linear elasticity problem<fnref n="2" /> $\text{div}(K(\xi) : \varepsilon(\mathbf{u}(\xi))) = 0$ over a domain $\Omega$ (represented in Figure 1(a)), where $\mathbf{u}(\xi) : \Omega \to \mathbb{R}^3$ is the displacement field and $\varepsilon(\mathbf{u}) = \frac{1}{2}(\nabla \mathbf{u} + \nabla \mathbf{u}^T) \in \mathbb{R}^{3 \times 3}$ is the strain tensor associated to the displacement field $\mathbf{u}$. The Hooke tensor $K(\xi)$ is such that

$$K(\xi) : \varepsilon(\mathbf{u}(\xi)) = \frac{E(\xi)}{1 + \nu} \left( \varepsilon(\mathbf{u}(\xi)) + \frac{\nu}{1 - 2\nu} \text{trace}(\varepsilon(\mathbf{u}(\xi)))I_3 \right),$$

<footnote n="2">²The authors thank Mathilde Chevreuil for having proposed this benchmark problem.</footnote>
**Algorithm 2** Alternate construction of $V_r$ and $W_k^Q$

**Input:** Error estimator $\Delta(\cdot)$, a samples set $\Xi$, maximum iteration $I$

1: Initialize $r, k = 0$ and the spaces $V_r = \{0\}$ and $W_k^Q = \{0\}$

2: **for** $i = 1$ to $I$ **do**

3: Find $\xi_i \in \arg\max_{\xi \in \Xi} \Delta(\xi)$

4: Compute a factorization of $A(\xi_i)$ and update the preconditioner if needed

5: **if** $i$ is even **then**

6: Solve $u(\xi_i) = A(\xi_i)^{-1}b(\xi_i)$

7: Update $V_{r+1} = V_r + \text{span}(u(\xi_i))$, and $r \leftarrow r + 1$

8: **else if** $i$ is odd **then**

9: **if** Full dual enrichment **then**

10: Solve $Q(\xi_i) = A(\xi_i)^{-*}L(\xi_i)^*$

11: Update $W_{k+l}^Q = W_k^Q + \text{range}(Q(\xi_i))$, and $k \leftarrow k + l$

12: **else if** Partial dual enrichment **then**

13: Find $z'$ according to (55) or (54)

14: Solve $y(\xi_i) = A(\xi_i)^{-*}(L(\xi_i)^*z')$

15: Update $W_{k+1}^Q = W_k^Q + \text{span}(y(\xi_i))$, and $k \leftarrow k + 1$

16: **end if**

17: **end if**

18: **end for**

where $\nu = 0.3$ is the Poisson coefficient and $E(\xi)$ is the Young modulus defined by $E(\xi) = 1_{\Omega_0} + \sum_{i=1}^6 \xi_i 1_{\Omega_i}$, $1_{\Omega_i}$ being the indicator function of the subdomain $\Omega_i$, see Figure 1(b). The components of $\xi = (\xi_1, \dots, \xi_6)$ are independent and log-uniformly distributed over $[10^{-1}, 10]$. We impose homogeneous Dirichlet boundary condition $u(\xi) = 0$ on $\Gamma_D$ (red lines), a unit vertical surface load on $\Gamma_L$ (green square), and a zero surface load on the complementary part of the boundary (see Figure 1(a)). We consider the Galerkin approximation $u^h(\xi)$ of $u(\xi)$ on a $\mathbb{P}_1$ finite element approximation space $V^h = \text{span}(\phi_i)_{i=1}^n \subset \{\mathbf{v} \in H^1(\Omega)^3 : \mathbf{v}|_{\Gamma_D} = 0\}$ of dimension $n = 8916$ associated to the mesh plotted in Figure 1(b). The vector $u(\xi) \in V = \mathbb{R}^n$ such that $u^h(\xi) = \sum_{i=1}^n u_i(\xi)\phi_i$ is the solution of the linear system $A(\xi)u(\xi) = b$ of size $n$, with

$$A(\xi) = A^{(0)} + \sum_{k=1}^{6} \xi_i A^{(k)}, \quad A_{i,j}^{(k)} = \int_{\Omega_k} \nabla \phi_i : K_0 : \nabla \phi_j \, d\Omega, \quad (56)$$

and $b_i = \int_{\Gamma_L} -\mathbf{e}_3 \cdot \phi_i \, d\Gamma$, where $K_0$ denotes the Hooke tensor with the Young modulus $E = 1$. The norm $\| \cdot \|_V$ on the space $V$ is chosen such that $\|v\|_V^2 = \langle A(\xi)v, v \rangle$, that means $R_V = A(\xi)$. We also consider the parameter-independent norm $\| \cdot \|_{V_0}$ defined by
**Algorithm 2** Alternate construction of $V_r$ and $W_k^Q$

**Input:** Error estimator $\Delta(\cdot)$, a samples set $\Xi$, maximum iteration $I$

1: Initialize $r, k = 0$ and the spaces $V_r = \{0\}$ and $W_k^Q = \{0\}$

2: **for** $i = 1$ to $I$ **do**

3: Find $\xi_i \in \arg\max_{\xi \in \Xi} \Delta(\xi)$

4: Compute a factorization of $A(\xi_i)$ and update the preconditioner if needed

5: **if** $i$ is even **then**

6: Solve $u(\xi_i) = A(\xi_i)^{-1}b(\xi_i)$

7: Update $V_{r+1} = V_r + \text{span}(u(\xi_i))$, and $r \leftarrow r + 1$

8: **else if** $i$ is odd **then**

9: **if** Full dual enrichment **then**

10: Solve $Q(\xi_i) = A(\xi_i)^{-*}L(\xi_i)^*$

11: Update $W_{k+l}^Q = W_k^Q + \text{range}(Q(\xi_i))$, and $k \leftarrow k + l$

12: **else if** Partial dual enrichment **then**

13: Find $z'$ according to (55) or (54)

14: Solve $y(\xi_i) = A(\xi_i)^{-*}(L(\xi_i)^*z')$

15: Update $W_{k+1}^Q = W_k^Q + \text{span}(y(\xi_i))$, and $k \leftarrow k + 1$

16: **end if**

17: **end if**

18: **end for**

where $\nu = 0.3$ is the Poisson coefficient and $E(\xi)$ is the Young modulus defined by $E(\xi) = 1_{\Omega_0} + \sum_{i=1}^6 \xi_i 1_{\Omega_i}$, $1_{\Omega_i}$ being the indicator function of the subdomain $\Omega_i$, see Figure 1(b). The components of $\xi = (\xi_1, \dots, \xi_6)$ are independent and log-uniformly distributed over $[10^{-1}, 10]$. We impose homogeneous Dirichlet boundary condition $u(\xi) = 0$ on $\Gamma_D$ (red lines), a unit vertical surface load on $\Gamma_L$ (green square), and a zero surface load on the complementary part of the boundary (see Figure 1(a)). We consider the Galerkin approximation $u^h(\xi)$ of $u(\xi)$ on a $\mathbb{P}_1$ finite element approximation space $V^h = \text{span}(\phi_i)_{i=1}^n \subset \{\mathbf{v} \in H^1(\Omega)^3 : \mathbf{v}|_{\Gamma_D} = 0\}$ of dimension $n = 8916$ associated to the mesh plotted in Figure 1(b). The vector $u(\xi) \in V = \mathbb{R}^n$ such that $u^h(\xi) = \sum_{i=1}^n u_i(\xi)\phi_i$ is the solution of the linear system $A(\xi)u(\xi) = b$ of size $n$, with

$$A(\xi) = A^{(0)} + \sum_{k=1}^{6} \xi_i A^{(k)}, \quad A_{i,j}^{(k)} = \int_{\Omega_k} \nabla \phi_i : K_0 : \nabla \phi_j \, d\Omega, \quad (56)$$

and $b_i = \int_{\Gamma_L} -\mathbf{e}_3 \cdot \phi_i \, d\Gamma$, where $K_0$ denotes the Hooke tensor with the Young modulus $E = 1$. The norm $\| \cdot \|_V$ on the space $V$ is chosen such that $\|v\|_V^2 = \langle A(\xi)v, v \rangle$, that means $R_V = A(\xi)$. We also consider the parameter-independent norm $\| \cdot \|_{V_0}$ defined by
$\|v\|_{V_0}^2 = \langle A(\xi_0)v, v \rangle$ with $\xi_0 = (1, \dots, 1)$. It corresponds to the norm induced by the operator associated with the Hooke tensor $K_0$ instead of $K(\xi)$.

Let us consider $s^h(\xi) = \mathbf{u}_{\Gamma}^h(\xi) \cdot \mathbf{e}_3$ which is the vertical displacement of the Galerkin approximation on the blue line $\Gamma$, see Figure 1(a). We can write $s^h(\xi) = \sum_{j=1}^l s_j(\xi)\psi_j$ where $\{\psi_j\}_{j=1}^l$ is a basis of the space $\{\mathbf{v}_{\Gamma}^h \cdot \mathbf{e}_3 : \mathbf{v}^h \in \mathbf{V}^h\}$ of dimension $l = 44$. Then there exists $L \in \mathbb{R}^{l \times n}$ such that

$$s(\xi) = Lu(\xi),$$

where $s(\xi) = (s_1(\xi), \dots, s_l(\xi)) \in Z = \mathbb{R}^l$ is the variable of interest. The norm $\|\cdot\|_Z$ is defined as the canonical norm of $\mathbb{R}^l$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Application 1: schematic representation of the problem and a realization of the solution.</figcaption></figure>
$\|v\|_{V_0}^2 = \langle A(\xi_0)v, v \rangle$ with $\xi_0 = (1, \dots, 1)$. It corresponds to the norm induced by the operator associated with the Hooke tensor $K_0$ instead of $K(\xi)$.

Let us consider $s^h(\xi) = \mathbf{u}_{\Gamma}^h(\xi) \cdot \mathbf{e}_3$ which is the vertical displacement of the Galerkin approximation on the blue line $\Gamma$, see Figure 1(a). We can write $s^h(\xi) = \sum_{j=1}^l s_j(\xi)\psi_j$ where $\{\psi_j\}_{j=1}^l$ is a basis of the space $\{\mathbf{v}_{\Gamma}^h \cdot \mathbf{e}_3 : \mathbf{v}^h \in \mathbf{V}^h\}$ of dimension $l = 44$. Then there exists $L \in \mathbb{R}^{l \times n}$ such that

$$s(\xi) = Lu(\xi),$$

where $s(\xi) = (s_1(\xi), \dots, s_l(\xi)) \in Z = \mathbb{R}^l$ is the variable of interest. The norm $\|\cdot\|_Z$ is defined as the canonical norm of $\mathbb{R}^l$.

<figure><img src="image_3.png" /><figcaption>Figure 1: Application 1: schematic representation of the problem and a realization of the solution.</figcaption></figure>
### 4.1.2 Application 2: a non symmetric problem

We consider the benchmark problem of the cooling of electronic components proposed in the OPUS project<fnref n="3" />. The equation to solve is an advection-diffusion equation over the domain $\Omega \subset \mathbb{R}^2$

$$ - \nabla \cdot (\kappa(\xi)\nabla T(\xi)) + \mathbf{c}(\xi) \cdot \nabla T(\xi) = f, \qquad (57) $$

whose solution $T(\xi) : \Omega \to \mathbb{R}$ is the temperature field. Here $\kappa$ and $\mathbf{c}$ denote respectively the diffusion coefficient and the advection field, which are parameter-dependent coefficients of the operator. The full description of this problem is given in [22]. Here, we only focus on the resulting algebraic parameter-dependent equation coming from stabilized finite element discretization of (57), that is $A(\xi)u(\xi) = b(\xi)$, where $u(\xi) \in \mathbb{R}^n$ are the coefficients of the finite element approximation $T^h = \sum_{i=1}^n u_i\varphi_i$ of $T$, and where $\xi = (\xi_1, \dots, \xi_4)$ is a 4-dimensional random vector. The space $V = \mathbb{R}^n$ with $n = 2.8 \times 10^4$ is endowed with the norm $\| \cdot \|_V = \| \cdot \|_{V_0}$ which corresponds to the $H^1(\Omega)$-norm<fnref n="4" />. The variable of interest $s(\xi) = (s_1(\xi), s_2(\xi))$ is the mean temperature of both electronic components, with

$$ s_1(\xi) = \frac{1}{|\Omega_{IC_1}|} \int_{\Omega_{IC_1}} T^h(\xi) d\Omega, \quad s_2(\xi) = \frac{1}{|\Omega_{IC_2}|} \int_{\Omega_{IC_2}} T^h(\xi) d\Omega, \qquad (58) $$

where $\Omega_{IC_i}$ ($i = 1, 2$) are two subdomains of $\Omega \subset \mathbb{R}^2$ (see [22, Fig.7]). Then we can write $s(\xi) = Lu(\xi)$ for an appropriate $L \in \mathbb{R}^{l \times n}$, with $l = 2$. Here we have $Z = \mathbb{R}^2$, which we equip with the canonical norm on $\mathbb{R}^2$.

## 4.2 Comparison of the projections methods

The goal of this section is to compare the projection methods proposed in Section 2 for the estimation of $s(\xi)$. Here the approximation spaces $V_r, W_k^Q$ and the test space $W_r$ are given. We denote by $V_r, W_k^Q$ and $W_r$ the matrices containing the basis vectors of the corresponding subspaces. In order to improve condition numbers of reduced systems of equations, these bases are orthogonalized using a Gram-Schmidt procedure.

### 4.2.1 Application 1

We first detail how we build $V_r, W_k^Q$ and $W_r$. The matrix $V_r$ contains $r = 20$ snapshots of the solution: $V_r = (u(\xi_1), \dots, u(\xi_{20}))$. The test space is $W_r = V_r$, which corresponds to a standard Galerkin projection method. The matrix $W_k^Q$ contains 2 snapshots of the dual

<footnote n="3">³See http://www.opus-project.fr</footnote>

<footnote n="4">⁴ It means that $\|v\|_{V_0} = \|v^h\|_{H^1(\Omega)}$ for all $v \in V$, where $v^h = \sum_{i=1}^n v_i\varphi_i$.</footnote>
### 4.1.2 Application 2: a non symmetric problem

We consider the benchmark problem of the cooling of electronic components proposed in the OPUS project<fnref n="3" />. The equation to solve is an advection-diffusion equation over the domain $\Omega \subset \mathbb{R}^2$

$$ - \nabla \cdot (\kappa(\xi)\nabla T(\xi)) + \mathbf{c}(\xi) \cdot \nabla T(\xi) = f, \qquad (57) $$

whose solution $T(\xi) : \Omega \to \mathbb{R}$ is the temperature field. Here $\kappa$ and $\mathbf{c}$ denote respectively the diffusion coefficient and the advection field, which are parameter-dependent coefficients of the operator. The full description of this problem is given in [22]. Here, we only focus on the resulting algebraic parameter-dependent equation coming from stabilized finite element discretization of (57), that is $A(\xi)u(\xi) = b(\xi)$, where $u(\xi) \in \mathbb{R}^n$ are the coefficients of the finite element approximation $T^h = \sum_{i=1}^n u_i\varphi_i$ of $T$, and where $\xi = (\xi_1, \dots, \xi_4)$ is a 4-dimensional random vector. The space $V = \mathbb{R}^n$ with $n = 2.8 \times 10^4$ is endowed with the norm $\| \cdot \|_V = \| \cdot \|_{V_0}$ which corresponds to the $H^1(\Omega)$-norm<fnref n="4" />. The variable of interest $s(\xi) = (s_1(\xi), s_2(\xi))$ is the mean temperature of both electronic components, with

$$ s_1(\xi) = \frac{1}{|\Omega_{IC_1}|} \int_{\Omega_{IC_1}} T^h(\xi) d\Omega, \quad s_2(\xi) = \frac{1}{|\Omega_{IC_2}|} \int_{\Omega_{IC_2}} T^h(\xi) d\Omega, \qquad (58) $$

where $\Omega_{IC_i}$ ($i = 1, 2$) are two subdomains of $\Omega \subset \mathbb{R}^2$ (see [22, Fig.7]). Then we can write $s(\xi) = Lu(\xi)$ for an appropriate $L \in \mathbb{R}^{l \times n}$, with $l = 2$. Here we have $Z = \mathbb{R}^2$, which we equip with the canonical norm on $\mathbb{R}^2$.

## 4.2 Comparison of the projections methods

The goal of this section is to compare the projection methods proposed in Section 2 for the estimation of $s(\xi)$. Here the approximation spaces $V_r, W_k^Q$ and the test space $W_r$ are given. We denote by $V_r, W_k^Q$ and $W_r$ the matrices containing the basis vectors of the corresponding subspaces. In order to improve condition numbers of reduced systems of equations, these bases are orthogonalized using a Gram-Schmidt procedure.

### 4.2.1 Application 1

We first detail how we build $V_r, W_k^Q$ and $W_r$. The matrix $V_r$ contains $r = 20$ snapshots of the solution: $V_r = (u(\xi_1), \dots, u(\xi_{20}))$. The test space is $W_r = V_r$, which corresponds to a standard Galerkin projection method. The matrix $W_k^Q$ contains 2 snapshots of the dual

<footnote n="3">³See http://www.opus-project.fr</footnote>

<footnote n="4">⁴ It means that $\|v\|_{V_0} = \|v^h\|_{H^1(\Omega)}$ for all $v \in V$, where $v^h = \sum_{i=1}^n v_i \varphi_i$.</footnote>
variable $Q(\xi) = A(\xi)^{-1}L^* \in \mathbb{R}^{n \times l}$. Then $k = 2l = 88$. Finally, according to (42) the matrix $\mathbf{T}_p = (\mathbf{W}_r, \mathbf{W}_k^Q)$ is the concatenation of the matrices $\mathbf{W}_r$ and $\mathbf{W}_k^Q$.

We consider a samples set $\Xi_t \subset \Xi$ of size $t = 10^4$. For each $\xi \in \Xi_t$ we compute the exact quantity of interest $s(\xi)$ and the approximation $\tilde{s}(\xi)$ by the following methods.

* **Primal only:** solve the linear system $(\mathbf{V}_r^T A(\xi) \mathbf{V}_r) U_r(\xi) = (\mathbf{V}_r^T b)$ of size $r$ and compute $\tilde{s}(\xi) = (L \mathbf{V}_r) U_r(\xi)$.

* **Dual only:** solve the linear system $((\mathbf{W}_k^Q)^T A(\xi) \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b)$ of size $k$ and compute $\tilde{s}(\xi) = (L \mathbf{W}_k^Q) Y_k(\xi)^5$.

* **Primal-dual:** solve the linear system of the *Primal only* method, solve the linear system $((\mathbf{W}_k^Q)^T A(\xi) \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b) - ((\mathbf{W}_k^Q)^T A(\xi) \mathbf{V}_r) U_r(\xi)$ of size $k$ and compute $\tilde{s}(\xi) = (L \mathbf{V}_r) U_r(\xi) + (L \mathbf{W}_k^Q) Y_k(\xi)$.

* **Saddle point:** According to Remark 2.15, solve the linear system

$$ (\mathbf{T}_p^T A(\xi) \mathbf{T}_p) Y_p(\xi) = (\mathbf{T}_p^T b) $$

of size $p = k + r$, and compute $\tilde{s}(\xi) = (L \mathbf{T}_p) Y_p(\xi)$.

The affine decomposition (56) of matrix $A(\xi)$ allows for a rapid solution of the reduced systems for any parameter $\xi$.

Figure 2 gives the probability density function (PDF), the $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated over the samples set $\Xi_t$. We see that the *primal-dual* method provides errors for the quantity of interest which correspond to the product of the errors of the *primal only* and *dual only* methods. This reflects the “squared effect”. Moreover the *saddle point* method provides errors that are almost 10 times lower than the *primal-dual* method. This impressive improvement can be explained by the fact that the proposed problem is “almost compliant”, in the sense that the primal and dual solutions are similar: the primal solution is associated to a vertical force on the green square of Figure 1(a), and the dual solution is associated to a vertical loading on $\Gamma$. To illustrate this, let us consider a “less compliant” application where the variable of interest is defined as the horizontal displacement (in the direction $e_2$, see figure 1(a)) of the solution on the blue line $\Gamma$, i.e. $s^h(\xi) = \mathbf{u}_{|\Gamma}(\xi) \cdot \mathbf{e}_2$ (instead of $s^h(\xi) = \mathbf{u}_{|\Gamma}(\xi) \cdot \mathbf{e}_3$). The results are given in Figure 3. For this new setting, we can draw similar conclusions but the *saddle point* method provides a solution which is “only” 2 times better (instead of 10 times) than the *primal-dual* method.

<footnote n="5">⁵The dual only method corresponds to the *primal-dual* method where we consider a zero primal approximation, i.e. $V_r = W_r = \{0\}$.</footnote>
variable $Q(\xi) = A(\xi)^{-1}L^* \in \mathbb{R}^{n \times l}$. Then $k = 2l = 88$. Finally, according to (42) the matrix $\mathbf{T}_p = (\mathbf{W}_r, \mathbf{W}_k^Q)$ is the concatenation of the matrices $\mathbf{W}_r$ and $\mathbf{W}_k^Q$.

We consider a samples set $\Xi_t \subset \Xi$ of size $t = 10^4$. For each $\xi \in \Xi_t$ we compute the exact quantity of interest $s(\xi)$ and the approximation $\tilde{s}(\xi)$ by the following methods.

* **Primal only:** solve the linear system $(\mathbf{V}_r^T A(\xi) \mathbf{V}_r) U_r(\xi) = (\mathbf{V}_r^T b)$ of size $r$ and compute $\tilde{s}(\xi) = (L \mathbf{V}_r) U_r(\xi)$.

* **Dual only:** solve the linear system $((\mathbf{W}_k^Q)^T A(\xi) \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b)$ of size $k$ and compute $\tilde{s}(\xi) = (L \mathbf{W}_k^Q) Y_k(\xi)^5$.

* **Primal-dual:** solve the linear system of the *Primal only* method, solve the linear system $((\mathbf{W}_k^Q)^T A(\xi) \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b) - ((\mathbf{W}_k^Q)^T A(\xi) \mathbf{V}_r) U_r(\xi)$ of size $k$ and compute $\tilde{s}(\xi) = (L \mathbf{V}_r) U_r(\xi) + (L \mathbf{W}_k^Q) Y_k(\xi)$.

* **Saddle point:** According to Remark 2.15, solve the linear system

$$ (\mathbf{T}_p^T A(\xi) \mathbf{T}_p) Y_p(\xi) = (\mathbf{T}_p^T b) $$

of size $p = k + r$, and compute $\tilde{s}(\xi) = (L \mathbf{T}_p) Y_p(\xi)$.

The affine decomposition (56) of matrix $A(\xi)$ allows for a rapid solution of the reduced systems for any parameter $\xi$.

Figure 2 gives the probability density function (PDF), the $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated over the samples set $\Xi_t$. We see that the *primal-dual* method provides errors for the quantity of interest which correspond to the product of the errors of the *primal only* and *dual only* methods. This reflects the “squared effect”. Moreover the *saddle point* method provides errors that are almost 10 times lower than the *primal-dual* method. This impressive improvement can be explained by the fact that the proposed problem is “almost compliant”, in the sense that the primal and dual solutions are similar: the primal solution is associated to a vertical force on the green square of Figure 1(a), and the dual solution is associated to a vertical loading on $\Gamma$. To illustrate this, let us consider a “less compliant” application where the variable of interest is defined as the horizontal displacement (in the direction $e_2$, see figure 1(a)) of the solution on the blue line $\Gamma$, i.e. $s^h(\xi) = \mathbf{u}_{|\Gamma}(\xi) \cdot \mathbf{e}_2$ (instead of $s^h(\xi) = \mathbf{u}_{|\Gamma}(\xi) \cdot \mathbf{e}_3$). The results are given in Figure 3. For this new setting, we can draw similar conclusions but the *saddle point* method provides a solution which is “only” 2 times better (instead of 10 times) than the *primal-dual* method.

<footnote n="5">⁵The dual only method corresponds to the *primal-dual* method where we consider a zero primal approximation, i.e. $V_r = W_r = \{0\}$.</footnote>
<figure><img src="image_4.png" /><figcaption>Figure 2: Application 1: Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Application 1 with a different variable of interest (“less compliant” case): Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$.</figcaption></figure>

Now we consider the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ associated to the primal-dual error estimate defined by (46) and to the saddle-point error estimate defined by (48). For the considered application, the coercivity constant $\alpha(\xi)$ can be obtained by the min-theta method [11, Proposition 2.35]. Figure 4 presents statistical information on $\eta(\xi)$: the PDF, the mean, the max-min ratio and the normalized standard deviation estimated on a samples set of size $10^4$. We first observe in Figure 4(a) that the effectivity index is always greater than 1: this illustrates the fact that the error estimates are certified. Moreover, the error estimate of the saddle point method is much better than the one of the primal-dual method. The max-min ratio and the standard deviation of the corresponding effectivity index
<figure><img src="image_4.png" /><figcaption>Figure 2: Application 1: Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Figure 3: Application 1 with a different variable of interest (“less compliant” case): Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$.</figcaption></figure>

Now we consider the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ associated to the primal-dual error estimate defined by (46) and to the saddle-point error estimate defined by (48). For the considered application, the coercivity constant $\alpha(\xi)$ can be obtained by the *min-theta* method [11, Proposition 2.35]. Figure 4 presents statistical information on $\eta(\xi)$: the PDF, the mean, the max-min ratio and the normalized standard deviation estimated on a samples set of size $10^4$. We first observe in Figure 4(a) that the effectivity index is always greater than 1: this illustrates the fact that the error estimates are certified. Moreover, the error estimate of the saddle point method is much better than the one of the primal-dual method. The max-min ratio and the standard deviation of the corresponding effectivity index
are much smaller and the mean value is much closer to one for the saddle point method.

<figure><img src="image_6.png" /><figcaption>Figure 4: Application 1: Probability density function, mean, min-max ratio and normalized standard deviation of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi)-\tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$. Here, $\Delta(\xi)$ is defined by (46) for the primal-dual method and by (48) for the saddle point method.</figcaption></figure>

### 4.2.2 Application 2

For this second application, $V_r = (u(\xi_1), \dots, u(\xi_{50}))$ contains 50 snapshots of the primal solution ($r = 50$), and $W_k^Q = (Q(\xi_1) \dots Q(\xi_{25}))$ contains 25 snapshots of the dual solution so that the dimension of $W_k^Q$ is $k = 25l = 50$. The test space $W_r$ is defined according to (44), where $P_m(\xi)$ is an interpolation of $A(\xi)^{-1}$ using $m$ interpolation points selected by a greedy procedure based on the residual $\|I - P_m(\xi)A(\xi)\|_F$ (where $\|\cdot\|_F$ denotes the matrix Frobenius norm), see [22]. The interpolation is defined by a Frobenius semi-norm projection (with positivity constraint) using a random matrix with 400 columns. The matrix associated to the test space is given by $W_r(\xi) = P_m^T(\xi)R_VV_r$.

Once again, we consider a samples set $\Xi_t$ of size $t = 10^4$. For any $\xi \in \Xi_t$ we compute the exact quantity of interest $s(\xi)$ and the approximation $\tilde{s}(\xi)$ by the following methods.

* **Primal only:** solve the linear system $(W_r^T(\xi)A(\xi)V_r)U_r(\xi) = W_r(\xi)^Tb$ of size $r$ and compute $\tilde{s}(\xi) = (LV_r)U_r(\xi)$.

* **Dual only:** solve the linear system

$$((W_k^Q)^T A(\xi) R_V^{-1} A(\xi)^* W_k^Q) Y_k(\xi) = (W_k^Q)^T b$$
are much smaller and the mean value is much closer to one for the saddle point method.

<figure><img src="image_6.png" /><figcaption>Figure 4: Application 1: Probability density function, mean, min-max ratio and normalized standard deviation of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi)-\tilde{s}(\xi)\|_Z$ estimated on a samples set of size $10^4$. Here, $\Delta(\xi)$ is defined by (46) for the primal-dual method and by (48) for the saddle point method.</figcaption></figure>

### 4.2.2 Application 2

For this second application, $V_r = (u(\xi_1), \dots, u(\xi_{50}))$ contains 50 snapshots of the primal solution ($r = 50$), and $W_k^Q = (Q(\xi_1) \dots Q(\xi_{25}))$ contains 25 snapshots of the dual solution so that the dimension of $W_k^Q$ is $k = 25l = 50$. The test space $W_r$ is defined according to (44), where $P_m(\xi)$ is an interpolation of $A(\xi)^{-1}$ using $m$ interpolation points selected by a greedy procedure based on the residual $\|I - P_m(\xi)A(\xi)\|_F$ (where $\|\cdot\|_F$ denotes the matrix Frobenius norm), see [22]. The interpolation is defined by a Frobenius semi-norm projection (with positivity constraint) using a random matrix with 400 columns. The matrix associated to the test space is given by $W_r(\xi) = P_m^T(\xi)R_VV_r$.

Once again, we consider a samples set $\Xi_t$ of size $t = 10^4$. For any $\xi \in \Xi_t$ we compute the exact quantity of interest $s(\xi)$ and the approximation $\tilde{s}(\xi)$ by the following methods.

* Primal only: solve the linear system $(W_r^T(\xi)A(\xi)V_r)U_r(\xi) = W_r(\xi)^Tb$ of size $r$ and compute $\tilde{s}(\xi) = (LV_r)U_r(\xi)$.

* Dual only: solve the linear system

$$((W_k^Q)^T A(\xi) R_V^{-1} A(\xi)^* W_k^Q) Y_k(\xi) = (W_k^Q)^T b$$
of size k and compute $\tilde{s}(\xi) = (LR_V^{-1}A(\xi)^*\mathbf{W}_k^Q)Y_k(\xi)$.

* **Primal-dual:** solve the linear system of the Primal only method, solve the linear system

$$((\mathbf{W}_k^Q)^T A(\xi) R_V^{-1} A(\xi)^* \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b) - ((\mathbf{W}_k^Q)^T A(\xi) \mathbf{V}_r) U_r(\xi)$$

of size k, and compute $\tilde{s}(\xi) = (LV_r)U_r(\xi) + (LR_V^{-1}A(\xi)^*\mathbf{W}_k^Q)Y_k(\xi)$.

* **Saddle point:** solve the linear system of size p + r

$$\begin{pmatrix} \mathbf{T}_p^T(\xi) A(\xi) R_V^{-1} A(\xi)^* \mathbf{T}_p(\xi) & \mathbf{T}_p^T(\xi) A(\xi) \mathbf{V}_r \\ (\mathbf{T}_p^T(\xi) A(\xi) \mathbf{V}_r)^T & 0 \end{pmatrix} \begin{pmatrix} Y_{r,p}(\xi) \\ U_{r,p}(\xi) \end{pmatrix} = \begin{pmatrix} \mathbf{T}_p(\xi)^T b \\ 0 \end{pmatrix}$$

with $\mathbf{T}_p(\xi) = (\mathbf{W}_r(\xi), \mathbf{W}_k^Q)$, and compute

$$\tilde{s}(\xi) = (LV_r)U_{r,p}(\xi) + (LR_V^{-1}A(\xi)^*\mathbf{T}_p(\xi))Y_{r,p}(\xi).$$

The numerical results are given in Figure 5. Once again, the saddle point method leads to the lowest error on the variable of interest. Also, we see that a good preconditioner (for example with $m = 30$) improves the accuracy for the saddle point method, the primal only method and the primal-dual method. However, this improvement is not really significant for the considered application: the errors are barely divided by 2 compared to the non preconditioned Galerkin projection ($m = 0$). In fact, the preconditioner improves the quality of the test space, and the choice $W_r = V_r$ (yielding the standard Galerkin projection) is sufficiently accurate for this example and for the chosen norm on $V$.

We discuss now the quality of the error estimate $\Delta(\xi)$ for the variable of interest. Since in this application the constant $\alpha(\xi)$ can not be easily computed, we consider surrogates for (46) and (49) using a preconditioner $P_m(\xi)$. We consider

$$\Delta(\xi) = \|P_m(\xi)(A(\xi)u_r(\xi) - b(\xi))\|_{V_0} \|L(\xi)^* - A(\xi)^*Q_k(\xi)\|_{Z' \to V_0'} \quad (59)$$

for the primal-dual method, and

$$\Delta(\xi) = \|P_m(\xi)(A(\xi)t_{r,p}(\xi) - b(\xi))\|_{V_0} \sup_{0 \neq z' \in Z'} \inf_{y \in T_p} \frac{\|L(\xi)^*z' - A(\xi)^*y\|_{V_0'}}{\|z'\|_{Z'}} \quad (60)$$

for the saddle point method. Figure 6 shows statistics of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ for different numbers $m$ of interpolation points for the preconditioner. We see that the max-min ratio and the normalized standard deviation are decreasing with $m$: this indicates an improvement of the error estimate. Furthermore, the mean value of
of size k and compute $\tilde{s}(\xi) = (LR_V^{-1}A(\xi)^*\mathbf{W}_k^Q)Y_k(\xi)$.

* **Primal-dual:** solve the linear system of the Primal only method, solve the linear system

$$((\mathbf{W}_k^Q)^T A(\xi) R_V^{-1} A(\xi)^* \mathbf{W}_k^Q) Y_k(\xi) = ((\mathbf{W}_k^Q)^T b) - ((\mathbf{W}_k^Q)^T A(\xi) \mathbf{V}_r) U_r(\xi)$$

of size k, and compute $\tilde{s}(\xi) = (LV_r)U_r(\xi) + (LR_V^{-1}A(\xi)^*\mathbf{W}_k^Q)Y_k(\xi)$.

* **Saddle point:** solve the linear system of size $p+r$

$$\begin{pmatrix} \mathbf{T}_p^T(\xi) A(\xi) R_V^{-1} A(\xi)^* \mathbf{T}_p(\xi) & \mathbf{T}_p^T(\xi) A(\xi) \mathbf{V}_r \\ (\mathbf{T}_p^T(\xi) A(\xi) \mathbf{V}_r)^T & 0 \end{pmatrix} \begin{pmatrix} Y_{r,p}(\xi) \\ U_{r,p}(\xi) \end{pmatrix} = \begin{pmatrix} \mathbf{T}_p(\xi)^T b \\ 0 \end{pmatrix}$$

with $\mathbf{T}_p(\xi) = (\mathbf{W}_r(\xi), \mathbf{W}_k^Q)$, and compute

$$\tilde{s}(\xi) = (LV_r)U_{r,p}(\xi) + (LR_V^{-1}A(\xi)^*\mathbf{T}_p(\xi))Y_{r,p}(\xi).$$

The numerical results are given in Figure 5. Once again, the saddle point method leads to the lowest error on the variable of interest. Also, we see that a good preconditioner (for example with $m = 30$) improves the accuracy for the saddle point method, the primal only method and the primal-dual method. However, this improvement is not really significant for the considered application: the errors are barely divided by 2 compared to the non preconditioned Galerkin projection ($m = 0$). In fact, the preconditioner improves the quality of the test space, and the choice $W_r = V_r$ (yielding the standard Galerkin projection) is sufficiently accurate for this example and for the chosen norm on $V$.

We discuss now the quality of the error estimate $\Delta(\xi)$ for the variable of interest. Since in this application the constant $\alpha(\xi)$ can not be easily computed, we consider surrogates for (46) and (49) using a preconditioner $P_m(\xi)$. We consider

$$\Delta(\xi) = \|P_m(\xi)(A(\xi)u_r(\xi) - b(\xi))\|_{V_0} \|L(\xi)^* - A(\xi)^*Q_k(\xi)\|_{Z' \to V_0'} \quad (59)$$

for the primal-dual method, and

$$\Delta(\xi) = \|P_m(\xi)(A(\xi)t_{r,p}(\xi) - b(\xi))\|_{V_0} \sup_{0 \neq z' \in Z'} \inf_{y \in T_p} \frac{\|L(\xi)^*z' - A(\xi)^*y\|_{V_0'}}{\|z'\|_{Z'}} \quad (60)$$

for the saddle point method. Figure 6 shows statistics of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ for different numbers $m$ of interpolation points for the preconditioner. We see that the max-min ratio and the normalized standard deviation are decreasing with $m$: this indicates an improvement of the error estimate. Furthermore, the mean value of
<figure><img src="image_7.png" /><figcaption>Figure 5: Application 2: Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated from a samples set of size $10^4$.</figcaption></figure>

$\eta(\xi)$ seems to converge (with $m$) to 19.5 for the primal-dual method, and to 13.8 for the saddle point method. In fact, with a good preconditioner, $\|P_m(\xi)(A(\xi)u_r(\xi) - b(\xi))\|_{V_0}$ (or $\|P_m(\xi)(A(\xi)t_{r,p}(\xi) - b(\xi))\|_{V_0}$) is expected to be a good approximation of the primal error $\|u(\xi) - u_r(\xi)\|_{V_0}$ (or $\|u(\xi) - t_{r,p}(\xi)\|_{V_0}$), but this does not ensure that the effectivity index $\eta(\xi)$ will converge to 1.

### 4.2.3 Partial conclusions and remarks

In both numerical examples, the saddle point method provides the most accurate estimation for the variable of interest. Let us note that the saddle point problem requires the solution of a dense linear system of size $(r+k)$ for the symmetric and coercive case, and of size $(2r+k)$ for the general case. When using Gauss elimination method for the solution of those systems, the complexity is either in $C(r+k)^3$ or $C(2r+k)^3$ (with $C = 2/3$), which is larger than the complexity of the primal-dual method $C(r^3+k^3)$. However, in the case where the primal and dual approximation spaces have the same dimension $r=k$, the saddle point method is only 4 times (in the symmetric and coercive case) or 13.5 times (in the general case) more expensive.

For the present applications, we showed that the preconditioner slightly improves the quality of the estimation $\tilde{s}(\xi)$, and of the error estimate $\Delta(\xi)$. Since the construction of the preconditioner yields a significant increase in computational and memory costs (see
<figure><img src="image_7.png" /><figcaption>Figure 5: Application 2: Probability density function, $L^\infty$ norm and $L^2$ norm of the error $\|s(\xi) - \tilde{s}(\xi)\|_Z$ estimated from a samples set of size $10^4$.</figcaption></figure>

$\eta(\xi)$ seems to converge (with $m$) to 19.5 for the primal-dual method, and to 13.8 for the saddle point method. In fact, with a good preconditioner, $\|P_m(\xi)(A(\xi)u_r(\xi) - b(\xi))\|_{V_0}$ (or $\|P_m(\xi)(A(\xi)t_{r,p}(\xi) - b(\xi))\|_{V_0}$) is expected to be a good approximation of the primal error $\|u(\xi) - u_r(\xi)\|_{V_0}$ (or $\|u(\xi) - t_{r,p}(\xi)\|_{V_0}$), but this does not ensure that the effectivity index $\eta(\xi)$ will converge to 1.

### 4.2.3 Partial conclusions and remarks

In both numerical examples, the saddle point method provides the most accurate estimation for the variable of interest. Let us note that the saddle point problem requires the solution of a dense linear system of size $(r+k)$ for the symmetric and coercive case, and of size $(2r+k)$ for the general case. When using Gauss elimination method for the solution of those systems, the complexity is either in $C(r+k)^3$ or $C(2r+k)^3$ (with $C = 2/3$), which is larger than the complexity of the primal-dual method $C(r^3+k^3)$. However, in the case where the primal and dual approximation spaces have the same dimension $r=k$, the saddle point method is only 4 times (in the symmetric and coercive case) or 13.5 times (in the general case) more expensive.

For the present applications, we showed that the preconditioner slightly improves the quality of the estimation $\tilde{s}(\xi)$, and of the error estimate $\Delta(\xi)$. Since the construction of the preconditioner yields a significant increase in computational and memory costs (see
<figure><img src="image_8.png" /><figcaption>(a) PDF of $\eta(\xi)$ for the primal-dual methods and the saddle point methods. Three different preconditioners $P_m(\xi)$ are used: $m = 0$ (dotted lines), $m = 10$ (dashed lines) and $m = 30$ (continuous lines)</figcaption></figure>

<table><caption>(b) Statistics of the effectivity index $\eta(\xi)$ for the primal-dual method and the saddle point method.</caption><thead><tr><th></th><th>E(η(ξ))</th><th>max η(ξ)<br>min η(ξ)</th><th>Var(η(ξ))<sup>1/2</sup><br>E(η(ξ))</th></tr></thead><tbody><tr><td>Primal-dual</td><td></td><td></td><td></td></tr><tr><td>m = 0</td><td>5.545</td><td>3.52 × 10<sup>3</sup></td><td>2.246</td></tr><tr><td>m = 5</td><td>16.03</td><td>8.68 × 10<sup>2</sup></td><td>1.920</td></tr><tr><td>m = 10</td><td>18.69</td><td>1.01 × 10<sup>3</sup></td><td>1.925</td></tr><tr><td>m = 20</td><td>19.20</td><td>5.77 × 10<sup>2</sup></td><td>1.504</td></tr><tr><td>m = 30</td><td>19.59</td><td>3.95 × 10<sup>2</sup></td><td>1.615</td></tr><tr><td>Saddle point</td><td></td><td></td><td></td></tr><tr><td>m = 0</td><td>4.726</td><td>6.93 × 10<sup>3</sup></td><td>3.597</td></tr><tr><td>m = 5</td><td>12.61</td><td>1.80 × 10<sup>2</sup></td><td>1.429</td></tr><tr><td>m = 10</td><td>13.27</td><td>1.72 × 10<sup>2</sup></td><td>1.160</td></tr><tr><td>m = 20</td><td>13.97</td><td>1.89 × 10<sup>2</sup></td><td>1.090</td></tr><tr><td>m = 30</td><td>13.84</td><td>2.17 × 10<sup>2</sup></td><td>1.113</td></tr></tbody></table>

Figure 6: Application 2: PDF, mean, max-min ratio and normalized standard deviation of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$. Here, $\Delta(\xi)$ is defined by (59) for the primal-dual method and by (60) for the saddle point method.

[22]), the preconditioning is not mandatory here. Nevertheless, these results revealed the important role of the test space $W_r(\xi)$ to reduce the projection error. The preconditioner used for constructing $W_r(\xi)$ can be improved, for example with a better selection of the interpolation point for $P_m(\xi)$, see Equation (44). Note also that alternative methods can be also applied for constructing $W_r(\xi)$, such as the subspace interpolation method proposed in [1].

## 4.3 Greedy construction of the reduced spaces

We now consider the greedy construction of the reduced spaces by Algorithms 1 or 2. For the two considered applications, we show the convergence of the error estimate with respect to the complexity of the offline and of the online phase. For the sake of simplicity, we measure the complexity of the offline phase with the number of operator factorizations (this corresponds to the number of iterations $I$ of Algorithms 1 and 2). Of course exact estimation of the offline complexity should take into account many other steps (for example, the computation of $\Delta(\xi)$, of the preconditioner, etc), but the operator factorization is considered, for large scale applications, as the main source of computation cost. For the online complexity, we only consider the computation cost for the solution of one reduced system, see Section 4.2.3. Here we do not take into account the complexity for assembling the reduced systems although
<figure><img src="image_8.png" /><figcaption>(a) PDF of $\eta(\xi)$ for the primal-dual methods and the saddle point methods. Three different preconditioners $P_m(\xi)$ are used: $m = 0$ (dotted lines), $m = 10$ (dashed lines) and $m = 30$ (continuous lines)</figcaption></figure>

<table><thead><tr><th></th><th>E(η(ξ))</th><th>max η(ξ)<br>min η(ξ)</th><th>Var(η(ξ))<sup>1/2</sup><br>E(η(ξ))</th></tr></thead><tbody><tr><td>m = 0</td><td>5.545</td><td>3.52 × 10<sup>3</sup></td><td>2.246</td></tr><tr><td>m = 5</td><td>16.03</td><td>8.68 × 10<sup>2</sup></td><td>1.920</td></tr><tr><td>m = 10</td><td>18.69</td><td>1.01 × 10<sup>3</sup></td><td>1.925</td></tr><tr><td>m = 20</td><td>19.20</td><td>5.77 × 10<sup>2</sup></td><td>1.504</td></tr><tr><td>m = 30</td><td>19.59</td><td>3.95 × 10<sup>2</sup></td><td>1.615</td></tr></tbody></table>

<table><caption>(b) Statistics of the effectivity index $\eta(\xi)$ for the primal-dual method and the saddle point method.</caption><thead><tr><th></th><th>E(η(ξ))</th><th>max η(ξ)<br>min η(ξ)</th><th>Var(η(ξ))<sup>1/2</sup><br>E(η(ξ))</th></tr></thead><tbody><tr><td>m = 0</td><td>4.726</td><td>6.93 × 10<sup>3</sup></td><td>3.597</td></tr><tr><td>m = 5</td><td>12.61</td><td>1.80 × 10<sup>2</sup></td><td>1.429</td></tr><tr><td>m = 10</td><td>13.27</td><td>1.72 × 10<sup>2</sup></td><td>1.160</td></tr><tr><td>m = 20</td><td>13.97</td><td>1.89 × 10<sup>2</sup></td><td>1.090</td></tr><tr><td>m = 30</td><td>13.84</td><td>2.17 × 10<sup>2</sup></td><td>1.113</td></tr></tbody></table>

Figure 6: Application 2: PDF, mean, max-min ratio and normalized standard deviation of the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$. Here, $\Delta(\xi)$ is defined by (59) for the primal-dual method and by (60) for the saddle point method.

[22]), the preconditioning is not mandatory here. Nevertheless, these results revealed the important role of the test space $W_r(\xi)$ to reduce the projection error. The preconditioner used for constructing $W_r(\xi)$ can be improved, for example with a better selection of the interpolation point for $P_m(\xi)$, see Equation (44). Note also that alternative methods can be also applied for constructing $W_r(\xi)$, such as the subspace interpolation method proposed in [1].

## 4.3 Greedy construction of the reduced spaces

We now consider the greedy construction of the reduced spaces by Algorithms 1 or 2. For the two considered applications, we show the convergence of the error estimate with respect to the complexity of the offline and of the online phase. For the sake of simplicity, we measure the complexity of the offline phase with the number of operator factorizations (this corresponds to the number of iterations $I$ of Algorithms 1 and 2). Of course exact estimation of the offline complexity should take into account many other steps (for example, the computation of $\Delta(\xi)$, of the preconditioner, etc), but the operator factorization is considered, for large scale applications, as the main source of computation cost. For the online complexity, we only consider the computation cost for the solution of one reduced system, see Section 4.2.3. Here we do not take into account the complexity for assembling the reduced systems although
it may be a significant part of the complexity for “not so reduced” systems of equations.

### 4.3.1 Application 1

Figure 7 shows the convergence of $\sup_{\xi} \Delta(\xi)$ with respect to the offline and online complexities (as defined above). In Figure 7(a), we see that the saddle point method (dashed lines) always provides lower values for the error estimate compared to the primal-dual method (continuous lines). However, as already mentioned, the saddle point method requires the solution of larger reduced systems during the *online* phase. Therefore, the primal-dual method can sometimes provide lower error estimates (see the blue and red curves of Figure 7(b)) for the same *online* complexity.

The simultaneous construction of $V_r$ and $W_k^Q$ with full dual enrichment (52) (green curves) yields a very fast convergence of the error estimate during the *offline* phase, see Figure 7(a)). But the rapid increase of $\dim(W_k^Q)$ leads to high *online* complexity, so that this strategy becomes non competitive during the *online* phase, see Figure 7(b).

We compare now the alternate and the simultaneous construction of $V_r$ and $W_k^Q$ with partial dual enrichment (53) (red and blue curves in Figure 7). The initial idea of the alternate construction is to build reduced spaces of better quality. Indeed, since the evaluation points of the primal solution are different from the one of the dual solution, the reduced spaces are expected to contain “complementary information” for the approximation of the variable of interest. In practice, we observe in Figure 7(a) that the alternate construction is (two times) more expensive during the *offline* phase, but the resulting error estimate behaves very similarly to the simultaneous strategy, see Figure 7(b). We conclude that the alternate strategy is not relevant for this application.

Furthermore, let us note that after iteration 50 of the greedy algorithm, the rate of convergence of the dashed red curve of Figure 7(a) (i.e. the simultaneous construction with partial dual enrichment using the saddle point method) rapidly increases. A possible explanation is that the dimension of the dual approximation space is large enough to reproduce correctly the dual variable, which requires a dimension higher than $l = 44$. The same observation can be done for the alternative strategy (the dashed blue curve) after iteration 100 (which corresponds to $\dim(W_k^Q) \ge 50$). Also, we note that the primal-dual method does not present this behavior.

### 4.3.2 Application 2

For the application 2, we first test Algorithms 1 and 2 with the use of a preconditioner (defined in Section 4.2.2). The interpolation points for the preconditioner are the ones where the solutions (primal and dual) have been computed, see Algorithms 1 and 2. The
it may be a significant part of the complexity for “not so reduced” systems of equations.

### 4.3.1 Application 1

Figure 7 shows the convergence of $\sup_{\xi} \Delta(\xi)$ with respect to the offline and online complexities (as defined above). In Figure 7(a), we see that the saddle point method (dashed lines) always provides lower values for the error estimate compared to the primal-dual method (continuous lines). However, as already mentioned, the saddle point method requires the solution of larger reduced systems during the *online* phase. Therefore, the primal-dual method can sometimes provide lower error estimates (see the blue and red curves of Figure 7(b)) for the same *online* complexity.

The simultaneous construction of $V_r$ and $W_k^Q$ with full dual enrichment (52) (green curves) yields a very fast convergence of the error estimate during the *offline* phase, see Figure 7(a)). But the rapid increase of $\dim(W_k^Q)$ leads to high *online* complexity, so that this strategy becomes non-competitive during the *online* phase, see Figure 7(b).

We compare now the alternate and the simultaneous construction of $V_r$ and $W_k^Q$ with partial dual enrichment (53) (red and blue curves in Figure 7). The initial idea of the alternate construction is to build reduced spaces of better quality. Indeed, since the evaluation points of the primal solution are different from the one of the dual solution, the reduced spaces are expected to contain “complementary information” for the approximation of the variable of interest. In practice, we observe in Figure 7(a) that the alternate construction is (two times) more expensive during the *offline* phase, but the resulting error estimate behaves very similarly to the simultaneous strategy, see Figure 7(b). We conclude that the alternate strategy is not relevant for this application.

Furthermore, let us note that after iteration 50 of the greedy algorithm, the rate of convergence of the dashed red curve of Figure 7(a) (i.e. the simultaneous construction with partial dual enrichment using the saddle point method) rapidly increases. A possible explanation is that the dimension of the dual approximation space is large enough to reproduce correctly the dual variable, which requires a dimension higher than $l = 44$. The same observation can be done for the alternative strategy (the dashed blue curve) after iteration 100 (which corresponds to $\dim(W_k^Q) \ge 50$). Also, we note that the primal-dual method does not present this behavior.

### 4.3.2 Application 2

For the application 2, we first test Algorithms 1 and 2 with the use of a preconditioner (defined in Section 4.2.2). The interpolation points for the preconditioner are the ones where the solutions (primal and dual) have been computed, see Algorithms 1 and 2. The
<figure><img src="image_9.png" /><figcaption>Figure 7: Application 1: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity (Figure 7(a)) and the online complexity (Figure 7(b)). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as reference.</figcaption></figure>

preconditioner is used for the definition of the test space $W_r(\xi)$, see equation (44), and for the error estimate $\Delta(\xi)$, see equation (59) for the primal-dual method and (60) for the saddle point method. The numerical results are given in Figure 8. We can draw the same conclusions as for application 1.

* In the offline phase, the saddle point method provides lower errors (Figure 8(a)). However, the corresponding reduced systems are larger, and we see that the primal-dual method provides lower errors for the same online complexity, see Figure 8(b). For this test case, the benefits (in term of accuracy) of the saddle point method does not compensate the additional online computational costs.

* The full dual enrichment yields a fast convergence during the offline phase, but the rapid increase of $W_k^Q$ is disadvantageous regarding the online complexity. However, since the dimension of the variable of interest is “only” $l = 2$, the full dual enrichment is still an acceptable strategy (compared to the previous application).

* Here, the alternate strategy (blue curves) seems to yield slightly better reduced spaces compared to the simultaneous strategy, see Figure 8(b). But this leads to higher offline costs, see Figure 8(a).

We also run numerical tests without using the preconditioner. In that case, we replace $P_m(\xi)$ by $R_V^{-1}$. Figure 9 shows numerical results which are very similar to those of Figure
<figure><img src="image_9.png" /><figcaption>(a) Maximum value of the error estimate $\Delta(\xi)$ with respect to the number of operator factorizations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Maximum value of the error estimate $\Delta(\xi)$ with respect to the complexity of solving one reduced system.</figcaption></figure>

Figure 7: Application 1: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity (Figure 7(a)) and the online complexity (Figure 7(b)). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as reference.

preconditioner is used for the definition of the test space $W_r(\xi)$, see equation (44), and for the error estimate $\Delta(\xi)$, see equation (59) for the primal-dual method and (60) for the saddle point method. The numerical results are given in Figure 8. We can draw the same conclusions as for application 1.

* In the offline phase, the saddle point method provides lower errors (Figure 8(a)). However, the corresponding reduced systems are larger, and we see that the primal-dual method provides lower errors for the same online complexity, see Figure 8(b). For this test case, the benefits (in term of accuracy) of the saddle point method does not compensate the additional online computational costs.

* The full dual enrichment yields a fast convergence during the offline phase, but the rapid increase of $W_k^Q$ is disadvantageous regarding the online complexity. However, since the dimension of the variable of interest is “only” $l = 2$, the full dual enrichment is still an acceptable strategy (compared to the previous application).

* Here, the alternate strategy (blue curves) seems to yield slightly better reduced spaces compared to the simultaneous strategy, see Figure 8(b). But this leads to higher offline costs, see Figure 8(a).

We also run numerical tests without using the preconditioner. In that case, we replace $P_m(\xi)$ by $R_V^{-1}$. Figure 9 shows numerical results which are very similar to those of Figure
8. To illustrate the benefits of using the preconditioner, let us consider the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ associated to the error estimate for the variable of interest. Figure 10 shows the confidence interval $I(p)$ of probability $p$ for $\eta(\xi)$ defined as the smallest interval which satisfies

$$ \mathbb{P}(\xi \in \Xi_t : \eta(\xi) \in I(p)) \geq p, $$

where $\mathbb{P}(A) = \#A/\#\Xi_t$ for $A \subset \Xi_t$. When using the preconditioner, we see in Figure 10 that the effectivity index is improved during the greedy iterations in the sense that the confidence intervals are getting smaller and smaller. Also, we note that after the iteration 15, the effectivity index is always above 1: this indicates that the error estimate tends to be certified. Furthermore, after iteration 20 we do not observe any further improvement, so that it seems not useful to continue enriching the preconditioner.

Let us finally note that the use of the preconditioner yields significant computational costs. Indeed, we have to store operator factorizations (in our current implementation of the method), and the computation of the interpolation of the inverse operator requires additional problems to solve (see [22]). For the present application, even if the effectivity index of the error estimate is improved, the benefits of using the preconditioner remains questionable.

<figure><img src="image_10.png" /><figcaption>Figure 8: Application 2 with preconditioner: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity (Figure 8(a)) and the online complexity (Figure 8(b)). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as references.</figcaption></figure>
8. To illustrate the benefits of using the preconditioner, let us consider the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_Z$ associated to the error estimate for the variable of interest. Figure 10 shows the confidence interval $I(p)$ of probability $p$ for $\eta(\xi)$ defined as the smallest interval which satisfies

$$ \mathbb{P}(\xi \in \Xi_t : \eta(\xi) \in I(p)) \geq p, $$

where $\mathbb{P}(A) = \#A/\#\Xi_t$ for $A \subset \Xi_t$. When using the preconditioner, we see in Figure 10 that the effectivity index is improved during the greedy iterations in the sense that the confidence intervals are getting smaller and smaller. Also, we note that after the iteration 15, the effectivity index is always above 1: this indicates that the error estimate tends to be certified. Furthermore, after iteration 20 we do not observe any further improvement, so that it seems not useful to continue enriching the preconditioner.

Let us finally note that the use of the preconditioner yields significant computational costs. Indeed, we have to store operator factorizations (in our current implementation of the method), and the computation of the interpolation of the inverse operator requires additional problems to solve (see [22]). For the present application, even if the effectivity index of the error estimate is improved, the benefits of using the preconditioner remains questionable.

<figure><img src="image_11.png" /><figcaption>Figure 8: Application 2 with preconditioner: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity (Figure 8(a)) and the online complexity (Figure 8(b)). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as references.</figcaption></figure>
<figure><img src="image_11.png" /><figcaption>(a) Maximum value of the error estimate $\Delta(\xi)$ with respect to the number of operator factorization.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(b) Maximum value of the error estimate $\Delta(\xi)$ with respect to the complexity for solving one reduced system.</figcaption></figure>

Figure 9: Application 2 without preconditioner: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity 9(a) and the online complexity 9(b). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as references.

<figure><img src="image_13.png" /><figcaption>(a) Without preconditioner.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(b) With preconditioner.</figcaption></figure>

Figure 10: Application 2: evolution with respect to the greedy iterations of the confidence interval $I(p)$ for the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_z$ for the saddle point method.

# 5 Conclusion

We have proposed and analyzed projection based methods for the estimation of vector-valued variables of interest in the context of parameter-dependent equations. This includes a generalization of the classical primal-dual method to the case of vector-valued variables of
<figure><img src="image_12.png" /><figcaption>(a) Maximum value of the error estimate $\Delta(\xi)$ with respect to the number of operator factorization.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Maximum value of the error estimate $\Delta(\xi)$ with respect to the complexity for solving one reduced system.</figcaption></figure>

Figure 9: Application 2 without preconditioner: error estimate $\sup_{\xi} \Delta(\xi)$ with respect to the offline complexity 9(a) and the online complexity 9(b). The continuous lines correspond to the primal-dual method, and the dashed lines correspond to the saddle point method. The primal only curves serve as references.

<figure><img src="image_14.png" /><figcaption>(a) Without preconditioner.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>(b) With preconditioner.</figcaption></figure>

Figure 10: Application 2: evolution with respect to the greedy iterations of the confidence interval $I(p)$ for the effectivity index $\eta(\xi) = \Delta(\xi)/\|s(\xi) - \tilde{s}(\xi)\|_z$ for the saddle point method.

# 5 Conclusion

We have proposed and analyzed projection based methods for the estimation of vector-valued variables of interest in the context of parameter-dependent equations. This includes a generalization of the classical primal-dual method to the case of vector-valued variables of
interest, and also a Petrov-Galerkin method based on a saddle point problem. Numerical results showed that the saddle point method always improves the quality of the approximation compared to the primal-dual method using the same reduced spaces. We have also derived computable error estimates and greedy algorithms for the goal-oriented construction of the reduced spaces. The performances of these approaches have been compared on numerical examples, with an analysis of both the offline complexity (construction of the reduced spaces) and the online complexity (evaluation of the reduced order model and estimation of the variable of interest for one instance of the parameter). This complexity analysis revealed that the saddle point method is preferable to the primal-dual method regarding the offline costs. However, in the situation where the reduction of the online costs matter more than the reduction of offline costs, then the primal-dual method seems to be a better option (at least for the considered applications). For the considered applications, the use of preconditioners allows the construction of better reduced test spaces and also better error estimates. Even if the additional computational costs for building the preconditioner is significant, this has demonstrated the importance of having a suitable test space and good residual based error estimates.

The proposed error estimates, which involve the use of Cauchy-Schwarz inequalities, are clearly not optimal. Extending probabilistic error bounds proposed in [12] to the case of vector-valued variables could improve these error estimates.

## References

[1] D. AMSALLEM AND C. FARHAT, *Interpolation method for adapting reduced-order models and application to aeroelasticity*. AIAA Journal, 46(7):1803-1813, 2008.

[2] U. BAUR, C. BEATTIE, P. BENNER AND S. GUGERCIN, *Interpolatory projection methods for parameterized model reduction*. SIAM J. Sci. Comput., 33(5):2489-2518, 2011.

[3] P. BENNER, V. MEHRMANN AND D.C. SORENSEN, *Dimension reduction of large-scale systems*. Springer, 45, 2005.

[4] P. BENNER, S. GUGERCIN AND K. WILLCOX, *Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems*. SIAM Review, Vol. 57, 2015.

[5] G. BERKOOZ, P. HOLMES AND J. LUMLEY, *The proper orthogonal decomposition in the analysis of turbulent flows*. Ann. Rev. Fluid Mech., 25:539-575, 1993.
interest, and also a Petrov-Galerkin method based on a saddle point problem. Numerical results showed that the saddle point method always improves the quality of the approximation compared to the primal-dual method using the same reduced spaces. We have also derived computable error estimates and greedy algorithms for the goal-oriented construction of the reduced spaces. The performances of these approaches have been compared on numerical examples, with an analysis of both the offline complexity (construction of the reduced spaces) and the online complexity (evaluation of the reduced order model and estimation of the variable of interest for one instance of the parameter). This complexity analysis revealed that the saddle point method is preferable to the primal-dual method regarding the offline costs. However, in the situation where the reduction of the online costs matter more than the reduction of offline costs, then the primal-dual method seems to be a better option (at least for the considered applications). For the considered applications, the use of preconditioners allows the construction of better reduced test spaces and also better error estimates. Even if the additional computational costs for building the preconditioner is significant, this has demonstrated the importance of having a suitable test space and good residual based error estimates.

The proposed error estimates, which involve the use of Cauchy-Schwarz inequalities, are clearly not optimal. Extending probabilistic error bounds proposed in [12] to the case of vector-valued variables could improve these error estimates.

## References

[1] D. AMSALLEM AND C. FARHAT, *Interpolation method for adapting reduced-order models and application to aeroelasticity*. AIAA Journal, 46(7):1803-1813, 2008.

[2] U. BAUR, C. BEATTIE, P. BENNER AND S. GUGERCIN, *Interpolatory projection methods for parameterized model reduction*. SIAM J. Sci. Comput., 33(5):2489-2518, 2011.

[3] P. BENNER, V. MEHRMANN AND D.C. SORENSEN, *Dimension reduction of large-scale systems*. Springer, 45, 2005.

[4] P. BENNER, S. GUGERCIN AND K. WILLCOX, *Survey of Projection-Based Model Reduction Methods for Parametric Dynamical Systems*. SIAM Review, Vol. 57, 2015.

[5] G. BERKOOZ, P. HOLMES AND J. LUMLEY, *The proper orthogonal decomposition in the analysis of turbulent flows*. Ann. Rev. Fluid Mech., 25:539-575, 1993.
[6] Y. CHEN, J. S. HESTHAVEN, Y. MADAY AND J. RODRGUEZ, *Certified Reduced Basis Methods and Output Bounds for the Harmonic Maxwell's Equations*. SIAM J. Sci. Comput., 32(2):970–996, 2010.

[7] N. CUONG, K. VEROY AND A.T. PATERA, *Certified real-time solution of parametrized partial differential equations*. Handbook of Materials Modeling, 1523–1558, 2005.

[8] W. DAHMEN, C. PLESKEN AND G. WELPER, *Double greedy algorithms: Reduced basis methods for transport dominated problems*. ESAIM Math. Model. Numer. Anal., 48(03):623–663, May 2013.

[9] A. ERN AND J.-L. GUERMOND *Theory and practice of finite elements*. Springer Science & Business Media, Vol. 159, 2013.

[10] M. A. GREPL AND A. T. PATERA, *A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations*. ESAIM Math. Model. Numer. Anal., 39(1):157–181, 2005.

[11] B. HAASDONK, *Reduced basis methods for parametrized PDEs–A tutorial introduction for stationary and instationary problems*. Reduced Order Modelling. Luminy Book series, 2014.

[12] A. JANON, M. NODET AND C. PRIEUR, *Goal-oriented error estimation for the reduced basis method, with application to sensitivity analysis*. Journal of Scientific Computing, 1–15, 2015.

[13] D. B P HUYNH, G. ROZZA, S. SEN AND A. T. PATERA, *A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants*. Comptes Rendus Math., 345(8):473–478, October 2007.

[14] K. KAHLBACHER AND S. VOLKWEIN, *Galerkin proper orthogonal decomposition methods for parameter dependent elliptic systems*. Discuss. Math. Differential Incl., 27:95–117, 2007.

[15] Y. MADAY, A. T. PATERA, AND D. V. ROVAS, *A blackbox reduced-basis output bound method for noncoercive linear problems*. Studies Math. Appl., vol. 31, pp. 533-569, 2002.

[16] N. A. PIERCE AND M. B. GILES, *Adjoint Recovery of Superconvergent Functionals from PDE Approximations*. SIAM Rev., 42(2):247–264, 2000.
[6] Y. CHEN, J. S. HESTHAVEN, Y. MADAY AND J. RODRGUEZ, *Certified Reduced Basis Methods and Output Bounds for the Harmonic Maxwell's Equations*. SIAM J. Sci. Comput., 32(2):970–996, 2010.

[7] N. CUONG, K. VEROY AND A.T. PATERA, *Certified real-time solution of parametrized partial differential equations*. Handbook of Materials Modeling, 1523–1558, 2005.

[8] W. DAHMEN, C. PLESKEN AND G. WELPER, *Double greedy algorithms: Reduced basis methods for transport dominated problems*. ESAIM Math. Model. Numer. Anal., 48(03):623–663, May 2013.

[9] A. ERN AND J.-L. GUERMOND *Theory and practice of finite elements*. Springer Science & Business Media, Vol. 159, 2013.

[10] M. A. GREPL AND A. T. PATERA, *A posteriori error bounds for reduced-basis approximations of parametrized parabolic partial differential equations*. ESAIM Math. Model. Numer. Anal., 39(1):157–181, 2005.

[11] B. HAASDONK, *Reduced basis methods for parametrized PDEs–A tutorial introduction for stationary and instationary problems*. Reduced Order Modelling. Luminy Book series, 2014.

[12] A. JANON, M. NODET AND C. PRIEUR, *Goal-oriented error estimation for the reduced basis method, with application to sensitivity analysis*. Journal of Scientific Computing, 1–15, 2015.

[13] D. B P HUYNH, G. ROZZA, S. SEN AND A. T. PATERA, *A successive constraint linear optimization method for lower bounds of parametric coercivity and inf-sup stability constants*. Comptes Rendus Math., 345(8):473–478, October 2007.

[14] K. KAHLBACHER AND S. VOLKWEIN, *Galerkin proper orthogonal decomposition methods for parameter dependent elliptic systems*. Discuss. Math. Differential Incl., 27:95–117, 2007.

[15] Y. MADAY, A. T. PATERA, AND D. V. ROVAS, *A blackbox reduced-basis output bound method for noncoercive linear problems*. Studies Math. Appl., vol. 31, pp. 533-569, 2002.

[16] N. A. PIERCE AND M. B. GILES, *Adjoint Recovery of Superconvergent Functionals from PDE Approximations*. SIAM Rev., 42(2):247–264, 2000.
[17] C. PRUD'HOMME, D. V. ROVAS, K. VEROY, L. MACHIELS, Y. MADAY, A. T. PATERA AND G. TURINICI, *Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods*. J. Fluids Eng., 124(1):70, 2002.

[18] A. QUARTERONI, G. ROZZA AND A. MANZONI, *Certified reduced basis approximation for parametrized partial differential equations and applications*. J. Math. Ind., 1(1):3, 2011.

[19] G. ROZZA, D. B. P. HUYNH AND A. T. PATERA, *Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics*. Arch. Comput. Methods Eng., 15(3):229-275, May 2008.

[20] G. ROZZA AND K. VEROY, *On the stability of the reduced basis method for Stokes equations in parametrized domains*. Comput. Methods Appl. Mech. Eng., 196(7):1244-1260, 2007.

[21] J. XU, AND L. ZIKATANOV *Some observations on Babuska and Brezzi theories*. Numerische Mathematik, 94(1):195-202, 2003

[22] O. ZAHM AND A. NOUY, *Interpolation of inverse operators for preconditioning parameter-dependent equations*. SIAM J. Sci. Comput., 2016.
[17] C. PRUD'HOMME, D. V. ROVAS, K. VEROY, L. MACHIELS, Y. MADAY, A. T. PATERA AND G. TURINICI, *Reliable Real-Time Solution of Parametrized Partial Differential Equations: Reduced-Basis Output Bound Methods*. J. Fluids Eng., 124(1):70, 2002.

[18] A. QUARTERONI, G. ROZZA AND A. MANZONI, *Certified reduced basis approximation for parametrized partial differential equations and applications*. J. Math. Ind., 1(1):3, 2011.

[19] G. ROZZA, D. B. P. HUYNH AND A. T. PATERA, *Reduced basis approximation and a posteriori error estimation for affinely parametrized elliptic coercive partial differential equations: Application to transport and continuum mechanics*. Arch. Comput. Methods Eng., 15(3):229-275, May 2008.

[20] G. ROZZA AND K. VEROY, *On the stability of the reduced basis method for Stokes equations in parametrized domains*. Comput. Methods Appl. Mech. Eng., 196(7):1244-1260, 2007.

[21] J. XU, AND L. ZIKATANOV *Some observations on Babuska and Brezzi theories*. Numerische Mathematik, 94(1):195-202, 2003

[22] O. ZAHM AND A. NOUY, *Interpolation of inverse operators for preconditioning parameter-dependent equations*. SIAM J. Sci. Comput., 2016.