Xiaojun Lei
email: leixiaojun19@gscaep.ac.cn

Tongxiang Gu
email: txgu@iapcm.ac.cn

Stef Graillat
email: stef.graillat@lip6.fr

Hao Jiang
email: haojiang@nudt.edu.cn

Jin Qi
email: qi_jin@iapcm.ac.cn

A Fast Parallel High-Precision Summation Algorithm Based on AccSumK

Keywords: Parallel algorithms, Accurate summation, Higher precision, Floating-point arithmetic

In this paper, we present a new parallel accurate algorithm called PAccSumK for computing summation of floating-point numbers. It is based on AccSumK algorithm. In the experiment, for the summation problems with large condition numbers, our algorithm outperforms in term of accuracy and computing time the PSumK algorithm. The reason is that our algorithm is based on a more accurate algorithm called AccSumK algorithm compared to the SumL algorithm used in PSumK. The proposed parallel algorithms in this paper is designed to compute a result as if computed internally in K-fold the working precision. Numerical results are presented showing the performance and the accuracy of our new parallel algorithm for calculating summation.

Introduction and Motivation

In this paper, we present a new parallel algorithm to compute a sum of floating-point numbers with high accuracy. Since summation is one of the most basic tasks in numerical analysis, getting fast and accurate algorithms is an important problem. Moreover, computing accurate summation has numerous applications in numerical analysis and scientific computing.

The authors (Rump, Ogita and Oishi) presented an accurate sum algorithm called AccSumK [START_REF] Rump | Accurate floating-point summation I: faithful rounding[END_REF][START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF] using so-called "error-free transformations". Given a vector of floating-point numbers with exact sum s, their algorithm calculates a K-fold faithful rounding of s. The algorithm is fast for mildly conditioned sums with slowly increasing computing time proportional to the logarithm of the condition number. The algorithm is fast, not only in terms of the number of floating-point operations but also in terms of measured computing time on a serial computer. But it is not a parallel algorithm. The authors (Yamanaka, Ogita, Rump and Oishi) presented a parallel algorithm PSumK [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF]. The PSumK algorithm is based on the SumL [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF] and SumK [START_REF] Ogita | Accurate sum and dot product[END_REF] algorithms. However, the accuracy of SumL depends on the condition number of the sum whereas the accuracy of AccSumK is independent of the condition number. For ill-conditioned problems, the AccSumK algorithm is more accurate than the SumL algorithm.

Based on the above observations on AccSumK and PSumK, we develop a parallelizing method to used AccSumK into a new parallel algorithm called PAccSumK.

The rest of the paper is organized as follows: In Section 2, we briefly review the error-free transformations and the algorithm AccSumK for accurate summation proposed in [START_REF] Rump | Accurate floating-point summation I: faithful rounding[END_REF][START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF]. In Section 3, we develop parallel algorithm PAccSumK, which is based on AccSumK algorithm. We propose an analysis of the error bounds of the proposed algorithm. This makes it possible to confirm that the proposed parallel algorithm achieves a result as if computed in K-fold the working precision. In Section 4, we present results of numerical experiments showing the performance and the accuracy of the proposed algorithm PAccSumK. Finally, we conclude the paper and propose some future work.

Throughout the paper, floating-point addition, subtraction and multiplication are counted as one floating-point operation (flop). Moreover, we use MATLAB-style notation for describing algorithms.

Accurate summation algorithm AccSumK

In this section, we use the same notation as in [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF][START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF] and we recall some results and algorithms described in [START_REF] Ogita | Accurate sum and dot product[END_REF][START_REF] Rump | Accurate floating-point summation I: faithful rounding[END_REF][START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF][START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF]. The set of floating-point numbers is denoted by F, and U denotes the set of subnormal floating-point numbers together with zero and the two normalized floating-point numbers of smallest nonzero magnitude. Throughout this paper, we assume floatingpoint arithmetic adhering to IEEE standard 754 [START_REF]ANSI/IEEE, IEEE Standard for Binary Floating Point Arithmetic[END_REF]. Let p = (p i ) ∈ F n and let fl(•) be the result of floating-point operations, where all operations inside parentheses are executed by ordinary floating-point arithmetic in roundingto-nearest. We denote by u the machine epsilon, and the underflow unit by η, that is the smallest positive (denormalized) floating-point number. In IEEE standard 754 double precision u = 2 -53 and η = 2 -1074 . Let us denote s and S by

s := n i=1 p i , S := n i=1 |p i |.
The condition number of the summation of the vector p is defined by cond(

n i=1 p i ) := S |s| , s = 0.
A computed result Res of the summation of p = (p i ) is said to be computed in K-fold precision if it satisfies [START_REF] Ogita | Accurate sum and dot product[END_REF] |Res -

n i=1 p i | ≤ u|s| + (ϕu) K n i=1 |p i |, (1)
with a moderate constant ϕ. The second term in the right-hand side of (1) reflects computation K-fold precision, and the first one rounding back into working precision.

Because "unit in the last place" (ulp) depends on the floating-point format and needs extra care in the underflow range. It is useful to use the "unit in the first place" (ufp) defined in [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF] by

0 = r ∈ R ⇒ ufp(r) := 2 log 2 |r| , (2)
where ufp(0) := 0. In order to explain faithful rounding, the definition of floating-point predecessor and successor is required. pred(r) := max{f ∈ F : f < r} and succ(r

) := min{f ∈ F : f < r}. Definition 1. [3] A floating-point number f ∈ F is called a faithful rounding of a real number r ∈ R if pred(f ) < r < succ(f ). (3)
We denote this by f ∈ (r). For r ∈ F this implies f = r.

Definition 2. [3] A sequence f 1 , . . . , f k ∈ F is called a (K-fold) faithful rounding of s ∈ R if f i ∈ (s - i-1 ν=1 f ν ) f or 1 ≤ i ≤ k. (4)
Following [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF], we define γ n as

γ n := nu 1 -nu , n ∈ N.
When using γ n , we implicitly assume that nu < 1. Let p = (p 1 , ..., p n ) T ∈ F n . Then it holds that [START_REF] Higham | Accuracy and Stability of Numerical Algorithms[END_REF] s := fl(

n i=1 p i ) ⇒ |s - n i=1 p i | ≤ γ n-1 n i=1 |p i |. (5)
Note that ( 5) is valid for any order of addition in the summation. Recently, it was shown that γ n can be replaced by nu, and the restriction on n can be removed [START_REF] Lange | Sharp estimates for perturbation errors in summations[END_REF][START_REF] Rump | Error bounds for computer arithmetics[END_REF][START_REF] Jeannerod | On relative errors of floating-point operations: Optimal bounds and applications[END_REF][START_REF] Lange | Error estimates for the summation of real numbers with application to floating-point summation[END_REF]. The algorithm AccSumK is based on the error-free transformations of addition and/or multiplication of two floating-point numbers.

We first recall a fast algorithm for adding two floating-point numbers. It requires only 3 flops. When |a| ≥ |b|, Algorithm 2.1 is an error-free transformation that satisfies Choose a split to extract a vector into a sum of higher order parts and a vector of lower order parts. Algorithm 2.2 [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF]. Error-free vector transformation extracting high order part.

x + y = a + b with x = fl(a + b), |y| ≤ u • ufp(x).
function [τ, p ] = ExtractVector(σ, p) τ = 0 for i = 1:n We let the input vector p be p (0) and the output vector p be p (m) . The main part of the summation algorithm is Algorithm 2.4. It transforms the input vector p (0) into two floating-point numbers τ 1 , τ 2 representing the high order part of the sum and a vector p (m) of low order parts. The transformation is error-free, i.e., s := p

q i = fl((σ + p i ) -σ) p i = fl(p i -q i ) τ = fl(τ + q i ) end for Algorithm 2.3 [3]. Computation of 2 log 2 |p| for p = 0. function L = NextPowerTwo(p) q = u -1 p L = fl(|(q + p) -q|) if L == 0 L = |p| end if Theorem 1.
(0) i = τ 1 + τ 2 + p (m) i . Moreover, fl(τ 1 + τ 2 + p (m)
i ) is a faithfully rounded result of s (see [START_REF] Rump | Accurate floating-point summation I: faithful rounding[END_REF]). Algorithm 2.4 [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF]. Transformation of vector p plus with check for zero.

function [τ 1 , τ 2 , p, σ] = Transform(p, ) µ = max(|p i |) if µ = 0, τ 1 = , τ 2 = 0, return, end if M = NextPowerTwo(length(p 0 )+2) σ = 2 M NextPowerTwo(µ) t = repeat t = t ; σ = σ [τ, p] = ExtractVector(σ, p) t = fl(t + τ ) if t = 0, [τ 1 , τ 2 , p] = Transform(p, 0); return; end if σ = fl(2 M uσ) until σ ≤ 1 2 u -1 η or |t | ≥ fl(2 2M +1 uσ) [τ 1 , τ 2 ] = FastTwoSum(t, τ )
The following Algorithm 2.5 establishes an approximation of K-fold accuracy.

Algorithm 2.5 [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF]. Error-free vector transformation including faithful rounding.

function

[res, R, p ] = TransformK(p, ) [τ 1 , τ 2 , p , σ] = Transform(p, ) res = fl(τ 1 + (τ 2 +( n i=1 p i ))) R = fl(τ 2 -(res -τ 1 ))
Algorithm 2.6 computes a result of K-fold accuracy, stored in a nonoverlapping result vector Res of length K.

Algorithm 2.6 [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF]. Accurate summation with faithful rounding and result of K-fold accuracy.

function

Res = AccSumK(p, K) p (0) = p, R 0 = 0 for k = 1 : K [Res k , R k , p (k) ] = TransformK(p (k-1) , R k-1 )
if Res k ∈ U, Res k+1...K = 0, return; end if end for Figure 1 illustrates an outline of AccSumK for K = 3. Theorem 2. [START_REF] Rump | Accurate floating-point summation part II: sign, K-Fold faithful and rounding to nearest[END_REF] Let p be a vector of n floating-point numbers. Define M := log 2 (n + 2) , and assume 2 2M u ≤ 1. Let 1 ≤ K ∈ N be given, and let Res be the results vector of Algorithm 2.6 ( AccSumK) applied to p and K.

Then, Res 1 , . . . , Res K is a faithful rounding of s. Furthermore, Abbreviate R es

s = k ν=1 Res ν + R k + n i=1 p (k) i f or1 ≤ k ≤ K, (6)
:= K k=1 Res k . Then Res k ∈ U for some 1 < k < K implies R es = s, and Res k ∈ U implies |s -R es| < 2u K ufp(s) ≤ 2u K |s|, (7)
and

|s -R es| < 2u K ufp(Res 1 ) ≤ 2u K |Res 1 |. ( 8
) Remark 2. [3] AccSumK needs (4m + 5K + 3)n + O(m + K) flops if the "repeat-until" loop in the first extraction TransformK is executed m times.
Next, we briefly review some algorithms used in the SumK and SumL algorithms. The following Algorithm 2.7 is proposed by Knuth. It needs 6 flops. Algorithm 2.7 [START_REF] Knuth | The Art of Computer Programming: Seminumerical Algorithms[END_REF]. Error-free transformation of the sum of two floating-point numbers.

function

[x, y] = TwoSum(a, b) x = fl(a + b) q = fl(x -a) y = fl(((a -(x -q)) + (b -q)))
Using the TwoSum algorithm, it can be extended to the error-free transformation of vector summation. Algorithm 2.8 [START_REF] Ogita | Accurate sum and dot product[END_REF]. Error-free vector transformation with respect to the summation.

function q = VecSum(p)

q 1 = p 1 for i = 2 : n [q i , q i-1 ] = TwoSum(p i , q i-1 ) end for
In [START_REF] Ogita | Accurate sum and dot product[END_REF], the authors (Ogita, Rump and Oishi) developed an algorithm SumK, which calculates the summation using VecSum iteratively. Algorithm 2.9 [START_REF] Ogita | Accurate sum and dot product[END_REF]. Summation n i=1 p i for p ∈ F n as in K-fold precision by (K -1)-fold error-free vector transformation.

function res = SumK(p, K)

for i = 1 : K -1 p = VecSum(p) end for res = fl( n i=1 p i )
The error bounds of the result res by SumK are given as follows [START_REF] Ogita | Accurate sum and dot product[END_REF]:

Theorem 3. [1]
Let res be the result obtained by SumK, then

|res -s| ≤ (u + 3γ 2 n-1 )|s| + γ K 2(n-1) S. (9)
Moreover, if s = 0, then | res -s s | ≤ u + 3γ 2 n-1 + γ K 2(n-1) cond( p i ). (10)
In [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF], the authors developed a summation algorithm SumL whose result is represented by an array of L floating-point numbers.

Algorithm 2.10 [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF]. Summation n i=1 p i for p ∈ F n as in L-fold precision whose result is represented by a sum of L floating-point numbers.

function q = SumL(p, L)

for k = 1 : L -1 p(1 : n -k + 1) = VecSum(p(1 : n -k + 1)) q k = p n-k+1 end for q L = fl( n-L+1 i=1 p i )
For later use, we present an error bound and inequality for SumL. Theorem 4. [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF] Let q be the result obtained by Algorithm 2.10 ( SumL), abbreviate q := L k=1 q k , then

|q -s| ≤ γ L n-1 S, (11)
Moreover, if s = 0, then | q -s s | ≤ γ L n-1 cond( p i ). ( 12
)
Algorithm 2.10 also satisfies the following inequality:

L k=1 |q k | ≤ (1 + γ 2(n-1) )S. ( 13
)
Remark 3. Let L = K, divide the right term of ( 11) by the right term of ( 7), we can get

γ K n-1 S 2u K |s| ≈ (n -1) K 2 cond( p i ). (14)
From ( 14) we know that the error bounds of AccSumK is smaller than that of SumL.

Parallel algorithm

In this section, we will develop parallel algorithm for calculating sum internally in K-fold working precision. It is based on the AccSumK algorithm, which is named PAccSumK.

Suppose the number of CPUs to be P ≥ 2 on a distributed memory system. Then the CPUs can be numbered id from 1 to P , so 1 ≤ id ≤ P .

PAccSumK is the same algorithm as PSumK [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF], but the difference is that we use the AccSumK algorithm instead of SumL and SumK.

Next, we present here the concrete algorithm of PAccSumK.

Algorithm 3.1 (PAccSumK). A parallel algorithm based on the AccSumK algorithm.

function

Res = PAccSumK(p, K, P ) c = n/P , c 1 = n -c(P -1) % parallel private (index1, index2) if id == 1 index1 = 1 : c 1 else index1 = c 1 + c(id -2) + 1 : c 1 + c(id -1) end index2 = K(id-1)+1 : K • id q(index2) = AccSumK(p(index1),K) % end parallel Res = AccSumK(q, K) Remark 4.
When the AccSumK algorithm in the parallel part of Algorithm 3.1 is replaced with the SumL algorithm, and the AccSumK algorithm in the serial part is replaced with the SumK algorithm, it is the PSumK algorithm. In PSumK, computational cost for the parallelized part is (6K -5)c flops and the rest requires (6K -5)(P K -1) flops [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF]. When the AccSumK algorithm in Algorithm 3.1 is replaced with the SumL algorithm, we call it PSumL. In PSumL, computational cost for the parallelized part is (6K -5)c flops and the rest requires (6K -5)P K flops. This is a unified framework, the parallel part and the serial part can use different summation algorithms, we focus on the PAccSumK algorithm.

In PAccSumK, computational cost for the parallelized part is (4m + 5K + 3)c + O(m + K) flops and the rest requires (4m + 5K + 3)(P K) + O(m + K) flops. Therefore, the theoretical parallel efficiency in terms of the flops count becomes as follows:

r 1 = (4m + 5K + 3)c + O(m + K) (4m + 5K + 3)c + O(m + K) + (4m + 5K + 3)(P K) + O(m + K) ≈ 1 1 + P K c .
If c P K, then r 1 approaches to one. We present the following theorem for PAccSumK.

: | R es -s| < 2u K |s| + (1 + 2u K )2u K S. ( 15
)
Moreover, if s = 0, then | R es -s s | < 2u K + (1 + 2u K )2u K cond( p i ). ( 16
)
Figure 3 illustrates an outline of PAccSumK for K = 3 and P = 3. Proof . From the definition of vectors p, p (id) , q and q (id) , we collect the following relations:

n i=1 p i = P id=1 c id i=1 p (id) i = s, (17)
n i=1 |p i | = P id=1 c id i=1 |p (id) i | = S, (18)
P K j=1 q j = P id=1 K k=1 q (id) k , (19)
P K j=1 |q j | = P id=1 K k=1 |q (id) k |. ( 20
) Then | R es -s| = | R es - P K j=1 q j + P K j=1 q j -s| ≤ | R es - P K j=1 q j | + | P K j=1 q j -s|. ( 21
)
Here, using Res=AccSumK(q, K) and ( 7) yields

| R es - P K j=1 q j | < 2u K | P K i=1 q i | ≤ 2u K (| P K i=1 q i -s| + |s|), (22)
Inserting ( 22) into (21), we have

| R es -s| < 2u K |s| + (1 + 2u K )| P K i=1 q i -s|. (23)
It follows by ( 17) and ( 19) that

| P K j=1 q j -s| = | P id=1 K k=1 q (id) k - P id=1 c id i=1 p (id) i | ≤ P id=1 | K k=1 q (id) k - c id i=1 p (id) i |. ( 24
)
Here, recalling that q (id) =AccSumK(p (id) , K) and using [START_REF]ANSI/IEEE, IEEE Standard for Binary Floating Point Arithmetic[END_REF], we have

| P j=1 q j -s| < P id=1 (2u K | c id i=1 p (id) i |) ≤ 2u K S. (25)
Inserting ( 25) into (23), we finally have

| R es -s| < 2u K |s| + (1 + 2u K )2u K S,
which proves [START_REF] Lange | Error estimates for the summation of real numbers with application to floating-point summation[END_REF] and divided by |s| also proves [START_REF] Demmel | Fast and accurate floating point summation with application to computational geometry[END_REF]. Theorem 5 and (1) mean the result Res is obtained by PAccSumK as if computed internally in K-fold working precision. In fact, we can construct p such that the right-hand side equation of (25) holds.

Remark 5. The PAccSumK algorithm does not share the same accuracy as AccSumK. Indeed, the accuracy of PAccSumK depends on the condition number of the problem and so is not a K-fold faithful rounding. Because Res is the K-fold faithful rounding of q, which will produce an error, q (id) is the K-fold faithful rounding of p (id) , and each process will also produce an error, Res is not the K-fold faithful rounding of p.

In order to compare PAccSumK and PSumK, we quote the following theorem for PSumK. Theorem 6. [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF] Let res be the result obtained by PSumK. Define c := n P . Suppose 2 max{c, P K}u < 1. Then the following inequality holds:

|res -s| ≤ (u + 3γ 2 P K-1 )|s| + φ 1 S, (26)
where

φ 1 := (1 + u + 3γ 2 P K-1 )γ K c-1 + (1 + γ 2(c-1) )γ K 2(P K-1) . Moreover, if s = 0, then | res -s s | ≤ u + 3γ 2 P K-1 + φ 1 cond( p i ). ( 27
)
The second term at the right end ( 15) is smaller than the second term at the right end (26), therefore, when K is the same, PAccSumK can solve problems with larger condition numbers than PSumK. The first term at the right end [START_REF] Lange | Error estimates for the summation of real numbers with application to floating-point summation[END_REF] is also smaller than the first term at the right end (26). The error of the algorithm is dominated by the second term.

In order to compare PAccSumK and PSumL, we present the following theorem for PSumL.

| R es -s| ≤ {γ L c-1 + (1 + γ 2(c-1) )γ L n-1 }S. (28)
Moreover, if s = 0, then | R es -s s | ≤ {γ L c-1 + (1 + γ 2(c-1) )γ L n-1 }cond( p i ). ( 29
)
Proof . We have

| R es -s| = | R es - P L j=1 q j + P L j=1 q j -s| ≤ | R es - P L j=1 q j | + | P L j=1 q j -s|. (30)
Here, using Res=SumL(q, L) and ( 11) yields

| R es - P L j=1 q j | ≤ γ L n-1 P L i=1 |q i |, (31)
Inserting (31) into (30), we have

| R es -s| ≤ | P L i=1 q i -s| + γ L n-1 P L i=1 |q i |. (32)
It follows by ( 17) and ( 19) that

| P L j=1 q j -s| = | P id=1 L k=1 q (id) k - P id=1 c id i=1 p (id) i | ≤ P id=1 | L k=1 q (id) k - c id i=1 p (id) i |. (33)
Here, recalling that q (id) =SumL(p (id) , L) and using [START_REF] Rump | Inversion of extremely ill-conditioned matrices in floatingpoint[END_REF], we have

| P L j=1 q j -s| ≤ P id=1 (γ L c id -1 c id i=1 |p (id) i |) ≤ γ L c-1 S. (34)
Furthermore, applying [START_REF] Rump | Error bounds for computer arithmetics[END_REF] to the right-hand side of (20) and using (18) yield

P L j=1 |q j | = P id=1 ( L j=1 |q (id) j |) ≤ P id=1 ((1 + γ 2(c id -1) ) c id i=1 |p (id) i |) ≤ (1 + γ 2(c-1) )S.
(35) Inserting ( 35) and ( 34) into (32), we finally have

| R es -s| ≤ {γ L c-1 + (1 + γ 2(c-1) )γ L n-1 }S,
which proves (28) and divided by |s| also proves (29). Let L = K, divide the right term in (29) by the second term on the right in [START_REF] Demmel | Fast and accurate floating point summation with application to computational geometry[END_REF].

γ K c-1 + (1 + γ 2(c-1) )γ K n-1 (1 + 2u K )2u K ≈ (c -1) K + (n -1) K 2 . (36)
From (36) we know that the error bounds of PAccSumK is smaller than that of PSumL. The error bounds of PAccSumK is less than PSumL because the error bounds of AccSumK is less than SumL. We can also see that the error bounds of PAccSumK has nothing to do with the number of sums n, while the error bounds of PSumL is related to n. In the example 4 of our numerical results section, we take n = 2002, K = 7, P = 6, therefore, the second term in the error bounds of PAccSumK is approximately 6.42 × 10 22 times smaller than the error bounds of PSumL.

Numerical results

The following parallel experiment is performed on Sugon HPC cluster with 424 compute nodes, consisting of two 14-core Intel(R) Xeon(R) Gold 6162 2.60GHz processors each (28 cores per node). Nodes are connected by Intel Omni-Path high-speed computing network. The peak performance of the whole system is 980Tflops/s, and the memory is 96GB. The OS used by the cluster is Redhat7.2. The MPI library used for this experiment is Intel MPI. Accuracy is evaluated by relative error e = |res-s|/|s|, where res is an estimate of s. For PAccSumK, the terms Res k are non-overlapping, then res is equal to the sum of the array elements calculated by the ordinary summation algorithm. For PSumL, we sort Res k by magnitude and sum in decreasing order [START_REF] Rump | Inversion of extremely ill-conditioned matrices in floatingpoint[END_REF]. If the output of the algorithm is a floating-point number, then res is equal to the output of the algorithm, s is the exact value calculated by the MPFR library [START_REF] Fousse | MPFR: A multiple-precision binary floating-point library with correct rounding[END_REF] or known. When not stated separately, PAccSumK always takes K = 3.

Example 1

We use Algorithm 4.2 in Yamanaka et al. [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF] to generate arbitrarily illconditioned sum data. First generating ill-conditioned dot product data, the ill-conditioned sum data length is 2n generated from the dot product data of length n, and then the algorithm TwoProduct is used to convert the illconditioned dot product data of length n to the ill-conditioned sum data of length 2n through error-free transformation. Finally, randomly disturbing 2n summation data can generate ill-conditioned sum data with different condition numbers. Figure 3 is the relative error image of AccSumK and SumK. SumK is the serial version of PSumK. In this example, we use K = 2 to complete on the MATLAB plateform. From Figure 3 we can see that the AccSumK algorithm is more accurate than SumK algorithm.

Example 2

We use Algorithm 6.1 in Ogita et al. [START_REF] Ogita | Accurate sum and dot product[END_REF] to generate arbitrarily illconditioned sum data. Convert the ill-conditioned dot product data into ill-conditioned sum data in the same way as in Example 1.

The following are some experimental results and analysis of experimental results.

Figure 4 is the relative error image of different condition numbers when 4 processes and 10 processes are used, and the summation scale is 200 and 2000 respectively. PSum is the result of parallel recursive summation algorithm [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF]. The test results show that PAccSumK has higher accuracy than PSum. For using different processes, the relative error results of PAccSumK are almost Figure 5 is the CPU time image of the algorithm under different condition numbers when using 4 processes and 10 processes, and the summation scale is 200 and 2000 respectively. The PAccSumK algorithm uses more CPU time than the PSum algorithm.

Taking the PSum algorithm as a comparison benchmark, the calculation time of the algorithm is averaged, and then the ratio is calculated. The calculation time ratio of the PAccSumK is shown in table 1.

Example 3

We use the same data generation method as the ReproBLAS library [START_REF] Ahrens | Efficient reproducible floating point summation and BLAS[END_REF], i.e. p i = sin(2.0 * π * (i/n -0.5)), compare the relative error and CPU time of different processes, and compare the relative error and CPU time in different input data sizes.

Figure 6 (a) uses the process number of 2, 4, 8, 10, 16 and 20 to calculate the relative error when n = 1 × 10 i , i = 3, . . . , 7, the purpose is to compare the relative errors when using different number of processes. It can be seen from the image that the PAccSumK algorithm uses different processes to obtain almost the same summation result on the same set of summation data, so the relative error curve overlap. PSum got inaccurate results, PAccSumK obtains high-precision results. Figures 6 (b) and 6 (c) use process numbers of 4, 8, 10, and 20 when n = 2 × 10 i , i = 3, . . . , 7 and n = 3 × 10 i , i = 3, . . . , 7 calculate relative error, the purpose is to compare the relative error of the same process number on different n, the same process number has slightly different curves on the sum data of different n. PSum got inaccurate results, PAccSumK obtains high-precision results.

In terms of CPU time, from figure 7 we can observe PAccSumK algorithm takes more time than the PSum algorithm. Taking the PSum algorithm as a comparison benchmark, the calculation time of the algorithm is averaged, and then the ratio is calculated. The calculation time ratio of the PAccSumK is shown in table 2.

Table 3 shows the calculation speed of PAccSumK using different processes on n = 1 × 10 i , i = 3, . . . , 7 and n = 2 × 10 i , i = 3, . . . , 7 scale sum data, n = 1 × 10 i , i = 3, . . . , 7 sum data is based on the number of processes as 2, n = 2 × 10 i , i = 3, . . . , 7 sum data is based on the number of processes 4 as a benchmark to calculate the acceleration when using other processes.

It can be seen from table 3 that the CPU time used by the PAccSumK algorithm is approximately inversely proportional to the number of processes used.

Example 4

We use the same ill-conditioned sum data as in Example 1.

The speed-up ratio R 1 and the parallel efficiency R are defined by R 1 := T 1 T P and R := R 1 P . where T P means the elapsed time in case of using P CPUs (1 ≤ P ≤ 8). From R 1 , we can see the scalability of the algorithms. If the floating-point operations are perfectly parallelized, then it holds that T P := C + T 1 T P , where C denotes the elapsed time for the parallelization overhead.

The following are some experimental results and analysis of experimental results. Next, we evaluate how tight the error bound for PAccSumK in Theorem 5 is in practice. To do this, we set n = 2002 and vary the condition number cond from 1 to 10 140 in Algorithm 4.2 [START_REF] Yamanaka | A parallel algorithm for accurate dot product[END_REF], the exact result of sum is equal to cond -1 . The error bounds [START_REF] Demmel | Fast and accurate floating point summation with application to computational geometry[END_REF] and true relative errors of the results obtained by PAccSumK for K = 2, 4, 8 in case of P = 6 are displayed in figure 11. In figure 11, the lines labeled 'exp.' for each K denote the experimental error corresponding to K, the lines labeled 'est.' for each K denote the error bounds (16) corresponding to K.

Finally, we evaluate the elapsed time and parallel efficiency of PAccSumK for n = 10 6 and n = 10 7 . For both cases we set cond = 1 • 10 100 in Algorithm 4.2 and vary K = 7, 8, 9 in PAccSumK. The general rule to choose K is based on [START_REF] Demmel | Fast and accurate floating point summation with application to computational geometry[END_REF] if the condition number and the expected accuracy are known. Here, we give an example to choose K. Suppose we calculate in double precision, we need the relative error between the approximate solution and the exact solution to meet the machine accuracy, in this example, cond = 1 • 10 100 , we can calculate K ≈ 7.29. Therefore, in this example, we choose K = 7. The elapsed times of PAccSumK for both cases are displayed in table 4. Moreover, figure 12 and 13 illustrate the scalability and the parallel efficiency of PAccSumK compared with AccSumK for n = 10 6 and n = 10 7 , respectively.

From table 4, figure 12 and 13, we can observe the following facts: • As the number of processes increases, the speed-up ratio increases approximately linearly, indicating that the PAccSumK algorithm has better scalability; • In both scales, K = 7 has the highest speed-up ratio and parallel efficiency.

Under the premise that the results can be obtained with high accuracy, it is best to take K = 7 in this example;

Conclusions and Future Work

In this paper, we proposed the parallel algorithm for accurate sum. We proved that our proposed algorithm achieved a result as if computed with Kfold the working precision. By the numerical experiments, we confirmed that the proposed parallel algorithm works effectively on our computer environment. For the summation problem with large condition number, PAccSumK performs better than PSumK and PSumL in terms of accuracy and time. In practical applications, large condition number sequences often appear in the case of massive cancellation, for example, in computational geometry [START_REF] Demmel | Fast and accurate floating point summation with application to computational geometry[END_REF]. In practice, the proposed algorithm frequently gives more accurate results than the theoretical error bounds.

In the future, we plan to replace the ordinary summation and dot product of linear equations solver with high-precision summation and dot product in order to obtain a high-precision linear equations solver. We also plan to use AccSumK algorithm to obtain a parallel high-precision reproducible summation algorithm.

Algorithm 2 .

 2 1 [6]. Compensated summation of two floating-point numbers. function [x, y] = FastTwoSum(a, b) x = fl(a + b) q = fl(x -a) y = fl(b -q)

[ 3 ]

 3 Let L be the result of Algorithm 2.3 ( NextPowerTwo) applied to a nonzero floating-point number p. If no overflow occurs, then L = 2 log 2 |p| . Remark 1. [3] L = 2 log 2 |p| is satisfied in the presence of underflow, and use Algorithm 2.3 (NextPowerTwo) so that only basic floating-point operations are necessary.

Figure 1 :

 1 Figure 1: Outline of AccSumK for K = 3

Theorem 5 .

 5 Let Res be the result obtained by Algorithm 3.1 ( PAccSumK), abbreviate R es := K k=1 Res k . Define c := n P , M := log 2 (c + 2) , and assume 2 2M u≤ 1. Let 1 ≤ K ∈ N be given. Then the following inequality holds

Figure 2 :

 2 Figure 2: Outline of PAccSumK for K = 3 and P = 3

Theorem 7 .

 7 Let Res be the result obtained by PSumL, abbreviate R es := L k=1 Res k . Define c := n P . Suppose 2 max{c, P K}u < 1. Then the following inequality holds:

Figure 3 :

 3 Figure 3: Relative error image of serial algorithm

Figure 4 :

 4 Figure 4: relative error images (left n = 200) (right n = 2000)

Figure 5 :

 5 Figure 5: CPU time (left n = 200) (right n = 2000)

Figure 6 :

 6 Figure 6: Relative error image (a) n = 1 × 10 i , i = 3, . . . , 7, (b) n = 2 × 10 i , i = 3, . . . , 7 and (c) n = 3 × 10 i , i = 3, . . . , 7.

Figure 7 :

 7 Figure 7: CPU time (left n = 1 × 10 i , i = 3, . . . , 7) (right n = 2 × 10 i , i = 3, . . . , 7)

Figure 8 (

 8 Figure 8 (a) is the relative error image under different condition numbers when using 4 processes and K = 7, figure 8 (b) is the relative error image under different condition numbers when using 6 processes and K = 7. This experiment uses the mpi intelmpi -2017.4.239. The test results show that PAccSumK is more accurate than the other three parallel summation algorithms.

Figure 8 : 6 Figure 9

 869 Figure 8: Relative error image (a) P = 4, (b) P = 6

Figure 9 : 1 Figure 10

 9110 Figure 9: Relative error image (a) mpi mpich2 -gnu -1.5, (b) mpi mpich -gnu -3.3.1 Figure 10 is the CPU time image when n = 10 6 , P = 4, and K = 7. PAccSumK algorithm is 1.7574 times slower than PSum algorithm. PAccSumK

Figure 10 :

 10 Figure 10: CPU time image

Figure 12 :

 12 Figure 12: Scalability (left) and parallel efficiency (right) of PAccSumK for n = 10 6

Figure 13 :

 13 Figure 13: Scalability (left) and parallel efficiency (right) of PAccSumK for n = 10 7

Table 1 :

 1

		4 processes calculation time ratio
			n = 200		n = 2000
	Algorithm PSum PAccSumK	PSum PAccSumK
	time ratio	1	2.222	1	2.298

Table 2 :

 2

			4 processes' calculation time ratio
		n = 1 × 10 i , i = 3, . . . , 7	n = 2 × 10 i , i = 3, . . . , 7
	Algorithm PSum	PAccSumK	PSum	PAccSumK
	time ratio	1	21.834	1	21.352

Table 3 :

 3 calculation time ratioPAccSumK n = 1 × 10 i , i = 3, • • • , 7 n = 2 × 10 i , i = 3, • • • , 7

	p = 2	1	-
	p = 4	0.510	1
	p = 8	0.256	0.518
	p = 10	0.207	0.403
	p = 16	0.146	-
	p = 20	0.130	0.254

Acknowledgment

To begin with, we would like to thank the reviewers for their thorough reading of the article as well as their valuable comments and suggestions. The second and fourth authors were supported in part by Science Challenge Project (TZ2016002), NSF of China (61472462, 11671049, 11601033) and the foundation of key laboratory of computational physics. The third author was supported by the NuSCAP (ANR-20-CE48-0014) project of the French National Agency for Research (ANR).


Mohamed-Raed Blel

Virginie Ehrlacher

Tony Lelièvre

Influence of sampling on the convergence rates of greedy algorithms for parameter-dependent random variables

 where the authors proposed a variance reduction technique for the computation of parameter-dependent expectations using a reduced basis paradigm. We study the effect of Monte-Carlo sampling on the theoretical properties of greedy algorithms. In particular, using concentration inequalities for the empirical measure in Wasserstein distance proved in [14], we provide sufficient conditions on the number of samples used for the computation of empirical variances at each iteration of the greedy procedure to guarantee that the resulting method algorithm is a weak greedy algorithm with high probability. These theoretical results are not fully practical and we therefore propose a heuristic procedure to choose the number of Monte-Carlo samples at each iteration, inspired from this theoretical study, which provides satisfactory results on several numerical test cases.

In other words, the subspace V n is a quasi-optimal subspace of dimension n for the approximation of random variables f µ (Z) for µ ∈ P in an L 2 norm sense.

However, in practice, variances cannot be computed exactly and have to be approximated by empirical means involving a finite number of samples of the random vector Z, which may be different from one iteration of the greedy algorithm to another. The main result of this article is to give theoretical lower bounds on the number of samples which have to be taken at each iteration of the greedy algorithm in order to guarantee that the resulting Monte-Carlo greedy algorithm enjoys quasi-optimality properties close to those of an ideal greedy algorithm with high probability.

, which plays a fundamental role for model-order reduction.

 on the mathematical analysis of greedy algorithms in Hilbert spaces. In Section 3, we present the Monte-Carlo greedy algorithm, which is the main focus of this article, our main theoretical result and its proof. This theoretical result does not yield a fully practical algorithm. To alleviate this difficulty, we propose in Section 4 a heuristic algorithm, inspired from the theoretical result, which provides satisfactory results on several test cases.

Introduction

The aim of this article is to provide a mathematical study of the algorithm proposed in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF] where the authors proposed a variance reduction technique for the computation of parameter-dependent expectations using a reduced basis paradigm.

More precisely, the problematic we are considering here is the following: let us denote by P ⊂ R m a set of parameter values. In several applications, it is of significant interest to be able to rapidly compute the expectation of a random variable of the form f µ (Z) for a large numbers of values of the parameter µ ∈ P, where Z is a random vector and where for all µ ∈ P, f µ is a real-valued function. In practice, such expectations may not be computable analytically and are approximated using empirical means involving a large number of random samples of the random vector Z. Variance reduction methods are commonly used in such contexts in order to reduce the computational cost of approximating these expectations by means of standard Monte-Carlo algorithms. Among these, control variates, which are chosen as approximations of the random variable f µ (Z) the expectation of which can be easily computed, can yield to interesting gains in terms of computational cost, provided that the variance of the difference between f µ (Z) and its approximation is small. The construction of efficient control variates for a given application is thus fundamental for the variance reduction technique to yield significant computational gains.

In [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF], the authors proposed a general algorithm in order to construct a control variate for f µ (Z) using a reduced basis paradim. More precisely, the approximation of f µ (Z) is constructed as a linear combination of f µ1 (Z), • • • , f µn (Z) for some small integer n ∈ N * and well-chosen values µ 1 , • • • , µ n ∈ P of the parameters. The choice of n and of the values of the parameters stems from an iterative procedure, called a greedy algorithm, which consists at iteration n ∈ N to compute

µ n+1 ∈ argmax µ∈P inf Zn∈Vn Var [f µ (Z) -Z n ] ,
where V n := Span {f µ1 (Z), • • • , f µn (Z)}. In the ideal (unpractical) case where variances can be exactly computed, the procedure boils down to a standard greedy algorithm in a Hilbert space [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]. It is now well-known [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF] that such a greedy procedure provides a quasi-optimal set of parameters µ 1 , • • • , µ n in the sense that the error sup

µ∈P inf Zn∈Vn Var [f µ (Z) -Z n ] = inf Zn∈Vn Var f µn+1 (Z) -Z n
is comparable to the so-called Kolmogorov n-width of the set {f µ (Z), µ ∈ P}, defined by

sup µ∈P inf W n vectorial subspace dimW n = n inf Zn∈Wn Var [f µ (Z) -Z n ] .
2 Motivation: greedy algorithms for reduced bases and variance reduction

Motivation: reduced basis control variate

The aim of this section is to present the motivation of our work, which aims at constructing control variates for reducing the variance of a Monte-Carlo estimator of the mean of parameter-dependent functions of random vectors.

Let us begin by introducing some notation. Let d ∈ N * , (Ω, F, P) be a probability space and Z a R d -valued random vector with associated probability measure ν. For all q ∈ N * , we denote by

L q ν (R d ) := f : R d → R, R d |f (x)| q dν(x) < +∞ .
Let C(R d ) denote the set of continuous real-valued functions defined on R d . Let p ∈ N * , P ⊂ R p be a set of parameter values, and for all µ ∈ P, let f µ be an element of C(R d ) ∩ L 2 ν (R d ). For all f, g ∈ C(R d ), any M ∈ N * and any collection Z := (Z k ) 1≤k≤M of random vectors of R d , we define the empirical averages:

E Z (f ) := 1 M M k=1 f (Z k ) , Cov Z (f, g) := E Z (f g) -E Z (f )E Z (g), Var Z (f ) := Cov Z (f, f ).
The aim of our work is to propose and analyse from a mathematical point of view a numerical method in order to efficiently construct control variates to reduce the variance of a Monte-Carlo estimator of E [f µ (Z)] for all µ ∈ P using a Reduced Basis paradigm [START_REF] Jan S Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF][START_REF] Barrault | An 'empirical interpolation'method: application to efficient reduced-basis discretization of partial differential equations[END_REF][START_REF] Quarteroni | Reduced basis methods for partial differential equations: an introduction[END_REF][START_REF] Devore | The theoretical foundation of reduced basis methods. Model Reduction and approximation: Theory and Algorithms[END_REF], which was originally proposed in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF]. be two independent collections of iid random vectors distributed according to the law of Z and independent of Z.

Let us assume that we have selected N values of parameters (µ 1 , µ 2 , ..., µ N ) ∈ P N for some N ∈ N * and assume that the empirical means (E Z ref (f µi )) 1≤i≤N have been computed in an offline phase.

In an online phase, for all µ ∈ P, we can build an approximation of E [f µ (Z)], using a control variate which reads as f µ (Z) for some function f µ : R d → R:

E [f µ (Z)] ≈ E Z ref (f µ ) + E Z small f µ -f µ . (1)
Remark 2.1. Let us point out that the statistical error between

E Z ref (f µ ) and E f µ (Z) is close to Var f µ (Z) M ref ,
whereas the error between

E Z small f µ -f µ and E (f µ -f µ )(Z) is of the order of Var f µ -f µ (Z) M small .
The aim of the Monte-Carlo greedy algorithm studied in this article is to give an approximation of

E [f µ (Z)]
with an error close to

Var[fµ(Z)] M ref
within a much smaller computational time than the one required by the computation of

E Z ref (f µ ).
In the method studied here, the control variate function f µ is constructed as follows:

f µ = N i=1 λ µ i f µi
where λ µ := (λ µ i ) 1≤i≤N ∈ R N is a solution of the linear system

Aλ µ = b µ (2)
where

A := (A ij ) 1≤i,j≤N ∈ R N ×N and b µ := (b µ i ) 1≤i≤N ∈ R N are defined as follows: for all 1 ≤ i, j ≤ N , A ij = Cov Z small (f µi , f µj ) and b µ i = Cov Z small (f µ , f µi ). (3)
Equivalently, the vector λ µ is a solution of the minimization problem

λ µ ∈ argmin λ:=(λi) 1≤i≤N ∈R N Var Z small f µ - N i=1 λ i f µi .
Let us point out that λ µ is a random vector which can be written as a deterministic function of Z small . In other words, λ µ is measurable with respect to Z small . Remarking that

E Z ref (f µ ) = N i=1 λ µ i E Z ref (f µi ), the computation of the approximation (1) of E [f µ (Z)
] thus requires the following steps:

• offline phase: Compute (E Z ref (f µi )) 1≤i≤N (N empirical means with M ref samples), (E Z small (f µi )) 1≤i≤N
(N empirical means with M small samples) and the matrix A (N 2 empirical covariances with M small samples).

• online phase: For all µ ∈ P, compute b µ (N empirical covariances with M small samples) and solve the linear system (2) to obtain λ µ . Then, compute the approximation (1

) of E [f µ (Z)] as E [f µ (Z)] ≈ N i=1 λ µ i E Z ref (f µi ) + E Z small (f µ ) - N i=1 λ µ i E Z small (f µi ) , (4)
which requires O(N ) elementary operations and the computation of one empirical mean with M small samples.

Naturally, the approximation of E [f µ (Z)] given by (1) can be interesting from a computational point of view in terms of variance reduction only if Var f µ (Z) -f µ (Z) is much smaller than Var [f µ (Z)]. The following question thus naturally arises: how can the set of parameters (µ 1 , µ 2 , .., µ N ) ∈ P N be chosen in the offline phase in order to ensure that Var f µ (Z) -f µ (Z) is as small as possible for any value of µ ∈ P?

Greedy algorithms stand as the state-of-the-art technique to construct such sets of snapshot parameters, enjoy very nice mathematical properties and are the backbone of the method proposed in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF] which we wish to analyze here. We present this family of algorithms and related existing theoretical convergence results in the next section.

Greedy algorithms for reduced basis

Let us recall here the results of [START_REF] Buffa | A priori convergence of the greedy algorithm for the parametrized reduced basis method[END_REF][START_REF] Binev | Convergence rates for greedy algorithms in reduced basis methods[END_REF][START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF] on the convergence rates of greedy algorithms for reduced bases, adapted to our context. Let us define

L 2 ν,0 (R d ) := g ∈ L 2 ν (R d ), R d g dν = 0 . It holds that L 2 ν,0 (R d
) is a Hilbert space, equipped with the scalar product •, • defined by

∀g 1 , g 2 ∈ L 2 ν,0 (R d ), g 1 , g 2 = R d g 1 g 2 dν = Cov [g 1 (Z), g 2 (Z)] .
The associated norm is denoted by • and is given by

∀g ∈ L 2 ν,0 (R d ), g = R d |g| 2 dν 1/2 = Var [g(Z)].
For all µ ∈ P, let us define

g µ := f µ -E [f µ (Z)] (5)
and let us denote by M := {g µ , µ ∈ P} [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF] so that M ⊂ L 2 ν,0 (R d ). Let us assume that M is a compact subset of L 2 ν,0 (R d ). For all n ∈ N * , we introduce the Kolmogorov n-width of the set M in L 2 ν,0 (R d ), defined by

d n (M) := inf V n ⊂ L 2 ν,0 (R d ) subspace, dim V n = n sup µ∈P inf gn∈Vn Var [g µ (Z) -g n (Z)] = inf V n ⊂ L 2 ν,0 (R d ) subspace, dim V n = n sup µ∈P inf gn∈Vn g µ -g n .
Let 0 < γ < 1 and consider the following weak greedy algorithm with parameter γ.

Weak-Greedy Algorithm:

Initialization: Find µ 1 ∈ P such that g µ1 2 ≥ γ 2 max µ∈P g µ 2 . ( 7
) Set V 1 := Span{g µ1 } and set n = 2. Iteration n ≥ 2: Find µ n ∈ P such that inf (λi) 1≤i≤n-1 ∈R n-1 g µn - n-1 i=1 λ i g µi 2 ≥ γ 2 max µ∈P inf (λi) 1≤i≤n-1 ∈R n-1 g µ - n-1 i=1 λ i g µi 2 , (8)
Set V n := V n-1 + Span{g µn } = Span{g µ1 , • • • , g µn }.
For all n ∈ N * , the error associated with the n-dimensional subspace V n given by the weak greedy algorithm is defined by

σ n (M) := max µ∈P inf (λi) 1≤i≤n ∈R n g µ - n i=1 λ i g µi .
The following result is then a direct corollary of the results proved in [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF]Corollary 3.3].

Theorem 2.2. For all n ∈ N * , σ n (M) ≤ √ 2γ -1 min 0≤m<n (d m (M))
n-m n . In particular, for all n ∈ N * ,

σ 2n (M) ≤ √ 2γ -1 d n (M).
This result indicates that the weak greedy algorithm provides a practical way to construct a quasi-optimal sequence (V n ) n∈N * of finite dimensional subspaces of L 2 ν,0 (R d ). Of course, the weak greedy algorithm introduced above cannot be implemented in practice since it requires at the n th iteration of the algorithm the computation of the exact variances of g µ (Z)-

n-1 i=1 λ i g µi (Z) for µ, µ 1 , • • • , µ n-1 ∈ P and λ 1 , • • • , λ n-1 ∈ R,
which is out of reach in our context. In practice, these quantities have to be approximated by Monte-Carlo estimators involving a finite number of samples of the random vector Z. The resulting greedy algorithm with Monte Carlo sampling is presented in Section 3. The mathematical analysis of this algorithm is the main purpose of the present article.

For the sake of simplicity, in the rest of the article, we assume that for all n ∈ N * , d n (M) > 0.

3 Greedy algorithm with Monte-Carlo sampling

Presentation of the algorithm

Let us begin by presenting the greedy algorithm with Monte Carlo sampling.

Let (M n ) n∈N * be a sequence of integers, which represents the number of samples used at iteration n. For all n ∈ N * , let Z n := (Z n k ) 1≤k≤Mn be a collection of random vectors such that (Z n k ) n≥1, 1≤k≤Mn are independent and identically distributed according to the law of Z, and independent of Z. Let ν,0 (R d ), we define

g 1 , g 2 Z 1:∞ := Cov g 1 (Z), g 2 (Z) Z 1:∞ and g 1 Z 1:∞ := Var g 1 (Z) Z 1:∞ .
Let us make here an important remark. Since Z 1:∞ is a collection of random vectors which are all independent of Z, it holds that, for all f, g ∈ L 2 ν,0 (R d ), almost surely,

f, g Z 1:∞ = Cov f (Z), g(Z) Z 1:∞ = Cov [f (Z), g(Z)] = f, g , g 2 Z 1:∞ = Var g(Z) Z 1:∞ = Var[g(Z)] = g 2 .
Hence, almost surely, •, • Z 1:∞ defines a scalar product on L 2 ν,0 (R d ), which is a Hilbert space when equipped with this scalar product, and • Z 1:∞ is the associated norm.

The greedy algorithm with Monte-Carlo sampling reads as follows:

MC-Greedy Algorithm: Initialization: Find µ 1 ∈ P such that, almost surely,

µ 1 ∈ argmax µ∈P Var Z 1 (g µ ) and g µ 1 = 0. (9) Set V 1 := Span{g µ 1 } and set n = 2.
Iteration n ≥ 2: Find µ n ∈ P such that, almost surely,

µ n ∈ argmax µ∈P inf (λi) 1≤i≤n-1 ∈R n-1 Var Z n g µ - n-1 i=1 λ i g µ i and g µ n / ∈ V n-1 . ( 10
) Set V n := V n-1 + Span{g µ n } = Span{g µ 1 , • • • , g µ n }.
Naturally, for all n ∈ N * , the parameter µ n and thus the finite-dimensional space V n are Z 1:n -measurable.

Let us first prove an auxiliary lemma.

Lemma 3.1. Almost surely, all the iterations of the MC-Greedy Algorithm are well-defined, in the sense that, for all n ∈ N * , there always exists at least one element µ n ∈ P such that (9) (when n = 1) or (10) (when n ≥ 2) is satisfied.

Proof of Lemma 3.1. Let us first consider the initialization step corresponding to n = 1. Two situations may a priori occur : either

max µ∈P Var Z 1 (g µ ) > 0 or max µ∈P Var Z 1 (g µ ) = 0. In the first case, choosing µ 1 ∈ argmax µ∈P Var Z 1 (g µ ) is sufficient to guarantee that g µ 1 = 0. Indeed, since M ⊂ C(R d ) (remember that f µ
is continuous for all µ ∈ P, and hence so is g µ ), the fact that Var Z 1 g µ 1 > 0 necessarily implies that

Var g µ 1 (Z) Z 1:∞ > 0 almost surely. Since Z 1:∞ is independent of Z and µ 1 is a Z 1:∞ measurable random
variable, this implies that almost surely g µ 1 = 0.

In the second case, it then holds that Var Z 1 (g µ ) = 0 for all µ ∈ P. Then, the assumption d 1 (M) > 0 implies that, almost surely, there exists at least one element µ 1 ∈ P such that g µ 1 = 0. In addition,

µ 1 ∈ argmax µ∈P Var Z 1 (g µ ).
Using similar arguments and the fact that d n (M) > 0 for all n ∈ N * , it is easy to see that, almost surely, all the iterations of the MC-Greedy algorithm are well-defined, in particular for n ≥ 2.

Remark 3.2. We stress on the fact that the practical implementation of the MC-greedy algorithm does not require the knowledge of the value of

E[f µ (Z)], even if g µ = f µ -E[f µ (Z)] for all µ ∈ P. Indeed, it holds that for all g ∈ C(R d ), all n ∈ N * and all C ∈ R, Var Z n (g) = Var Z n (g + C). Thus, for all µ ∈ P, n ∈ N * and λ := (λ i ) 1≤i≤n-1 ∈ R n-1 , Var Z 1 (g µ ) = Var Z 1 (f µ )
and

Var Z n g µ - n-1 i=1 λ i g µ i = Var Z n f µ - n-1 i=1 λ i f µ i .
Thus, the MC-greedy algorithm naturally makes sense with a view to the construction of a reduced basis control variate for variance reduction as explained in Section 2.1.

Remark 3.3. In practice, a discrete subset P trial ⊂ P has to be introduced. The optimization problems [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF] and [START_REF] Devore | Greedy algorithms for reduced bases in banach spaces[END_REF] have to be replaced respectively by

µ 1 ∈ argmax µ∈P trial Var Z 1 (g µ ) and g µ 1 = 0, and
µ n ∈ argmax µ∈P trial inf (λi) 1≤i≤n-1 ∈R n-1 Var Z n g µ - n-1 i=1 λ i g µ i and g µ n / ∈ V n-1 .
The influence of the choice of the set P trial on the mathematical properties of the MC-greedy algorithm is an important question which we do not address in our analysis for the sake of simplicity. For related discussion, we refer the reader to the work [START_REF] Cohen | Reduced basis greedy selection using random training sets[END_REF], where the authors study the mathematical properties of a greedy algorithm where the set P trial depends on the iteration n of the greedy algorithm and is randomly chosen according to appropriate probability distributions defined on the set of parameters P.

For all n ∈ N * , we also define

σ n-1 (M) := max µ∈P inf (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g µ i (Z) Z 1:∞ , (11)
σ n-1 (M) := inf (λi) 1≤i≤n-1 ∈R n-1 Var g µ n (Z) - n-1 i=1 λ i g µ i (Z) Z 1:∞ (12)
Let us point out here that σ n-1 (M) is a random variable which is measurable with respect to Z

1:(n-1)
whereas µ n and σ n-1 (M) are measurable with respect to Z 1:n .

Main theoretical result

The aim of this section is to study the effect of Monte-Carlo sampling on the convergence of such a greedy algorithm. We consider here the probability space (Ω, A(Z 1:∞ ), P) the probabilty space where A(Z 1:∞ ) denotes the set of events that are measurable with respect to Z 1:∞ . We prove, under appropriate assumptions on the probability density ν and on the set of functions M = {g µ , µ ∈ P}, that for all 0 < γ < 1, there exist explicit conditions on the sequence (M n ) n∈N * so that, with high probability, the MC-greedy algorithm is actually a weak greedy algorithm with parameter γ. More precisely, under this set of assumptions, we prove that, with high probability, it holds that for all n ∈ N * ,

σ n-1 (M) ≥ γ σ n-1 (M).
Let us now present the set of assumptions we make on ν and on the set M = {g µ , µ ∈ P} for our main result to hold.

From now on, we make the following assumption on the probability distribution ν.

Assumption (A):

The probability law ν is such that there exist α > 1 and β > 0 such that R e β|x| α dν(x) < +∞.

Let us denote by L the set of Lipschitz functions of R d and for all f ∈ L, let us denote by f L its Lipschitz constant. In the sequel, we denote by φ : R * + → R * + the function defined by

∀κ ∈ R * + , φ(κ) :=    κ 2 1 κ≤1 + κ α 1 κ>1 if d = 1, (κ/ log(2 + 1/κ) 2 )1 κ≤1 + κ α 1 κ>1 if d = 2, κ d 1 κ≤1 + κ α 1 κ>1 if d ≥ 3. (13)
A key ingredient in our analysis is the use of concentration inequalities in the Wasserstein-1 distance between a probability distribution and its empirical measure proved in [START_REF] Bolley | Quantitative concentration inequalities for empirical measures on non-compact spaces[END_REF][START_REF] Fournier | On the rate of convergence in wasserstein distance of the empirical measure[END_REF]. Let us recall here a direct corollary of Theorem 2 of [START_REF] Fournier | On the rate of convergence in wasserstein distance of the empirical measure[END_REF], which is the backbone of our analysis.

Corollary 3.4. Let us assume that ν satisfies assumption (A). Then, there exist positive constants c, C depending only on ν, d, α and β, such that, for all M ∈ N * , all Z := (Z k ) 1≤k≤M iid random vectors distributed according to ν and all κ > 0, it holds that

P T 1 Z ≥ κ ≤ Ce -cM φ(κ) ,
where

T 1 Z := sup f ∈L; f L ≤1 |E[f (Z)] -E Z (f )| .
Remark 3.5. We would like to mention here that other concentration inequalities are stated in Theorem 2 of [START_REF] Fournier | On the rate of convergence in wasserstein distance of the empirical measure[END_REF] under different sets of assumptions than (A) on the probability law ν. In particular, weaker concentration inequalities may be obtained when ν only has some finite polynomial moments. Our analysis can then be easily adapted to these different settings but we restrict ourselves here to a framework where ν satisfies Assumption (A) for the sake of clarity.

We finally make the following set of assumptions on M defined in [START_REF] Boyaval | A variance reduction method for parametrized stochastic differential equations using the reduced basis paradigm[END_REF].

Assumption (B):

The set M satisfies the four conditions:

(B1) M is a compact subset of L 2 ν,0 (R d ) and let K 2 := sup µ∈P g µ < ∞;

(B2) M ⊂ L and K L := sup µ∈P g µ L < +∞;

(B3) M ⊂ L ∞ (R d ) and K ∞ := sup µ∈P g µ L ∞ < +∞; (B4) for all n ∈ N * , d n (M) > 0.
Before presenting our main result, we need to introduce some additional notation. Using Lemma 3.1, we can almost surely define the sequence (g n ) n∈N * as the orthonormal family of L 2 ν,0 (R d ) obtained by a Gram-Schmidt orthonormalization procedure (for • Z 1:∞ ) from the family (g µ n ) n∈N * . More precisely, we define

g 1 := g µ 1 Var g µ 1 (Z) Z 1:∞ . Moreover, for all n ≥ 2, et λ n := λ n i 1≤i≤n-1
∈ R n-1 be a solution to the minimization problem

λ n ∈ argmin λ:=(λi) 1≤i≤n-1 ∈R n-1 Var g µ n (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ .
Then it holds that

g n := g µ n - n-1 i=1 λ n i g i Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ .
As a consequence, it always holds that

V n = Span g µ 1 , • • • , g µ n = Span {g 1 , • • • , g n }. Moreover, g n is Z 1:n -measurable.
We are now in position to state our main result, the proof of which is postponed to Section 3.3.

Theorem 3.6. Let 0 < δ < 1 and (δ n ) n∈N * ⊂ (0, 1) N * be a sequence of numbers satisfying

Π n∈N * (1 -δ n ) ≥ 1 -δ.
Let us assume that M satisfies assumption (B) and that ν satisfies assumption (A). Let C, c > 0 be the constants defined in Corollary 3.4. For all n ∈ N * , let

K n ∞ := max (K ∞ , g 1 L ∞ , • • • , g n L ∞ ) and K n L := max (K L , g 1 L , • • • , g n L ) . (14)
Let us assume that there exists 0 < γ < 1 such that for all n ∈ N * , M n ∈ N * is a Z 1:(n-1) measurable random variable which satisfies almost surely the following condition:

∀n ≥ 1, M n ≥ -ln δ n C 1 cφ (κ n-1 ) , ( 15
)
where κ n-1 is a deterministic function of Z 1:(n-1) , defined by

κ 0 := 1 -γ 2 σ 0 (M) 2 8K ∞ K L ; (16)
and

∀n ≥ 2, κ n-1 := min 1 2(n-1) , (1-γ 2 ) σn-1(M) 2 n(9K 2 2 +4) 6K n-1 ∞ K n-1 L . ( 17
)
Then, for all n ∈ N * , it holds that

P σ n-1 (M) ≥ γ σ n-1 (M) Z 1:(n-1) ≥ 1 -δ n . ( 18
)
As a consequence, denoting by G n the event σ n-1 (M) ≥ γ σ n-1 (M) for all n ∈ N * , it holds that

P n∈N * G n ≥ 1 -δ. (19)
Thus, it then holds that the MC-greedy algorithm is a weak greedy algorithm with parameter γ and norm • Z 1:∞ with probability at least 1 -δ.

We state here a direct corollary of Theorem 3.6, the proof of which is given below.

Corollary 3.7. Under the assumptions of Theorem 3.6, with probability 1 -δ, it holds that for all n ∈ N * ,

σ n (M) ≤ √ 2γ -1 min 1≤m<n (d m (M)) n-m n . ( 20
)
In particular, with probability 1 -δ, it holds that

∀n ∈ N * , σ 2n (M) ≤ √ 2γ -1 d n (M). (21)
Proof. With probability 1 -δ, the MC-greedy algorithm is a weak greedy algorithm with parameter γ and norm • Z 1:∞ . Thus, since for all n ∈ N * , µ n is a Z 1:∞ measurable random variable, if such an event is realized, using Theorem 2.2, it holds that for all n ∈ N *

σ n (M) ≤ √ 2γ -1 min 1≤m<n d Z 1:∞ m (M) n-m n
, where for all n ∈ N * ,

d Z 1:∞ n (M) := inf V n ⊂ L 2 ν,0 (R d ) subspace, dim V n = n sup µ∈P inf gn∈Vn Var g µ (Z) -g n (Z) Z 1:∞ = inf V n ⊂ L 2 ν,0 (R d ) subspace, dim V n = n sup µ∈P inf gn∈Vn Var [g µ (Z) -g n (Z)] = d n (M).
Hence, we obtain (20), and (21) as a consequence.

Some remarks are in order here.

Remark 3.8. Note that, since the random variables

K n-1 ∞ , K n-1
L and σ n-1 (M) are measurable with respect to Z 1:(n-1) , κ n-1 is also measurable with respect to Z 1:(n-1) .

Remark 3.9. A natural question is then the following: can Theorem 3.6 be used (at least in principle) to design a constructive strategy to choose a number of samples M n , so that the MC-greedy algorithm can be guaranteed to be a weak greedy algorithm with parameter γ? This can indeed be done in principle using the following remark: for all n ∈ N * , the quantity σ n-1 (M) defined by ( 11) cannot be computed in practice since variances cannot be computed exactly for any parameter µ ∈ P. However, almost surely, it holds that σ n-1 (M) defined by ( 12) satisfies σ n-1 (M) ≤ σ n-1 (M). Let us recall that σ n-1 (M) depends on Z 1:n , whereas σ n-1 (M) only depends on Z 1:(n-1) . Since φ is an increasing function, this implies that, if the sequence (M n ) n∈N * satisfies condition [START_REF] Jan S Hesthaven | Certified reduced basis methods for parametrized partial differential equations[END_REF] where σ 0 (M) is replaced by σ 0 (M) in ( 16) and σ n-1 (M) is replaced by σ n-1 (M) in [START_REF] Smetana | Randomized residual-based error estimators for the proper generalized decomposition approximation of parametrized problems[END_REF], the assumptions of Theorem 3.6 are satisfied. Besides, it is reasonable to expect in this case that σ n-1 (M) should provide a reasonable approximation of σ n-1 (M) since, from Theorem 3.6, σ n-1 (M) ≥ γ σ n-1 (M) with high probability.

Unfortunately, we will see that such an approach is not viable in practice, because it leads to much too large values of M n for small values of n for the MC-greedy algorithm to be interesting with a view to the variance reduction method explained in Section 2.1. That is why in Section 4, we will present numerical results with heuristic ways to choose values of (M n ) n∈N * which are not theoretically guaranteeed, but which nevertheless yield satisfactory numerical results in several test cases.

Proof of Theorem 3.6

The aim of this section is to prove Theorem 3.6. For all n ∈ N * , we denote by G n the event σ n (M) ≥ γ σ n (M).

Let us begin by proving some intermediate results which will be used later. We first need the following auxiliary lemma. Lemma 3.10. Let n ∈ N * . Then, almost surely,

sup f Z 1:∞ -measurable random function such that f L ≤ 1 almost surely E f (Z)|Z 1:∞ -E Z n (f ) = sup f ∈L; f L ≤1 |E [f (Z)] -E Z n (f )| .
Proof. On the one hand, it is obvious to check that

sup f ∈L; f L ≤1 |E [f (Z)] -E Z n (f )| ≤ sup f Z 1:∞ -measurable random function such that f L ≤ 1 almost surely E f (Z)|Z 1:∞ -E Z n (f ) .
On the other hand, for any Z 1:∞ -measurable random function f such that f L ≤ 1 almost surely, it holds that, almost surely, since Z

1:∞ is independent of Z, E[f (Z)|Z 1:∞ ] = E Z [f (Z)]
where the index Z in E Z indicates that the expectation is only taken with respect to Z, and thus

E f (Z)|Z 1:∞ -E Z n (f ) ≤ sup f ∈L; f L ≤1 |E [f (Z)] -E Z n (f )| .
Hence the result.

We start by considering the case of the initialization of the MC-greedy algorithm.

Lemma 3.11. Let 0 < γ < 1. Then, it holds that almost surely,

P Var g µ 1 (Z) Z 1:∞ ≥ γ 2 max µ∈P Var g µ (Z) Z 1:∞ ≥ 1 -δ 1 . ( 22
)
As a consequence, P [G 1 ] ≥ 1 -δ 1 and ( 18) holds for n = 1.

Proof. Let µ 1 ∈ P such that

σ 0 (M) 2 = max µ∈P Var g µ (Z) Z 1:∞ = Var g µ1 (Z)|Z 1:∞ .
Inequality (22) holds provided that

P Var g µ1 (Z)|Z 1:∞ -Var g µ 1 (Z) Z 1:∞ > σ 0 (M) 2 ≤ δ 1 , with := 1 -γ 2 . Almost surely, since µ 1 ∈ argmax µ∈P Var Z 1 (g µ ), it holds that Var g µ1 (Z)|Z 1:∞ -Var g µ 1 (Z) Z 1:∞ = Var g µ1 (Z)|Z 1:∞ -Var Z 1 (g µ1 ) + Var Z 1 (g µ1 ) -Var Z 1 g µ 1 + Var Z 1 g µ 1 -Var g µ 1 (Z) Z 1:∞ ≤ Var g µ1 (Z)|Z 1:∞ -Var Z 1 (g µ1 ) + Var Z 1 g µ 1 -Var g µ 1 (Z) Z 1:∞ = E |g µ1 | 2 (Z)|Z 1:∞ -E Z 1 |g µ1 | 2 + E Z 1 (g µ1 ) 2 -E g µ1 (Z) Z 1:∞ 2 -E |g µ 1 | 2 (Z)|Z 1:∞ + E Z 1 |g µ 1 | 2 -E Z 1 g µ 1 2 + E g µ 1 (Z)|Z 1:∞ 2 ≤ E |g µ1 | 2 (Z)|Z 1:∞ -E Z 1 |g µ1 | 2 + 2K ∞ E Z 1 (g µ1 ) -E [g µ1 (Z)] + E |g µ 1 | 2 (Z)|Z 1:∞ -E Z 1 |g µ 1 | 2 + 2K ∞ E Z 1 g µ 1 -E g µ 1 (Z)|Z 1:∞ , ≤ 2K ∞ K L × E |g µ1 | 2 2K ∞ K L (Z)|Z 1:∞ -E Z 1 |g µ1 | 2 2K ∞ K L + E |g µ 1 | 2 2K ∞ K L (Z)|Z 1:∞ -E Z 1 |g µ 1 | 2 2K ∞ K L + E Z 1 g µ 1 K L -E g µ 1 K L (Z)|Z 1:∞ + E Z 1 g µ1 K L -E g µ1 K L (Z)|Z 1:∞ . It holds that for all µ ∈ P, |g µ | 2 L ≤ 2K ∞ K L .
Indeed, for all x, y ∈ R d , we have

||g µ | 2 (x) -|g µ | 2 (y)| = |(g µ (x) + g µ (y))(g µ (x) -g µ (y))| ≤ 2K ∞ K L |x -y|.
Thus, almost surely, it holds that

Var g µ1 (Z)|Z 1:∞ -Var g µ 1 (Z) Z 1:∞ ≤ 8K ∞ K L sup f ∈L; f L ≤1 E f (Z)|Z 1:∞ -E Z 1 (f ) .
Then, using Lemma 3.10, we obtain that, almost surely,

Var g µ1 (Z)|Z 1:∞ -Var g µ 1 (Z) Z 1:∞ ≤ 8K ∞ K L sup f ∈L; f L ≤1 E [f (Z)] -E Z 1 (f ) .
Thus, using Theorem 3.4, the assumption on M 1 and the definition of κ 0 , we obtain that

P sup f ∈L; f L ≤1 E [f (Z)] -E Z 1 (f ) ≥ κ 0 ≤ Ce -cφ(κ0) ≤ δ 1 .
Hence the desired result.

We now turn to the case of the n th iteration of the algorithm, with n ≥ 2, that we analyze in the next two lemmas. Lemma 3.12. Let n ≥ 2. Let us denote by

M n-1 := M ∪ {g 1 , • • • , g n-1 }. ( 23
)
Then, for all > 0, it holds that, almost surely,

P sup g,h∈M n-1 Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) ≥ Z 1:(n-1) ≤ Ce -cMnφ 6K n-1 ∞ K n-1 L , where K n-1 L and K n-1
∞ are defined by [START_REF] Fournier | On the rate of convergence in wasserstein distance of the empirical measure[END_REF].

Proof. For all g, h ∈ M n-1 , it holds that, almost surely,

Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) ≤ E g(Z)h(Z) Z 1:∞ -E Z n (gh) + K n-1 ∞ E g(Z) Z 1:∞ -E Z n (g) + E h(Z) Z 1:∞ -E Z n (h) ≤ 2K n-1 ∞ K n-1 L E gh 2K n-1 ∞ K n-1 L (Z) Z 1:∞ -E Z n gh 2K n-1 ∞ K n-1 L + E g 2K n-1 L (Z) Z 1:∞ -E Z n g 2K n-1 L + E h 2K n-1 L (Z) Z 1:∞ -E Z n h 2K n-1 L . For all g, h ∈ M n-1 , it holds that gh 2K n-1 ∞ K n-1 L L ≤ 1 and g 2K n-1 L L ≤ 1.
This implies that, almost surely,

sup g,h∈M n-1 Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) ≤ 6K n-1 ∞ K n-1 L sup f ∈L, f L ≤1 E f (Z) Z 1:∞ -E Z n (f ) .
Using Lemma 3.10, this yields that, almost surely,

sup g,h∈M n-1 Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) ≤ 6K n-1 ∞ K n-1 L sup f ∈L, f L ≤1 |E [f (Z)] -E Z n (f )| .
We finally obtain, using Corollary 3.4, that

P sup g,h∈M n-1 Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) > Z 1:(n-1) ≤ P sup f ∈L, f L ≤1 |E [f (Z)] -E Z n (f )| > 6K n-1 ∞ K n-1 L Z 1:(n-1) ≤ P sup f ∈L, f L ≤1 |E [f (Z)] -E Z n (f )| > 6K n-1 ∞ K n-1 L Z 1:(n-1) ≤ Ce -cMnφ 6K n-1 ∞ K n-1 L .
Hence the result.

Lemma 3.13. Let 0 < γ < 1 and n ≥ 2. Then, it holds that almost surely

P G n |Z 1:(n-1) ≥ 1 -δ n .
Proof. Since (g 1 , • • • , g n ) forms a basis of V n-1 , for all µ ∈ P, there exists one unique minimizer to

min (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ . Let λ n := λ n i 1≤i≤n-1
∈ R n-1 be the unique minimizer of

λ n := argmin λ:=(λi) 1≤i≤n-1 Var g µ n (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ . ( 24
)
As a consequence, it holds that

σ n-1 (M) = Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞
where σ n-1 (M) is defined by [START_REF] Devore | The theoretical foundation of reduced basis methods. Model Reduction and approximation: Theory and Algorithms[END_REF]. Let µ n ∈ P such that

µ n ∈ argmax µ∈P min (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ , so that σ n-1 (M) = min (λi) 1≤i≤n-1 ∈R n-1 Var g µn (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ ,
where σ n-1 (M) is defined in [START_REF] Devore | Nonlinear approximation[END_REF].

Let λ n := λ n i 1≤i≤n-1 ∈ R n-1 the unique minimizer of

λ n := argmin λ:=(λi) 1≤i≤n-1 Var g µn (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ , (25)
so that

σ n-1 (M) = Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ .
The event G n holds if and only if

Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≥ γ 2 Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ . ( 26
)
Let us begin by pointing out that, since

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≥ Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ , if the inequality Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≤ (1 -γ 2 )Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ (27)
holds, then ( 26) is necessarily statisfied. The rest of the proof consists in estimating the probability that ( 27) is realized. To this aim, as a first step, we are going to prove an upper bound on

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ as a function of η := sup g,h∈Mn-1 Cov g(Z), h(Z)| Z 1:∞ -Cov Z n (g, h) , (28)
which is the quantity estimated in Lemma 3.12. More precisely, let us now prove that

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≤ n 2 + K 2 + √ n -1η 1 -(n -1)η 2 + K 2 2 η. ( 29
)
It holds that for all 1 ≤ i ≤ n -1, from ( 24) and ( 25),

λ n i = Cov g µn (Z), g i (Z)|Z 1:∞ and λ n i = Cov g µ n (Z), g i (Z)|Z 1:∞ ,
and it then holds that, almost surely,

max λ n 2 , λ n 2 ≤ max g µn , g µ n ≤ K 2 , (30)
where

• 2 denotes the Euclidean norm of R n-1 . Let now λ n,n := λ n,n i 1≤i≤n-1 ∈ R n-1 be a minimizer of λ n,n := argmin λ:=(λi) 1≤i≤n-1 ∈R n-1 Var Z n g µn - n-1 i=1 λ i g i , and λ n,n := λ n,n i 1≤i≤n-1 ∈ R n-1 be a minimizer of λ n,n := argmin λ:=(λi) 1≤i≤n-1 ∈R n-1 Var Z n g µ n - n-1 i=1 λ i g i
It then holds that for all 1 ≤ i ≤ n -1, λ n,n and λ n,n are solution to the linear systems

A n λ n,n = b n and A n λ n,n = b n ,
where

A n := A n ij 1≤i,j≤n-1 ∈ R (n-1)×(n-1) , b n := b n i 1≤i≤n-1 , b n := b n i 1≤i≤n-1 ∈ R n-1 are defined as follows: for all 1 ≤ i, j ≤ n -1, A n ij = Cov Z n g i , g j , b n i = Cov Z n (g µn , g i ) and b n i = Cov Z n g µ n , g i .
Then, it holds that, almost surely,

max 1≤i≤n-1 b n i -λ n i , b n i -λ n i ≤ η, which implies that max b n 2 , b n 2 ≤ K 2 + √ n -1η.
Moreover, we have max 1≤i,j≤n-1

A n ij -δ ij ≤ η, which yields that for all ξ ∈ R n-1 , ( 1
-(n -1)η) ξ 2 2 ≤ ξ T A n ξ ≤ (1 + (n -1)η) ξ 2 2 .
Assume for now that η(n -1) < 1, this implies that, for all ξ ∈ R n-1 ,

(A n ) -1 ξ 2 ≤ 1 1 -(n -1)η ξ 2 . ( 31
)
Using (31), we obtain that

max λ n,n 2 , λ n,n 2 ≤ K 2 + √ n -1η 1 -(n -1)η . ( 32
)
We then have,

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ = Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µn (Z) - n-1 i=1 λ n,n i g i (Z) Z 1:∞ + Var g µn (Z) - n-1 i=1 λ n,n i g i (Z) Z 1:∞ -Var Z n g µn - n-1 i=1 λ n,n i g i + Var Z n g µn - n-1 i=1 λ n,n i g i -Var Z n g µ n - n-1 i=1 λ n,n i g i + Var Z n g µ n - n-1 i=1 λ n,n i g i -Var Z n g µ n - n-1 i=1 λ n i g i + Var Z n g µ n - n-1 i=1 λ n i g i -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ .
Using the fact that

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µn (Z) - n-1 i=1 λ n,n i g i (Z) Z 1:∞ ≤ 0, Var Z n g µn - n-1 i=1 λ n,n i g i -Var Z n g µ n - n-1 i=1 λ n,n i g i ≤ 0, Var Z n g µ n - n-1 i=1 λ n,n i g i -Var Z n g µ n - n-1 i=1 λ n i g i ≤ 0,
from the definition of λ n , λ n,n , λ n,n , µ n , we obtain that

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≤ Var g µn (Z) - n-1 i=1 λ n,n i g i (Z) Z 1:∞ -Var Z n g µn - n-1 i=1 λ n,n i g i + Var Z n g µ n - n-1 i=1 λ n i g i -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ , = Var g µn (Z)| Z 1:∞ -Var Z n (g µn ) -2 n-1 i=1 λ n,n i Cov g µn (Z), g i (Z)| Z 1:∞ -Cov Z n (g µn , g i ) + n-1 i,j=1 λ n,n i λ n,n j Cov g i (Z), g j (Z) Z 1:∞ -Cov Z n g i , g j + Var Z n g µ n -Var g µ n (Z) Z 1:∞ -2 n-1 i=1 λ n i Cov Z n g µ n , g i -Cov g i (Z), g j (Z) Z 1:∞ + n-1 i,j=1 λ n i λ n j Cov Z n g i , g j -Cov g µ n (Z), g i (Z) Z 1:∞ .
Now, using the definition of M n-1 given in (23), we obtain that

Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≤   1 + 2 n-1 i=1 | λ n,n i | + n-1 i,j=1 | λ n,n i || λ n,n j | + 1 + 2 n-1 i=1 |λ n i | + n-1 i,j=1 |λ n i ||λ n j |   × sup g,h∈Mn-1 Cov g(Z), h(Z)| Z 1:∞ -Cov Z n (g, h) .
Since sup

g,h∈Mn-1 Cov g(Z), h(Z)| Z 1:∞ -Cov Z n (g, h) = η, we then have, almost surely, Var g µn (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ -Var g µ n (Z) - n-1 i=1 λ n i g i (Z) Z 1:∞ ≤   1 + n-1 i=1 | λ n,n i | 2 + 1 + n-1 i=1 |λ n i | 2   sup g,h∈Mn-1 Cov g(Z), h(Z)| Z 1:∞ -Cov Z n (g, h) ≤ n 2 + n-1 i=1 | λ n,n i | 2 + n-1 i=1 |λ n i | 2 η ≤ n 2 + λ n,n 2 2 + λ n 2 2 η.
Finally, using (30) and (32), we obtain (29), i.e.

σ n-1 (M) -σ n-1 (M) ≤ n 2 + K 2 + √ n -1η 1 -(n -1)η 2 + K 2 2 η.
Let us now evaluate the probability, conditioned to the knowledge of Z 1:∞ , that

n 2 + K 2 + √ n -1η 1 -(n -1)η 2 + K 2 2 η ≤ (1 -γ 2 ) σ n-1 (M). ( 33
)
If η is chosen to be smaller that 1 2(n-1) , then it holds that

2 + K 2 + √ n -1η 1 -(n -1)η 2 + K 2 2 ≤ 2 + (2K 2 + 1) 2 + K 2 2 ≤ 9K 2 2 + 4.
A sufficient condition for (33) to hold is then to ensure that η ≤ with

:= min 1 2(n -1) , (1 -γ 2 ) σ 2 n-1 (M) n (9K 2 2 + 4)
, Then, it holds that

P G n Z 1:∞ = P σ n-1 (M) 2 ≥ γ 2 σ n-1 (M) 2 Z 1:∞ = P σ n-1 (M) 2 -σ n-1 (M) 2 ≤ (1 -γ 2 ) σ n-1 (M) 2 Z 1:∞ ≥ P n 2 + K 2 + √ n -1η 1 -(n -1)η 2 + K 2 2 η ≤ (1 -γ 2 ) σ n-1 (M) Z 1:∞ ≥ P η ≤ Z 1:∞ .
Thus, using the definition of η given by (28) and applying Lemma 3.12, we then obtain that

P G n |Z 1:(n-1) ≥ P η ≤ Z 1:∞ = P sup g,h∈M n-1 Cov g(Z), h(Z) Z 1:∞ -Cov Z n (g, h) ≤ Z 1:(n-1) ≥ 1 -δ n , since Ce -cMnφ(κn-1) ≤ δ n , with κ n-1 := min 1 2(n-1) , (1-γ 2 ) σ 2 n-1 (M) n(9K 2 2 +4) 6K n-1 ∞ K n-1 L ,
which yields the desired result.

We are now in position to end the proof of Theorem 3.6.

Proof of Theorem 3.6. Collecting Lemma 3.11 and Lemma 3.13, we obtain (18) for all n ∈ N * . Let us now prove (19). Let us first prove by recursion that for all n ∈ N * ,

P n k=1 G k ≥ Π n k=1 (1 -δ k ). (34)
Using Lemma 3.11, it holds that (34) is true for n = 1. Now we turn to the proof of the recursion. Let n ∈ N * . For any event Z, we denote by 1 Z the random variable which is equal to 1 if Z is realized and 0 if not. Using the fact that n k=1 G k is measurable with respect to Z 1:n , it holds that

P n+1 k=1 G k = E 1 n+1 k=1 G k = E E 1 Gn+1 1 n k=1 G k Z 1:n = E E 1 Gn+1 Z 1:n 1 n k=1 G k = E P G n+1 Z 1:n 1 n k=1 G k .
Now using Lemma 3.13, it holds that almost surely P G n+1 Z 1:n ≥ 1 -δ n+1 . Hence, it holds that

P n+1 k=1 G k ≥ (1 -δ n+1 )E 1 n k=1 G k = (1 -δ n+1 )P n k=1 G k .
The recursion is thus proved, together with (34), which implies (19).

If n∈N G n is realised, it then holds that the MC-greedy algorithm is a weak greedy algorithm with parameter γ and norm

• Z 1:∞ = Var •|Z 1:∞ .

Numerical results

The aim of this section is to compare several procedures to choose the value of the sequence (M n ) n∈N * in the MC-greedy algorithm presented in Section 3.1.

Three numerical procedures

As mentioned in Remark 3.9, it is possible to design a constructive way to define a sequence of numbers of samples (M n ) n∈N * which satisfies assumptions of Theorem 3.6, and thus which guarantees that the corresponding MC-greedy algorithm is a weak-greedy algorithm with high probability. Unfortunately, it appears that such a procedure cannot be used in practice to design a variance reduction method since the values of the sequence (M n ) n∈N * increases too sharply. The objective of this section is to propose a heuristic procedure to choose a sequence (M n ) n∈N * for an MC-greedy algorithm. This heuristic procedure appears to yield a reduced basis approximation f µ of f µ which provides very satisfactory results in terms of variance reduction, at least on the different test cases presented below.

We use here the same notation as in Section 2.1 and consider

M ref ∈ N * such that M ref 1.
The idea of this heuristic method is the following: assume that the sequence (M n ) n∈N * can be chosen so that for all n ∈ N * , the inequality

inf (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ - inf (λi) 1≤i≤n-1 ∈R n-1 Var Z n g µ - n-1 i=1 λ i g i ≤ (1 -γ 2 ) inf (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ (35)
holds for all µ ∈ P. Then, it can easily be checked that such an MC-greedy algorithm is a weak greedy algorithm with parameter γ. Of course, such an algorithm could not be of any use for variance reduction since it would imply the computation of

inf (λi) 1≤i≤n-1 ∈R n-1 Var g µ (Z) - n-1 i=1 λ i g i (Z) Z 1:∞ (or an approximation of this quantity of the form inf (λi) 1≤i≤n-1 ∈R n-1 Var Z ref g µ - n-1 i=1 λ i g i with Z ref = Z ref k 1≤k≤M ref a collection
of iid random variables with the same law as Z and independent of Z) for all µ ∈ P.

The idea of the heuristic procedure is then to check if the inequality (35) holds, only for the value µ = µ n . If the inequality holds, the value of M n+1 is chosen to be equal to M n for the next iteration. Otherwise, the value of M n is increased and the n th iteration of the MC-greedy algorithm is performed again. This heuristic procedure leads to the Heuristic MC-greedy algorithm (or HMC-greedy algorithm), see Algorithm 1. Notice that we introduce here P trial a finite subset of P, which is classically called the trial set of parameters in reduced basis methods.

For the sake of comparison, we introduce two other algorithms, which cannot be implemented in practice, but which will allow us to compare the performance of the HMC-greedy algorithm with ideal procedures. The first method, called SHMC-greedy algorithm and also presented in Algorithm 1 as a variant, consists in designing the sequence (M n ) n∈N * in order to ensure that the inequality (35) is satisfied for all µ ∈ P trial (and not only for µ n ). The second one consists in performing an ideal MC-greedy algorithm, called herefater IMC-greedy algorithm, see Algorithm 2, where all the variances and expectations are evaluated using M ref number of samples of the random variable Z at each iteration of the MC-greedy algorithm.

Let us comment on the termination criterion

Var Z ref f µ (S)H n-1 M ref > Var Z ref f µ (S)H n-1 -f µ (S)H n-1 M n-1
introduced in line 11 of the (S)HMC-greedy algorithm. Recall that, for µ = µ

(S)H n-1 , the expectation E f µ (S)H n-1
(Z) is approximated after n -1 iterations of the greedy algorithm by the control variate formula (see ( 1))

E Z ref (f µ (S)H n-1 ) + E Z n-1 f µ (S)H n-1 -f µ (S)H n-1 . ( 36
)
This criterion ensures that the iterative scheme ends when the statistical error associated with the second term in (36) becomes smaller than the statistical error of the first term (see Remark 2.1).

Figure 1 illustrates the evolution of the values of M n as a function of n for the HMC and SHMC algorithms. In comparison, in Figure 5, the relative error e (S)H n (µ) is plotted as a function of µ for n = 0, 5, 10, 20, 30. In particular, we observe that this error remains lower than 1% as soon as n ≥ 10 on P. Naturally, this error is larger for µ ∈ P test \ P. Finally, to illustrate the gain of our proposed method in terms of variance reduction, we plot on Figure 6 the value of the number of random Monte-Carlo samples M M C (µ) that would have been necessary to compute an approximation of the mean of f µ (Z) with a standard Monte-Carlo method with the same level of accuracy than the one given by the HMC-algorithm after n = 30 iterations. In this case, let us point out that M n = 349. More precisely, we compute M M C (µ) by the follwoing formula:

M M C (µ) = Var Z ref (f µ ) × M n Var Z n f µ - n i=1 λ µ i f µ i . (43)
Figure 6 illustrates that, for all µ ∈ P, the classical Monte Carlo method would have required a number of samples M M C (µ) in the range 10 6 ≤ M M C (µ) ≤ 10 12 in order to obtain the same level of statistical error. Thus, we see that the HMC-algorithm significantly improves the efficiency of the computation of the expectation of f µ (Z) with respect to a standard Monte-Carlo algorithm.

Second test case

In this example, we consider a second family of one-dimensional functions where P = [0, 1] is the set of parameter values. More precisely, we consider the family of functions (f µ ) µ∈P such that, for all µ in P:

∀x ∈ [0, 1], f µ (x) :=    √ x + 0.1 if x ∈ [0, µ] 1 2 (µ + 0.1) -1 2 x - 1 2 (µ + 0.1) -1 2 µ + (µ + 0.1) 1 2 if x ∈ [µ, 1] (44)
Let us point out that for all µ ∈ P, f µ is a C 1 function on [0, 1]. In this case, it is known [START_REF] Ehrlacher | Nonlinear model reduction on metric spaces. application to one-dimensional conservative pdes in wasserstein spaces[END_REF] that there exists a constant c > 0 such that d n (M) ≤ cn -2 for all n ∈ N * .

Let Z be a random variable with probability measure ν = U(0, 1).

In this example, M 1 = 10, M ref = 10 5 , γ = 0.9 and the trial set P trial is chosen to be a set of 300 random parameter values which were uniformly sampled in P. In this test case, we osbserve a similar behaviour of the (S)HMC-algorithm as in the first test case.

Figure 11 illustrates the computational gain brought by the HMC algorithm after n = 70 iterations (so that M n = 3109) with respect to the classical Monte Carlo method. Indeed, the quantity M M C (µ) defined in (43) is observed to vary in this case between 10 12 and 10 18 .

Two-dimensional heat equation

Let Z 1 and Z 2 be two independent real-valued random varibales with probability laws respectively U(0. We introduce a conform triangular mesh T of the domain D as represented on the left-hand side plot of Figure 13 and denote by V h := {u ∈ C (D) , u| T ∈ P 1 ∀T ∈ T , u| ∂D = 0} , the standard P 1 finite element space associated to this mesh.

For µ ∈ P and z ∈ (0.5, 2) × R, we define u µ,z h ∈ V h the unique solution to In this example, M ref = 10 5 , M 1 = 800 and γ = 0.9. Figure 14 illustrates the evolution of the values of M n as a function of n for the HMC and SHMC algorithms. It is to be noted here that the quantities θ H n , θ SH n and θ I n are very close: the quality of approximation of the reduced spaces V H n or V SH n is very close to the quality of approximation of the reduced space V I n given by an ideal greedy algorithm. Figure 17 shows the value of M M C (µ) given by (43), knowing that M n = 12800 after n = 7 iterations of the HMC algorithm. We observe that in this case 10 14 ≤ M M C (µ) ≤ 10 20 , which shows the huge computational gain brought by the HMC algorithm with respect to a standard Monte-Carlo method in this test case.

a µ,z (u µ,z h , v) = b(v), ∀v ∈ V h , (45)

  More precisely, let M small , M ref ∈ N * and assume that M ref M small . Let Z ref := Z ref k 1≤k≤M ref and Z small := Z small k 1≤k≤M small

  For any random functions g 1 , g 2 with values in L 2

Figure 1 :

 1 Figure 1: Evolution of M n as a function of n for the HMC and SHMC-greedy algorithms in test case 1.

Figure 2 :

 2 Figure 2: Evolution of θ H n (µ H n ), θ SH n (µ SH n ), θ H n , θ SH nas a function of n in test case 1.

Figure 2 ()

 2 Figure2illustrates the fact that at each iteration n ∈ N * , for the (S)HMC-algorithm, the value of the selected parameter µ (S)H n is relevant since we observe numerically that θ (S)H (µ(S)H n) is very close to

Figure 3 :

 3 Figure 3: Evolution of β H n , β SH n , θ H n , θ SH n , θ I n as a function of n in test case 1.

Figure 3

 3 Figure 3 illustrates the fact that the value of the number of samples M n chosen at each iteration n ∈ N * enables to compute empirical variances that are close to exact variances since the values of β (S)H n are very close to the θ (S)H n for the (S)HMC-algorithm.

Figure 4

 4 Figure 4: θ H n (µ), θ SH n (µ), β H n (µ), β SH n (µ) as a function of µ for n = 0, 5, 10, 20, 30 on P test = [0, 4].

Figure 5

 5 Figure 5: e H n (µ) and e SH n (µ) as a function of µ for n = 0, 5, 10, 20, 30 on P test = [0, 4].

Figure 6 :

 6 Figure 6: M M C (µ) as a function of µ ∈ P test = [0, 4].

Figure 7 :

 7 Figure 7: Evolution of M n as a function of n for the HMC and SHMC-greedy algorithms.

Figure 11 :Figure 12 :

 1112 Figure 11: M M C (µ) as a function of µ ∈ P test = [0, 1].

5 , 2 )

 52 and N (0, 1)) and let Z = (Z 1 , Z 2 ). Let D = (0, 2) 2 , P := [0, 10]. The trial set P trial is constructed by selecting 50 random values uniformly distributed in P.For all µ ∈ P and z := (z 1 , z 2 ) ∈ (0, 5, 2) × R, we introduceD µ,z : y) → D µ,z 11 (x, y) 0 0 D µ,z 22 (x, y) where ∀(x, y) ∈ D, D µ,z11 (x, y) = 13 + µ sin(2πx/z 1 ) + 0.5z 2 and D µ,z 22 (x, y) = 13 + µ sin(2πy/z 1 ) + 0.5z 2 .

where ∀v, w ∈ H 1 0Figure 13 :

 113 Figure 13: Left: mesh T (the triangle T 1 is highlighted in red color); Center: u µ,z h for µ = 9 and z = (1, 0); Right: u µ,z h for µ = 9 and z = (1.777, 0.2062).

Figure 14 :

 14 Figure 14: Evolution of M n as a function of n for the HMC and SHMC-greedy algorithms in test case 3.

Figure 15 :

 15 Figure 15: Evolution of θ H n(µ H n ), θ SH n (µ SH n ), θ H n , θ SHn as a function of n in test case 3.

Figure 16 :

 16 Figure 16: Evolution of β H n , β SH n , θ H n , θ SH n , θ I n as a function of n in test case 3.

Figure 17 :

 17 Figure 17: M M C (µ) as a function of µ ∈ P test = [0, 12].

Figure 18

 18 Figure 18: θ H n (µ), θ SH n (µ), θ I n (µ), β H n (µ), β SH n (µ) as a function of µ for n = 5 in test case 3.

Figure 19

 19 Figure 19: e H n (µ) and e SH n (µ) as a function of µ for n = 5 and n = 8 in test case 3.

Acknowledgements

The authors acknowledges financial support from the Mohammed VI Polytechnic University for funding Raed Blel's PhD thesis. Virginie Ehrlacher acknowledges support from project COMODO (ANR-19-CE46-0002). This work is partially supported by a funding from the European ResearchCouncil (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement No 810367), project EMC2.

Definitions of quantities of interest

For each of the test cases presented below, we plot different quantities of interest which we define here.

For all n ∈ N * , we denote by

) the set of parameter values selected after n iterations of the HMC-greedy (respectively SHMC-greedy and IMCgreedy) algorithm. We also denote by V

and

For all µ ∈ P and n ∈ N * , we define for the three algorithms presented in Section 4.1,

and

In what follows, we denote by θ H n (µ) and θ H n (respectively by θ SH n (µ), θ SH n , θ I n (µ) and θ I n ) the quantities defined by (37) and (38) obtained with the HMC-greedy (respectively the SHMC-greedy and IMC-greedy) algorithm. Note that, by definition of the IMC-greedy algorithm, θ I n = θ I n (µ I n ). A second quantity of interest for the HMC-greedy and the SHMC-greedy algorithms is given, for all n ∈ N * and µ ∈ P, by

and

In the sequel, we denote by β H n (respectively by β SH n ) the quantity defined by (40) obtained with the HMC-greedy (respectively the SHMC-greedy) algorithm.

Let us point out that when P trial = P and when M ref = ∞, it holds that θ n = σ n-1 (M) and θ n (µ n ) = σ n-1 (M) where σ n-1 (M) and σ n-1 (M) are defined respectively in [START_REF] Devore | Nonlinear approximation[END_REF] and [START_REF] Devore | The theoretical foundation of reduced basis methods. Model Reduction and approximation: Theory and Algorithms[END_REF].

We finally wish to evaluate the error made on the approximation of E[f µ (Z)] obtained by using the variance reduction method based on these MC-greedy algorithm. More precisely, this approximation is computed as output: N ∈ N * size of the reduced basis, µ

be a collection of M ref iid random variables with the same law as Z and independent of Z.

be a collection of M1 iid random variables with the same law as Z and independent of Z and Z ref .

5

Compute µ

).

8

• HMC case:

. Set n = 2 and Mn = M1.

11 while

Set R n = 1.

Let Z n := (Z n k ) 1≤k≤Mn be a collection of Mn iid random variables with the same law as Z and independent of Z and Z ref .

16

Compute µ

• SHMC case: For all µ ∈ P trial , compute θ SH n (µ) = min

, bn := max(1.1, rn) and set Mn = bnMn + 1 .

Algorithm 2: IMC-greedy algorithm input : > 0, P trial trial set of parameters (finite subset of P), M ref ∈ N * (high fidelity sampling number). output: N ∈ N * size of the reduced basis,

be a collection of M ref iid random variables with the same law as Z and independent of Z.

and µ i is equal to µ H i , µ SH i or µ I i depending on the chosen algorithm (remember formula [START_REF] Balabanov | Randomized linear algebra for model reduction. part i: Galerkin methods and error estimation[END_REF]). This quantity has to be compared with the approximation obtained with a standard Monte-Carlo with M ref samples, i.e. E Z ref (f µ ). To this aim, for all n ∈ N * and µ ∈ P, we define

In what follows, we denote by e H n (µ) (respectively e SH n (µ)) the quantity defined by (41) obtained by the HMC-greedy (respectively the SHMC-greedy) algorithm.

Explicit one-dimensional functions

We consider in this section two sets of one-dimensional explicit functions. The motivation for considering these two simple test cases is that the decays of the Kolmogorov n-widths of the associated sets M are known. Let P = [0, 3] be the set of parameter values. We consider in this first test case the family of functions (f µ ) µ∈P such that f µ (x) = f (x -µ) for all µ in P and x ∈ R. Let Z be a real-valued random variable with probability measure ν = U(0, 5). In this case, it is known [START_REF] Devore | Nonlinear approximation[END_REF] that there exists a constant c > 0 such that d n (M) ≥ cn -1/2 for all n ∈ N * .

First test case

In this example, M 1 = 10, M ref = 10 5 , γ = 0.9 and the trial set P trial is chosen to be a set of 300 random parameter values which were uniformly sampled in P.