Rémi Stroh

Julien Bect

Séverine Demeyer

Nicolas Fischer

Damien Marquis

Emmanuel Vazquez

Sequential design of multi-fidelity computer experiments: maximizing the rate of stepwise uncertainty reduction

Keywords: Multi-fidelity, Computer experiments, Sequential design of experiments, Gaussian process emulator, Meta-model, Surrogate model, Stochastic simulator

This article deals with the sequential design of experiments for (deterministic or stochastic) multi-fidelity numerical simulators, that is, simulators that offer control over the accuracy of simulation of the physical phenomenon or system under study. Accurate simulations usually entail a high computational effort, while coarse simulations are obtained at a lower cost. The cost can be measured, e.g., by the run time of the simulator or the financial cost of the computing resources. In this setting, simulation results obtained at several levels of fidelity can be combined in order to estimate quantities of interest (the optimal value of the output, the probability that the output exceeds a given threshold. . . ) in an efficient manner. We propose a new Bayesian sequential strategy called Maximal Rate of Stepwise Uncertainty Reduction (MR-SUR), that selects additional simulations to be performed by maximizing the ratio between the expected reduction of uncertainty and the cost of simulation. This generic strategy unifies several existing methods, and provides a principled approach to develop new ones. We assess its performance on several examples, including a computationally intensive problem of fire safety analysis where the quantity of interest is the probability of exceeding a tenability threshold during a building fire.

Introduction

In the domain of computer experiments, multi-fidelity refers to the idea of combining results from numerical simulations at different levels of accuracy: high-fidelity simulations correspond to more accurate but, in general, more expensive computations. As a representative example of a multi-fidelity simulator, consider the case of a partial differential equation (PDE) solver based on a finite element or finite volume method: the accuracy of the numerical solution depends, among other things, on the fineness of the discretization.

High-fidelity results are obtained when the mesh size is small. Conceptually, a simulator is viewed in this article as a black box with inputs and outputs. The parameter that controls the level of accuracy/fidelity-the mesh size in the case of a PDE solver-is one of the inputs of this black box, alongside others such as design or control variables and environmental variables (see, e.g., [START_REF] Santner | The Design and Analysis of Computer Experiments[END_REF]. Examples of multi-fidelity simulators can be found in virtually all areas of engineering and science, including aeronautics (Forrester et al., 2007), fire safety [START_REF] Demeyer | Surrogate model based sequential sampling estimation of conformance probability for computationally expensive systems: application to fire safety science[END_REF], electromechanics [START_REF] Hage Hassan | Radial output space mapping for electromechanical systems design[END_REF], electromagnetism [START_REF] Koziel | Robust variable-fidelity optimization of microwave filters using co-kriging and trust regions[END_REF], haemodynamics [START_REF] Perdikaris | Model inversion via multi-fidelity Bayesian optimization: a new paradigm for parameter estimation in haemodynamics, and beyond[END_REF]) and many more.

The cost of a simulation is often quantified by the run time of the simulator, or the financial cost of the computing resources, but other measures of cost can be more relevant in some contexts-for instance energy consumption or carbon emission in a green computing perspective (see, e.g., [START_REF] Kansal | Cloud load balancing techniques: A step towards green computing[END_REF]. We use henceforth the word "cost" as a generic term to refer to whichever measure of cost is relevant, in the application at hand, to express the available budget of simulations.

When the objective is to estimate a particular quantity of interest (QoI), such as the optimal value of the design variables (optimization problem) or the probability that the outputs belong to a prescribed "safe region" (reliability problem), multi-fidelity makes it possible to obtain a good approximation of the QoI with a computational effort lower than what would have been necessary if high-fidelity simulations only had been carried out. This cost reduction is achieved through the joint use of multi-fidelity models, which allow simulation results obtained at different levels of fidelity to be combined, and multi-fidelity designs of experiments (DoE). See Giselle [START_REF] Fernández-Godino | Review of multi-fidelity models[END_REF] for a review that covers both aspects. This article addresses the problem of constructing, sequentially, a multi-fidelity DoE targeting a given QoI.

We adopt a Bayesian point of view following the line of research initiated by [START_REF] Sacks | Design and analysis of computer experiments[END_REF], where prior belief about the simulator is modeled using a Gaussian process (see Santner et al., 2018, and references therein). The Bayesian approach provides a rich framework for the construction of sequential DoE, which has been abundantly relied upon in previous works dealing with the case of single-fidelity simulators, where the cost of a simulation is assumed to be independent of the value of the input variables (see, e.g., [START_REF] Kushner | A new method of locating the maximum point of an arbitrary multipeak curve in the presence of noise[END_REF][START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF][START_REF] Ranjan | Sequential experiment design for contour estimation from complex computer codes[END_REF][START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Picheny | Adaptive designs of experiments for accurate approximation of a target region[END_REF][START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]Chevalier et al., 2014). In this framework, sequential designs are usually constructed by means of a sampling criterion-also called infill criterion, acquisition function or merit function. The value of the criterion indicates whether a particular point in the input space is promising or not. The expected improvement (EI) criterion [START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF] is a popular example of such a sampling criterion. The extension of the Bayesian approach to sequential DoE in a multi-fidelity setting is based on two ingredients: 1) the construction of prior models for simulators with adjustable fidelity; 2) the construction of sampling criteria that take the variable cost of simulations into account.

For the case of deterministic multi-fidelity simulators, Gaussian process-based models have already been proposed in the literature [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF][START_REF] Qian | Building surrogate models based on detailed and approximate simulations[END_REF][START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF][START_REF] Gratiet | Multi-fidelity Gaussian process regression for computer experiments[END_REF]Le Gratiet and Garnier, 2014;Picheny and Ginsbourger, 2013;[START_REF] Tuo | Surrogate modeling of computer experiments with different mesh densities[END_REF]. Extensions to combinations of computer experiments and physical experiments [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF] and to stochastic simulators (Stroh et al., 2017a) have been considered as well.

Sampling criteria for single-fidelity sequential designs do not reflect a crucial feature of multi-fidelity simulators: the cost of a run depends on the value of the inputs (in particular on the one that controls the fidelity of simulation). Various methods that take into account a variable cost of simulation have been proposed for particular cases: single-objective unconstrained optimization [START_REF] Huang | Sequential kriging optimization using multiple-fidelity evaluations[END_REF][START_REF] Swersky | Multi-task Bayesian optimization[END_REF][START_REF] He | Optimization of multi-fidelity computer experiments via the EQIE criterion[END_REF] and global approximation [START_REF] Xiong | Sequential design and analysis of high-accuracy and low-accuracy computer codes[END_REF][START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF] notably, and more recently reliability analysis and contour estimation [START_REF] Marques | Contour location via entropy reduction leveraging multiple information sources[END_REF][START_REF] Marques | Multifidelity method for locating aeroelastic flutter boundaries[END_REF][START_REF] Yi | An active-learning method based on multi-fidelity kriging model for structural reliability analysis[END_REF].

In this article, we provide a general principled methodology to construct sequential DoE for multi-fidelity simulators and, more generally, for simulators where the cost of a simulation depends on the value of the inputs. The methodology is applicable to any QoI, and builds on the Stepwise Uncertainty Reduction (SUR) principle (see, e.g., [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF]Chevalier et al., 2014;Bect et al., 2019, and references therein), which unifies many of the aforementioned sequential DoE for the fixed-cost case.

More precisely, for the variable-cost case, we propose the Maximal Rate of Stepwise Uncertainty Reduction (MR-SUR) principle, which consists in constructing a sequential design by maximizing, at each step, the ratio between the expected reduction of uncertainty (to be defined more precisely later on) and the simulation cost.

The article is organized as follows. Section 2 reviews Gaussian process modeling for deterministic simulators, and discusses some possible extensions to (normally distributed) stochastic simulators. Section 3 first reviews both existing methods of sequential design for multi-fidelity simulators and the SUR principle for fixed-cost simulators, and then presents the MR-SUR principle and its relations with some existing sequential DoE. Finally, Section 4 illustrates the method and assesses its performance through several academic examples, including a computationally intensive problem of fire safety analysis where the quantity of interest is the probability of exceeding a tenability threshold during a building fire.

Gaussian process models for multi-fidelity

We consider a computer simulator with input variables u ∈ U ⊂ R d and one or several scalar outputs, which are generally obtained after some post-processing steps (e.g., an aerodynamic drag in a CFD model). Moreover, we consider that the accuracy, or fidelity, of the computer simulation can be tuned using a parameter δ that ranges in a discrete or continuous set T. For instance, δ can be the mesh size in a finite element or finite volume method. Such a parameter will be called fidelity parameter and can be viewed as an additional input of the simulator. We denote by x = (u, δ) ∈ X the aggregated vector of inputs, with X = U × T, and we assume from now on that the output is scalar.

The auto-regressive model for deterministic simulators

The so-called auto-regressive model of [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] assumes a deterministic simulator with a finite number S of levels of increasing fidelity. Let δ 1 , . . . , δ S denote the corresponding values of the fidelity parameter and set T = {δ 1 , . . . , δ S }. The simulator is then modeled by a Gaussian process ξ on X = U × T, defined through an auto-regressive relationship between successive levels:

     ξ(u, δ 1 ) = η 1 (u), ξ(u, δ s ) = ρ s-1 ξ(u, δ s-1 ) + η s (u), 1 < s ≤ S, (1)
where η 1 , . . . , η S are S mutually independent Gaussian processes, and (ρ s ) 1≤s<S ∈ R S-1 .

The model has been used in numerous applications, where the actual number S of levels is most often two (see, e.g., [START_REF] Qian | Building surrogate models based on detailed and approximate simulations[END_REF]Forrester et al., 2007;[START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF][START_REF] Kuya | Multifidelity surrogate modeling of experimental and computational aerodynamic data sets[END_REF][START_REF] Brooks | Multi-fidelity design optimisation of a transonic compressor rotor[END_REF][START_REF] Goh | Prediction and computer model calibration using outputs from multifidelity simulators[END_REF]Le Gratiet and Garnier, 2014;[START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF][START_REF] Elsayed | Optimization of the cyclone separator geometry for minimum pressure drop using co-kriging[END_REF][START_REF] Demeyer | Surrogate model based sequential sampling estimation of conformance probability for computationally expensive systems: application to fire safety science[END_REF], occasionally three (Perdikaris and Karniadakis, 2016, Section 3.2). In practice, the Gaussian processes η s are chosen among a family of Gaussian processes indexed by parameters such as correlation lengths, regularity parameters, etc., which are estimated from data (simulation results), by maximum likelihood for instance (see, e.g., [START_REF] Stein | Interpolation of Spatial Data: Some Theory for Kriging[END_REF]). Since the processes η s are assumed independent, there must be enough simulation results at each level of fidelity, even at the possibly very expensive highest fidelity levels, to obtain good estimates of the parameters-which explains perhaps why this model is typically used with a small number of levels.

Remark 1. Some authors introduce more flexibility in the model by letting the coefficients ρ s in Equation (1) depend on u as well, for instance linearly [START_REF] Qian | Building surrogate models based on detailed and approximate simulations[END_REF]Le Gratiet and Garnier, 2014) or using another Gaussian process in a hierarchical manner [START_REF] Qian | Bayesian hierarchical modeling for integrating low-accuracy and high-accuracy experiments[END_REF].

The additive model for deterministic simulators

Another approach to building Gaussian process models for deterministic multi-fidelity simulators, which readily applies to the case where T contains a continuum of levels of fidelity or a large number of ordered discrete levels of fidelity, has been proposed by Picheny and Ginsbourger (2013) and [START_REF] Tuo | Surrogate modeling of computer experiments with different mesh densities[END_REF].

Assuming for simplicity that T = [0, ∞), with δ = 0 corresponding as in [START_REF] Tuo | Surrogate modeling of computer experiments with different mesh densities[END_REF] to the highest-often unreachable-level of fidelity, a Gaussian process ξ over the product space X = U × T is defined in this approach as the sum of two independent parts:

ξ(u, δ) = ξ 0 (u) + ε(u, δ), ( 2
)
where ξ 0 and ε are mutually independent Gaussian processes, and ε has zero mean and goes to zero in the mean-square sense when δ → 0. In other words, var (ε(u, δ)) → 0 when δ → 0, for all u: as a consequence, ξ is a non-stationary Gaussian process on X ⊂ R d+1 . Under this decomposition, ξ 0 represents the "ideal" version of the simulator, while ε represents numerical error.

In both articles, ξ 0 is then assumed to be stationary, whereas the covariance function of ε is multiplicatively separable: for all u, u ′ ∈ U and δ, δ ′ ∈ T,

cov (ε(u, δ), ε(u ′ , δ ′ )) = r(δ, δ ′ ) k(u, u ′ ), ( 3
)
where k is a stationary covariance function on U, and r is a non-stationary covariance function on T such that r(δ, δ) → 0 when δ → 0. As an example of a suitable choice for r, consider the Brownian-type model proposed by [START_REF] Tuo | Surrogate modeling of computer experiments with different mesh densities[END_REF]:

r(δ, δ ′ ) = min{δ, δ ′ } L , ( 4
)
with L a positive parameter. Other choices are of course possible.

Extension to stochastic simulators

We turn now to the case of stochastic simulators, that is, simulators whose output is stochastic, as happens for instance when the computer program relies on a Monte Carlo method (see, e.g., [START_REF] Cochet | Capabilities overview of the MORET 5 Monte Carlo code[END_REF]. Extending the multi-fidelity Bayesian methodology of Sections 2.1 and 2.2 to stochastic simulators is not straightforward in general, since the output at a given input point x i = (u i , δ i ) ∈ X is now a random variable Z i , the distribution of which is in general unknown and different at each point in X. (Several runs at the same input point yield independent and identically distributed responses.)

We focus in this section on the simpler case where the output Z i can be assumed to be normally distributed:

Z i | ξ, λ ∼ N (ξ(x i ), λ(x i )), (5)
with mean ξ(x i ) and variance λ(x i ) possibly depending on the input point. In this setting, we propose to extend the multi-fidelity models of previous sections using independent prior distributions for ξ and λ, with either the autoregressive model of Section 2.1 or the additive model of Section 2.2 as a prior for ξ. Then, since the variance λ must have positive values, we suggest modeling log(λ) by a Gaussian process λ, following [START_REF] Goldberg | Regression with inputdependent noise: A Gaussian process treatment[END_REF], [START_REF] Kersting | Most likely heteroscedastic Gaussian process regression[END_REF], [START_REF] Boukouvalas | Learning heteroscedastic Gaussian processes for complex datasets[END_REF] and others.

Under this type of model, the inference task-estimating the parameters of the Gaussian process models for ξ and λ, and then computing posterior distributions-becomes more difficult since neither ξ nor λ are directly observable. [START_REF] Goldberg | Regression with inputdependent noise: A Gaussian process treatment[END_REF] take a fully Bayesian approach using a time-consuming Monte-Carlo method. Other authors suggest optimization-based approaches, that simultaneously produce estimates the parameters of of both Gaussian processes and of the unobserved log-variances. In particular, [START_REF] Kersting | Most likely heteroscedastic Gaussian process regression[END_REF] and [START_REF] Boukouvalas | Learning heteroscedastic Gaussian processes for complex datasets[END_REF] propose a method called most likely heteroscedastic GP, stemming from the Expectation-Maximization (EM) algorithm (see also [START_REF] Marrel | Global sensitivity analysis of stochastic computer models with joint metamodels[END_REF], for a similar algorithm), while [START_REF] Binois | Practical heteroscedastic Gaussian process modeling for large simulation experiments[END_REF] use a more sophisticated joint maximization procedure with relaxation to obtain the joint MAP (maximum a posterior) estimator.

For the numerical experiments of this article (Sections 4.3 and 4.4) we will take a simpler route, assuming that the variance λ depends only on the fidelity level δ-which is approximately true in the two examples we shall consider. In this setting, as long as the number of fidelity levels of interest is not too large, the value of the variance at these levels can be simply estimated jointly with the other parameters of the model; a general-purpose log-normal prior for the vector of variances is proposed by [START_REF] Stroh | Integrating hyperparameter uncertainties in a multi-fidelity Bayesian model for the estimation of a probability of failure[END_REF].

3 Sequential design of experiment for multi-fidelity

Existing methods

In the literature of multi-fidelity, a variety of sequential design algorithms have been proposed. (See Supplementary Material for a review of nonsequential multi-fidelity designs, which can be used as initial designs for sequential ones.)

For instance, Forrester et al. (2007) suggest using the auto-regressive model of [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] and a standard single-level sequential design at the highest level of fidelity to select input variables u ∈ U for the next experiment. Then, simulations at all levels of fidelity are run for the selected u. Building on Forrester et al. (2007), [START_REF] Kuya | Multifidelity surrogate modeling of experimental and computational aerodynamic data sets[END_REF] suggest a two-stage method: run a large number of simulations at the low-fidelity level, and then use a sequential design strategy to select simulations at the high-fidelity level. In a different spirit, [START_REF] Xiong | Sequential design and analysis of high-accuracy and low-accuracy computer codes[END_REF] use Nested Latin Hypercube Sampling (NLHS) and suggest to double the number of simulations when going from level δ (s) to level δ (s+1) , until some cross-validation-based criterion is satisfied.

More interestingly in the context of this article, some methods have been proposed

that explicitly take into account the simulation cost. This is typically achieved by crafting a sampling criterion that takes the form of a ratio between a term which measures the interest of a simulation at (u, δ), and the cost of the simulation [START_REF] Huang | Sequential kriging optimization using multiple-fidelity evaluations[END_REF][START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF][START_REF] He | Optimization of multi-fidelity computer experiments via the EQIE criterion[END_REF][START_REF] Marques | Contour location via entropy reduction leveraging multiple information sources[END_REF][START_REF] Marques | Multifidelity method for locating aeroelastic flutter boundaries[END_REF][START_REF] Yi | An active-learning method based on multi-fidelity kriging model for structural reliability analysis[END_REF].

For instance, [START_REF] He | Optimization of multi-fidelity computer experiments via the EQIE criterion[END_REF] propose a global optimization method using the Expected Quantile Improvement (EQI) of Picheny et al. (2013) and the multi-fidelity model of [START_REF] Tuo | Surrogate modeling of computer experiments with different mesh densities[END_REF], and build a new sampling criterion corresponding to the ratio between the EQI sampling criterion and the simulation cost.

Outside the multi-fidelity literature, a similar idea has been proposed by [START_REF] Johnson | An information theory approach to diagnosis[END_REF] to design sequential testing procedures, by [START_REF] Swersky | Multi-task Bayesian optimization[END_REF] for multi-task optimization and by [START_REF] Poloczek | Multi-information source optimization[END_REF] for multi-information source optimization. In the first two cases, the numerator of the criterion is the expected reduction of the entropy of the QoI, while in the third one it is equal to the Knowledge Gradient criterion [START_REF] Frazier | The knowledge-gradient policy for correlated normal beliefs[END_REF].

In this article, we propose a general methodology to build such sequential designs, which is not tied to a particular kind of model or QoI. The key idea is to measure the potential of a particular design point using the SUR framework, recalled in Section 3.2. The methodology itself, that we call MR-SUR, is presented in Section 3.3.

Stepwise uncertainty reduction

We recall here the principle of SUR strategies, introduced in the design of computer experiments by Vazquez and co-authors (see, e.g., [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF][START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF].

Given a Bayesian model of a simulator and an unknown QoI Q, that is, a particular feature of the simulator that we want to estimate, a SUR strategy is a Bayesian method for the construction of a sequence of evaluation locations X 1 , X 2 , . . . ∈ X at which observations of the simulator will be taken in order to reduce the uncertainty on Q. (In this section, X denotes a generic input space, not necessarily of the form X = U × T.)

The starting point of the construction of a SUR strategy is the definition of a statistic H n measuring the residual uncertainty about Q given past observations Z 1 , . . . , Z n . Many choices for H n are possible for any particular problem, but a natural requirement [START_REF] Bect | A supermartingale approach to Gaussian process based sequential design of experiments[END_REF] is that H n should be decreasing on average when n increases. For instance,

H n can be the posterior Shannon entropy of Q for a discrete QoI, or the posterior (total)

variance for a scalar (or vector-valued) QoI. If Q is a function defined on X, as will be the case in Section 4, a possible choice is

H n = E n Q -Q n 2 µ = X var n (Q(x)) µ(dx), (6)
where µ denotes a measure X,

h 2 µ = X h(x) 2 µ(dx), E n (resp. var n ) is the posterior expectation (resp. variance) given Z 1 , . . . , Z n , and Q n (x) = E n (Q(x)).
Then, given past observations, X n+1 is chosen by minimizing the expectation of the future residual uncertainty:

X n+1 = argmin x∈X J n (x), with J n (x) = E n (H n+1 | X n+1 = x) , ( 7
)
where the expectation is with respect to the outcome Z n+1 of a new simulation at x ∈ X.

Similarly, a batch of q evaluation points can be chosen by minimizing a multi-point criterion:

J n,q (x 1 , . . . , x q ) = E n (H n+q | X n+1 = x 1 , . . . , X n+q = x q ) (8)
(see, e.g., Chevalier et al., 2014;[START_REF] Azzimonti | Adaptive design of experiments for conservative estimation of excursion sets[END_REF].

Example. Assume a stochastic multi-fidelity simulator defined over X = U × T as in Section 2.3, and consider the functional QoI defined on X by

Q(x) = P Z x > z crit ξ, λ = Φ   ξ(x) -z crit λ(x)   , ( 9
)
where Z x denotes the outcome of a new simulation at x, z crit ∈ R is a given threshold, and Φ the cdf of the standard normal distribution. Pick some reference level δ ref ∈ T and consider the residual uncertainty

H n = U var n Q(u, δ ref ) du, (10)
which is a special case of Equation ( 6) with µ equal to Lebesgue's measure on U at fixed δ = δ ref .

Then, using computations similar to those of Chevalier et al. (2014), it can be proved that

J n (x) = U Φ 2 a n (x ′ ), a n (x ′ ); k n (x ′ , x ′ ) v n (x ′ ) -Φ 2 a n (x ′ ), a n (x ′ ); k n (x, x ′ ) 2 v n (x)v n (x ′ ) du ′ , ( 11
)
where

x ′ = (u ′ , δ ref ), m n (resp. k n ) denotes the posterior mean (resp. covariance) of ξ, v n (x) = λ(x) + k n (x, x), a n (x) = (m n (x) -z crit ) / v n (x)
, and Φ 2 (•, • ; ρ) is the cdf of the standard bivariate normal distribution with correlation ρ. (For tractability, the variance function λ is assumed to be known in the computation of the criterion. In practice, the estimated variance function is plugged in the expression, and the integral over U is approximated using a Monte Carlo method.)

Remark 2. See Supplementary Material for a proof of Equation ( 11), in a more general form which also allows for (possibly parallel) batches of evaluations as in ( 8), and integration with respect to a generic measure µ as in Equation ( 6).

Remark 3. In the special case λ ≡ 0 (deterministic simulator), corresponding to Q(x) = 1 ξ(x)>z crit , the criterion (11) has been proposed by [START_REF] Bect | Sequential design of computer experiments for the estimation of a probability of failure[END_REF] and computed by Chevalier et al. (2014). The general case, to the best of our knowledge, is new.

The reader is referred to, e.g., [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF]

Maximum rate of stepwise uncertainty reduction

The proposed Maximum Rate of Stepwise Uncertainty Reduction (MR-SUR) strategy builds on the SUR strategy presented in Section 3.2. The goal is to achieve a balance between the (expected) reduction of uncertainty brought by new observations on the one hand, and the cost of these observations, usually measured by their computation time, on the other hand. Denote by C : X → R + the cost of an observation of the simulator, which in general depends mainly on the fidelity level δ ∈ T, but can also depend on the input u ∈ U in some applications. The MR-SUR strategy is then given by

X n+1 = argmax x∈X H n -J n (x) C(x) = argmax x∈X G n (x) C(x) , ( 12
)
where G n (x) = H n -J n (x) is the expected uncertainty reduction associated to a future observation at x ∈ X. This strategy boils down to a SUR strategy when C is constant.

A few special cases of MR-SUR strategies, adapted to particular models and estimation goals, have been proposed earlier in the literature. To the best of our knowledge, the oldest example is the sequential testing method of [START_REF] Johnson | An information theory approach to diagnosis[END_REF], where H n is the posterior to the expected improvement [START_REF] Mockus | The application of Bayesian methods for seeking the extremum[END_REF][START_REF] Jones | Efficient global optimization of expensive black-box functions[END_REF], IAGO [START_REF] Villemonteix | An informational approach to the global optimization of expensive-to-evaluate functions[END_REF] and KG [START_REF] Frazier | The knowledge-gradient policy for correlated normal beliefs[END_REF] algorithms-see [START_REF] Bect | A supermartingale approach to Gaussian process based sequential design of experiments[END_REF] for a presentation of EI and KG in the SUR framework. The first sequential design procedure of Le [START_REF] Le Gratiet | Cokriging-based sequential design strategies using fast cross-validation techniques for multi-fidelity computer codes[END_REF] can also be seen as an approximate MR-SUR strategy for the approximation problem, where H n is the posterior integrated prediction variance. More recently, the CLoVER algorithm of [START_REF] Marques | Contour location via entropy reduction leveraging multiple information sources[END_REF][START_REF] Marques | Multifidelity method for locating aeroelastic flutter boundaries[END_REF] provides an illustration of the MR-SUR principle for the problem of contour estimation.

Remark 4. Assuming that the cost depends only on the fidelity level δ, the MR-SUR principle extends very naturally to batches of evaluations with a common cost C(δ), replacing the expected uncertainty J n (x) in the numerator by J n,q (x 1 , . . . , x q ), where q is the size of the batch (as in Equation ( 8)). The application of this idea to (synchronous) parallel computations-assuming that C(δ) is the computation time-is left for future work.

To illustrate the MR-SUR principle, let us consider a simple simulated example, with ξ a Gaussian process on

X = U × T = [-0.5, 0.5] × [0, 1] such that ξ | m ∼ GP(m, k),
m ∼ U(R), and k as in Section 2.2:

k((u, δ), (u ′ , δ ′ )) = σ 2 0 M ν 0 |u -u ′ | ρ 0 + σ 2 0 G min {δ, δ ′ } L M νε |u -u ′ | ρ ε ,
where 9)-( 10), with z crit = 0. Note that the level of highest fidelity δ = 0 is not observable in practice. A NLHS of size n = 12 + 6 + 6 + 3 on the levels δ = 1, 1/2, 1/5, 1/10 is taken as the observed DoE, and the outputs Z 1 , . . . , Z n are simulated according to Equation (5) with constant variance λ = 0.4 2 . We compute the functions J n , G n and C over a 100 × 100 regular grid on U × T, to obtain Figures 1a and1b.

Observe on Figure 1a that, for each cost value (corresponding to a fixed level of fidelity), there is a range of points that yield different values of the expected uncertainty reduction.

Good observation points lie on the Pareto front (in solid black line), that is, the set of points for which there is no larger expected uncertainty reduction at lower cost. The MR-SUR strategy selects an observation location that correspond to the maximum of the "slope" of this Pareto front.

Figure 2 shows the sequence of Pareto fronts as more observation points are added in the design using Equation ( 11). The horizontal axis now represents the total cost of the design, including the cost of the candidate point, so that the left-ends of the Pareto fronts are shifted (the top-left curve corresponds to the situation depicted on Figure 1a). This figure illustrates how the MR-SUR strategy progressively reduces the uncertainty about the QoI, by selecting at each step a pair (u, δ) that makes a trade-off between cost and expected uncertainty reduction.

Remark 5. In the numerical experiments of Section 4, it will be assumed that the cost function C(x) is known beforehand. When it is not, it is usually estimated prior to the sequential design phase using a combination of expert knowledge and pilot runs (as described, e.g., by [START_REF] He | Optimization of multi-fidelity computer experiments via the EQIE criterion[END_REF]. Alternatively, it could be estimated during the sequential design phase, using for instance (as suggested by [START_REF] Swersky | Multi-task Bayesian optimization[END_REF] a second GP model for x → log C(x).

Numerical results

total cost

J n (x)

Setup of the experiments

In each example, we consider a multi-fidelity simulator for which the simulation cost C depends on δ alone, and is assumed to be known. Some common features of all three numerical experiments are presented in this section.

Initial DoE.

A nested Latin hypercube sample (NLHS) is used as an initial design.

More specifically, we use the algorithm developed by [START_REF] Qian | Nested latin hypercube designs[END_REF], with an additional maximin optimization at each level to obtain better space-filling properties (see Stroh, 2018, Section 2.2.3 for details).

GP modeling. In each example, a multi-fidelity GP model of the type described in Section 2 is used. The posterior distribution of the parameters is initially sampled using an adaptive Metropolis-Hastings algorithm [START_REF] Haario | An adaptive Metropolis algorithm[END_REF] and then updated at each iteration by sequential Monte Carlo (see, e.g., [START_REF] Chopin | A sequential particle filter method for static models[END_REF]. More details about the particular GP model that is used, and the prior distribution on the parameters, are provided inside each example section.

Optimization of the sampling criterion. At each iteration of a SUR or MR-SUR strategy, a new simulation point is selected according to Equation (7) or Equation ( 12). This step involves an optimization of the SUR or MR-SUR criterion, which is carried out in the experiments of this article using a simple two-step approach: the criterion is first optimized by exhaustive search on a regular grid over U×T, and then a local optimization is performed starting from the best point in the grid. Other approaches have been proposed in the literature, that would be more efficient in higher-dimensional problems (see, e.g., [START_REF] Feliot | A Bayesian approach to constrained singleand multi-objective optimization[END_REF].

Other computational details. All integrals are approximated by Monte-Carlo methods.

SUR and MR-SUR criteria are evaluated using the Maximum A Posteriori (MAP) estimator of the parameters-obtained by local optimization from the best point in the MCMC/SMC sample-in a plug-in manner. (A fully Bayesian approach could be considered in principle, but would lead to much higher computational complexity.)

A one-dimensional example

Consider as a first (toy) example the two-level deterministic simulator defined for u ∈ [0; 1] and δ ∈ {1, 2} by the analytical formulas (Forrester et al., 2007)

     f 1 (u) = f (u, 1) = 0.5 (6u -2) 2 sin(12u -4) + 10 (u -0.5), f 2 (u) = f (u, 2) = (6u -2) 2 sin(12u -4) + 10, (13)
and assume that computing f 2 -hereafter referred to as the "high fidelity" function-is four times more costly than computing f 1 : C(2) = 1 and C(1) = 1 4 . Note that the two functions are related by f 2 (u) = 2 f 1 (u) -20(u -1), which makes them perfect candidates for the autoregressive model presented in Section 2.1. The goal in this example is to estimate the set

Γ = f 2 > z crit = u ∈ U, f 2 (u) ≥ z crit
with z crit = 10. The performance of MR-SUR for this task will be compared with that of SUR strategies operating at the low-fidelity level only (LF-SUR) or at the high-fidelity level only (HF-SUR).

In this experiment, all three sequential strategies start with the same multi-fidelity initial design, and use the same Gaussian process prior and measure of uncertainty H n .

The initial design consists of six observations at the low-fidelity level and three at the high-fidelity level, for a total of n = 9 observations, corresponding to an initial budget of 6 × 1 4 + 3 × 1 = 4.5 cost units. A supplementary budget of 9.0 cost units is assumed to be available for the sequential design. The autoregressive model of Section 2.1 is used, with Matérn 5/2 covariance functions and weakly informative priors on the parameters (see Supplementary Material for details). The uncertainty on Γ is quantified using the uncertainty measure (10). In the special case of a deterministic simulator, we have (cf.

Remark 3) where p n (u) = P n (ξ 2 (u) ≥ z crit ) is the posterior mean of 1 ξ 2 (u)≥z crit , and

H n = 1 0 p n (u) (1 -p n (u)) du,
p n (u) (1 -p n (u))
its posterior variance.

The experiment is repeated R = 60 times-the simulator is deterministic, but randomness in the result comes from both the initial DoE and the use of a Monte Carlo procedure to sample from the posterior of the parameters. Figure 3a presents the evolution of the median estimation error, defined as MedErr n = median 1≤r≤R p (r) n -1 f 2 >z crit with • the L 2 -norm on U, as function of the cost c n = i≤n C(δ i ). First, it appears clearly that highfidelity evaluations are needed: the LF-SUR strategy achieves no significant error reduction with respect to the initial design. Second, we observe that the combination of low-and high-fidelity evaluations chosen by the MR-SUR strategy is more efficient, on average, than a purely high-fidelity sequential design. The actual number of evaluations on each level is summarized, for the 60 repetitions, in the table in Figure 3b: the MR-SUR strategy tends to use between two and five high-fidelity evaluations. (The recommendation of [START_REF] Xiong | Sequential design and analysis of high-accuracy and low-accuracy computer codes[END_REF]-observing the low-fidelity level twice as many times as the high-fidelity one-would correspond here to six high-fidelity evaluations.)

Random damped harmonic oscillator

We now assess the performance of MR-SUR on an example proposed by [START_REF] Au | Estimation of small failure probabilities in high dimensions by subset simulation[END_REF]. We consider a random damped harmonic oscillator, whose displacement Y is the Level δ 1 s 0.5 s 1/3 s 0.25 s 0. solution of the stochastic ordinary differential equation

Ÿ (t) + 2ζω 0 Ẏ (t) + ω 2 0 Y (t) = W (t), t ∈ [0, t end ], Ẏ (0) = 0, Y (0) = 0 , ( 14
)
where ω 0 is the resonance frequency of the oscillator, ζ is a damping coefficient, W is a Gaussian white noise and t end = 30 s. The solution of Equation ( 14) is approximated using an exponential Euler scheme with time step δ > 0 (more details can be found in the Supplementary Material): we denote by

Y (δ) k
the resulting approximation of Y at time steps

t k = kδ, k ∈ N, k ≤ K δ = ⌊t end /δ⌋.
We will be interested in the maximal log-displacement

max t≤t end log |Y (t)|, that we approximate by Z(ω 0 , ζ, δ) = max k≤K δ log |Y (δ)
k | . We view the random variable Z(ω 0 , ζ, δ) as the response of a multi-fidelity stochastic simulator, where δ controls the level of fidelity and u = (ω 0 , ζ) ∈ R 2 + is the vector of input variables. In this problem, the QoI is the function Q defined by A good approximation of the output distributions is obtained if we assume

Q(ω 0 , ζ) = P(Z(ω 0 , ζ, δ ref ) > z crit ), 0 ≤ ω 0 ≤ 30, 0 ≤ ζ ≤
Z(ω 0 , ζ, δ) | ξ, λ ∼ N (ξ(ω 0 , ζ, δ ref ), λ(δ))
, where the variance only depends on the fidelity level. This assumption makes it possible to write

Q(ω 0 , ζ) = Φ   ξ(ω 0 , ζ, δ ref ) -z crit λ(δ ref )   .
The mean function ξ is modeled by the additive Gaussian process model ( 2)-( 4) of Section 2.2, where the variance λ is log-Gaussian as in Section 2.3 and the prior distributions for the parameters are set as in [START_REF] Stroh | Integrating hyperparameter uncertainties in a multi-fidelity Bayesian model for the estimation of a probability of failure[END_REF]. The posterior mean

Q n (ω 0 , ζ) = E n (Q(ω 0 , ζ)) is used to estimate Q.
Remark 6. In this particular application, it was found that the value c = 0.99 proposed by [START_REF] Stroh | Integrating hyperparameter uncertainties in a multi-fidelity Bayesian model for the estimation of a probability of failure[END_REF], for the correlation between logarithms of noise variances, gave too much influence to the levels of lowest fidelity. The value c = 0.9 was used instead. In this example we consider S = 10 levels, and the initial design is an NLHS on the five first levels. The different levels of fidelity, their costs, and the initial design are summarized in Table 1. The total cost of the initial design is 9.88. The total simulation budget, taking into account the initial budget, is set to 20. We also use a very high simulation budget to compute a reference value for Q, which will be used to assess the estimation error.

We compare the MR-SUR strategy, using the integrated posterior variance (10) as the uncertainty measure, to five different SUR strategies based on the same uncertainty measure. All of them start with the same multi-fidelity initial design (see Table 1). Each SUR strategy corresponds to sampling on only one of the five highest-fidelity levels, using this particular fidelity level both in the definition of the uncertainty measure (10) and for the prediction of the QoI. The experiment is repeated 200 times with different initial designs, and the strategies are compared using the median estimation error as in Section 4.2. A regular grid of size 100 × 100 is used to discretize the integral in Equation ( 10) and the L 2 norm of the error. The result is shown in Figure 4.

Comparing first MR-SUR with the fixed-level SUR strategy targetting the same level (namely, δ = δ ref = 0.01), it can be seen that MR-SUR makes a good use of the additional degree of freedom that multi-fidelity brings. Indeed, MR-SUR achieves superior performance for any budget between 10 and 16, and then performs similarly to the fixed-level strategy. Figure 5 shows some examples of the sequence of levels selected by MR-SUR :

the design typically explores the level δ = 0.1 first, and then gradually moves to levels of higher fidelity.

Turning now to the results of the other fixed-level SUR strategies, it is interesting to see that the results depends very much on the fidelity level that is chosen. MR-SUR is by far superior to the fixed-level strategies at δ = 0.16 or δ = 0.17, for any budget, cost cost cost but is outperformed by the fixed-level strategy at δ = 0.05 and, for budgets larger than approximately 14, also by the one at δ = 0.02. The reason for this lies in two facts: First, the function (ω 0 , ζ) → P(Z(ω 0 , ζ, δ) > z crit ) at these levels is close enough to the highfidelity one. Second, the variance of the noise is cheaper to estimate at these levels. The MR-SUR strategy, as implemented here, fails to figure this out automatically since the uncertainty about the variance of the output is not taken into account.

δ δ δ 1/4 1/4 1/4 1/5 1/5 1/

A fire safety example

In this section, we illustrate the MR-SUR strategy on a fire safety application. The goal is to assess the safety of a 20 m × 12 m × 16 m parallelepiped-shaped storage facility, with two 2 m × 1 m open doors and two 2 m × 2 m open windows. The propagation of smoke and heat is simulated using Fire Dynamics Simulator (FDS; see [START_REF] Mcgrattan | Fire dynamics simulator (version 5), user's guide[END_REF], a state-of-the-art CFD software for fire engineering, which solves the transport equations using finite difference methods. The fire is located at the center of the room, and burns polyurethane. To assess fire safety, the values of several physical quantities are compared against regulatory thresholds-in this illustration, we focus on one of them only, called visibility and hereafter denoted by V . According to the ISO 13571 standard (2012), visibility must remain greater than z crit = 5 m to ensure safety during an evacuation.

Our FDS-based simulator will be treated as a stochastic simulator with nine input variables: three environmental variables (external temperature T ext , atmospheric pressure P atm and ambient temperature T amb ) denoted by u e ∈ R 3 , five "scenario variables" (fire growth rate α, fire area A f , maximal heat release rate Q′′ h , total released energy per unit area q f d and soot yield Y soot ) denoted by u s ∈ R 5 , and finally the size δ of the spatial discretization mesh, which plays the role of a fidelity parameter. The reader is referred to Stroh et (2017a) and [START_REF] Stroh | Sequential design of numerical experiments in multi-fidelity : Application to a fire simulator[END_REF] for more details on the application.

Remark 7. Strictly speaking, the simulator considered in this section is actually deterministic, since the seed of the random number generator is fixed by the FDS software; cf. Stroh et al. (2017a) for details.

In this example, our objective is to estimate the probability that V becomes less than z crit in a particular fire scenario, defined by α = 0.1057 kW Four levels of fidelity will be considered for running simulations: the reference level δ ref = 20 cm, and three levels of lower fidelity (δ = 50 cm, 33 cm and 25 cm). Table 2 shows the correspondence between levels and computation times, assuming for simplicity that the computation time does not depend on the input variables (the fire area A f and the total released energy per unit area q f d , which together control the amount A f q f d of combustible material to be burnt, actually have some influence on the computation time, which is neglected). Four independent initial NLHS designs of size n = 130, distributed across the first three levels of fidelity as shown in Table 2, are available from previous studies. The normalized cost of each initial design is 9.89 (i.e., 20 days). A reference value for Q has also been obtained from 150 Monte Carlo simulations, distributed on the highest fidelity level using P ue . This reference value has a normalized cost of 150 (i.e., about ten months).

• s -2 , A f = 14 m 2 , Q′′ h = 460 kW • m -2 , q fq = 450 MJ • m -2 ,
We run the MR-SUR strategy starting from our four initial designs, using a supplementary budget of 24 for each run (about 48.7 days). The underlying Bayesian model is the same as in Section 4.3, the QoI Q is estimated using the posterior mean Qn = E n (Q), and the measure of uncertainty is the integrated posterior variance which is a special case of Equation ( 6). The corresponding SUR criterion is similar to (11), with an integral over the environmental variables only. The result is shown in Figure 6.

We can see on Figure 6a that the estimations are initially, in three out of four cases, incompatible with the Monte Carlo one, but tend to get closer to the reference value when more simulations are carried out using the MR-SUR strategy. Figure 6b shows the measure of uncertainty as a function of the cost: the uncertainty is large at the beginning of the sequential design and rapidly becomes smaller, as expected, as the MR-SUR strategy proceeds. (Note that the cost of the whole design is approximately 9.9 + 24 = 33.9, which is must cheaper than the cost of 150 of the Monte Carlo reference.)

Conclusion

The main contribution of this article is to unify and extend several methods of the literature of Bayesian sequential design of experiments for multi-fidelity numerical simulators. The unification that we propose is cast in the framework of Stepwise Uncertainty Reduction (SUR) strategies: when the accuracy of computer simulations can be chosen by the user, a natural extension of SUR strategies is to consider sampling criteria built as the ratio between the expected reduction of uncertainty and the simulation cost. We call this approach Maximal Rate of Stepwise Uncertainty Reduction (MR-SUR). It can be applied to deterministic or stochastic simulators. Our numerical experiments show that the MR-SUR approach typically provides estimations which, for a given computational cost, are never much worse, and often better, than the best SUR strategy using a single level of fidelity.

Further work directions could be considered in the future. For instance, there is no explicit ingredient in MR-SUR strategies that tells the procedure to "learn" the model, and in particular, to learn the correlations between the levels of fidelity, and the variance of the output in the stochastic case. It seems to us that this would be important, particularly when simulations are very expensive and the simulation budget is very limited, as in our fire safety application. Using a fully Bayesian approach would somehow answer this problem,

as the uncertainty about the model would be propagated to the uncertainty about the QoI.

Another useful extension would be to address the case where the cost function is not known beforehand. As already mentioned, when sufficient prior knowledge about this function is available, a few pilot runs at several levels of fidelity are typically used to estimate the cost. In situations where such knowledge is lacking, the Bayesian approach can in principle be extended to this function as well (as suggested by [START_REF] Swersky | Multi-task Bayesian optimization[END_REF], which can then be learned on the fly. This idea, however, raises the question of how to adapt the MR-SUR strategy when the cost function is uncertain.

Finally, a third important research direction would be to address the design of parallel multi-fidelity simulations. As mentioned earlier in the article, leveraging batch-sequential SUR criteria makes it possible in principle to address the case of parallel batches of simulations with the same cost and run time-but the computational savings that can be obtained using such designs remain to be evaluated. Moreover, it seems even more promising to allow parallel simulations with different accuracies, costs and/or run times; but what would be a principled approach of (asynchronous) resource allocation in such a setting?

SM-3 A new SUR criterion SM-3.1 Criterion definition and result statement

Let ξ = (ξ(x)) x∈X denote a Gaussian process prior for the mean response of a stochastic simulator with Gaussian responses-i.e., a simulator which produces random responses

Z | ξ ∼ N (ξ(x), λ(x)) ,
where x ∈ X denotes the vector of inputs of the simulator and λ : X → [0, +∞) is a known variance function. Assume that the quantity of interest is the probability function

α : X → [0, 1] defined by α(x) = P Z x > z crit ,
where z crit ∈ R is a given threshold and Z x denotes a (future) response of the simulator with x as the input. More explicitely, we have

α(x) = Φ   ξ(x) -z crit λ(x)   . (SM1)
Let α n (x) denote the posterior mean of α(x) given n responses Z 1 , . . . , Z n of the simulator at design points X 1 , . . . , X n ∈ X (possibly selected in a sequential manner). Let µ denote a positive, bounded measure on X, and consider the measure of uncertainty H n defined by

H n = E n (α(x) -α(x)) 2 µ(dx) = var n (α(x)) µ(dx).
The following result provides a tractable expression for the corresponding SUR criterion for a batch x = ( xl ) 1≤l≤q ∈ X q of candidate points. The criterion presented in the article corresponds to the fully sequential case (q = 1) and to a particular choice of the measure µ.

Proposition 1. Let m n (resp. k n ) denote the posterior mean (resp. the posterior covariance) of ξ given n observations. Let

J n ( x) = E H n+q X n+1 = x1 , . . . , X n+q = xq and G n ( x) = H n -J n ( x). Then, J n ( x) = Φ 2 (u n (x), a n (x); r n (x)) -Φ 2 (a n (x), a n (x); r n (x; x)) µ(dx)
and

G n ( x) = Φ 2 (a n (x), a n (x); r n (x; x)) -Φ(a n (x)) 2 µ(dx),
with Φ the cdf of the standard normal distribution, Φ 2 (•, • ; ρ) the cdf of the standard bivariate normal distribution with correlation ρ, and

a n (x) = (m n (x) -z crit ) / v n (x), v n (x) = k n (x, x) + λ(x), r n (x) = k n (x, x)/v n (x), r n (x; x) = ν n (x, x; x)/v n (x), ν n (x, x ′ ; x) = k n (x, x ′ ) -k n+q (x, x ′ ; x) = k n ( x, x) T K n ( x, x) -1 k n ( x, x ′ ), k n ( x, x) = (k n ( xl , x)) 1≤l≤q , K n ( x, x) = (k n ( xl , x l ′ ) + λ( xl ) • δ l=l ′ )) 1≤l,l ′ ≤q .
Remark SM1. In the expressions of J n ( x) and G n ( x), the only part which depends on the future design x is Φ 2 (a n (x), a n (x); r n (x; x)) µ(dx), which must be maximized. Lemma 2. Let W ∼ N (m, K) be a d-dimensional Gaussian vector. Then, for any mean vector m ′ and covariance matrix K ′ ,

SM-3.2 A useful identity

E Φ d (W ; m ′ , K ′ ) = Φ d (m; m ′ , K + K ′ ).
Proof. Let W ′ ∼ N (m ′ , K ′ ) be independent from W . Then

E Φ d (W, m ′ , K ′ ) = E P W ′ ≤ W W = P(W ′ ≤ W ),
and, using that W ′ -W = W ′′m with W ′′ ∼ N (m ′ , K + K ′ ),

P(W ′ ≤ W ) = P(W ′′ ≤ m) = Φ d (m; m ′ , K + K ′ ). Corollary 3. Let m, m ′ ∈ R, v, v ′ ∈ [0, +∞) and W ∼ N (m, v). Then E Φ W -m ′ √ v ′ = Φ m -m ′ √ v + v ′ and E   Φ W -m ′ √ v ′ 2   = Φ 2 m -m ′ √ v + v ′ , m -m ′ √ v + v ′ ; v v + v ′ .

SM-3.3 Proof of Proposition 1

Proof. Recall from Equation (SM1) that Let us now compute separately the expectation with respect to P n of the two terms in the right-hand side of Equation (SM3). For the first term we have

E n E n+q α(x) 2 = E n α(x) 2 = E n   Φ   ξ(x) -z crit λ(x)
√ v + v ′ = m n (x) -z crit √ v n = a n (x), v v + v ′ = ν n (x, x; x) v n (x) = r n (x).
Combining Equations (SM3)-(SM5) and integrating on X with respect to µ yields the desired expression for J n (x). Similarly, combining Equation (SM2) with q = 0 and Equation (SM4) we have

H n =
Φ 2 (a n (x), a n (x); r n (x)) -Φ (a n (x)) 2 µ(dx)

and the expression of G n ( x) follows.

SM-4 Additional information regarding the examples SM-4.1 One-dimensional example

This section provides additional information regarding the Bayesian model used in the "One-dimensional example" presented in Section 4.2 of the article.

The model used in this example is the one proposed by [START_REF] Kennedy | Predicting the output from a complex computer code when fast approximations are available[END_REF] and reviewed in Section 2.1 of the article, with S = 2 levels. The two independent Gaussian processes η 1 and η 2 are stationary processes with unknown constant means and Matérn covariance functions with regularity 5/2:

η s ∼ GP m s , σ 2 s M 5/2 (a s (• -•)) , s ∈ {1, 2}.
Independent priors are used for all the remaining parameters of the model:

• Improper uniform prior distributions on R are used for the means m s .

• The parameters of the covariance function follow log-normal distributions: log(σ 2 s ) ∼ N (2 log(0.2), log(100) 2 ), log(a s ) ∼ N (log(2), log(10) 2 ).

• Finally, the regression term between the two levels follows a normal prior distribution:

ρ ∼ N (1, 2 2 ).

SM-4.2 Random damped harmonic oscillator

This section provides additional information regarding the explicit exponential Euler scheme used in the "Random damped harmonic oscillator" example (Section 4.3 of the article).

Consider a stochastic differential equation dY t = AY t dt + b dB t , where (Y t ) is a vector-valued stochastic process, (B t ) a Wiener process, b a vector and A a matrix. Let δ a time step, we would like to approximate (Y t ) by a finite difference method:

Y (δ) k
≈ Y kδ . The explicit exponential Euler scheme is a finite difference method proposed by [START_REF] Jentzen | Overcoming the order barrier in the numerical approximation of stochastic partial differential equations with additive space-time noise[END_REF] to ensure stability when approximating a stochastic equation.

The method is to apply recursively the formula

Y (δ) k+1 = exp (Aδ) Y (δ) k + √ 2πSδ • b • U
with S the spectral intensity of the Brownian motion, and U a normal random vector.

In particular, for the application of Section 4.2, we have

d   Y t Ẏt   =   0 1 -ω 2 0 -2ζω 0     Y t Ẏt   dt +   0 1   dB t .
Consequently, the approximation with a finite time-step δ is

  Y (δ) k+1 Ẏ (δ) k+1   = exp   0 δ -ω 2 0 δ -2ζω 0 δ       Y (δ) k Ẏ (δ) k   +   0 √ 2πSδu   
 , u ∼ N (0, 1).

  , Picheny et al. (2010), Chevalier et al. (2014), and Bect et al. (2019) for other examples of SUR criteria.

  entropy of the location of faulty component in an electronic equipment-with a discrete distribution over all possible fault locations as the underlying model. More recently, Snoek et al. (2012), Swersky et al. (2013), and Poloczek et al. (2017) have proposed Bayesian optimization procedures of the MR-SUR type, for unconstrained global optimization problems with variable-cost evaluations, corresponding respectively, when the cost is constant,

  U(R) denotes the (improper) uniform distribution on R and M ν the Matérn correlation function with regularity parameter ν. The values m = 0, σ 0 = 1, G = 4, L = 2, ν 0 = ν ε = 5/2, ρ 0 = 0.3, ρ ε = 0.1 are used in the simulations, and all the parameters except m are assumed to be known in this experiment. The cost function is defined as C(u, δ) = 1/δ. The QoI Q and the measure of residual uncertainty H n are defined as in Equations (

Figure 1 :

 1 Figure 1: Illustration of the MR-SUR criterion. (a) Representation of possible designs in the (C, J) plane. (b) Representation in the (C, G/C) plane. Each gray point corresponds to one particular input point x = (u, δ), the solid line indicates Pareto-optimal points, and the square is the design selected by the MR-SUR criterion.

Figure 2 :

 2 Figure 2: Sequence of Pareto fronts in (C, J)-space, as a function of the total cost of the design, in a particular run of the MR-SUR algorithm. On each curve, the square indicates the design selected by the MR-SUR criterion. A mixed black line connects the expected uncertainty at this point with the actual uncertainty obtained after the simulation has been run.

  MR-SUR: number of LF/HF eval.

Figure 3 :

 3 Figure 3: The one-dimensional experiment. (a) Median estimation error as a function of the cost. Lightgray disks: LF-SUR; dark-gray squares: HF-SUR; black diamonds: MR-SUR. (b) Number of LF/HF evaluations in the MR-SUR strategy. The last column indicates how many times, in the 60 repetitions of the experiment, a given combination appears (recall that HF evaluations are four times as costly as LF ones). For instance, 10 out of 60 repetitions, the sequential design was composed of 20 LF evaluations and 4 HF evaluations.

  1, where δ ref = 0.01s denotes the level of highest fidelity and z crit = -3 is a given critical threshold. The computational cost of Z is an affine function of 1/δ: C(δ) = a/δ + b. After normalization to have C(δ ref ) = 1, the coefficients are a = 9.792 ×10 -3 and b = 2.08 ×10 -2 .

Figure 4 :

 4 Figure 4: Median estimation error as a function of the cost for the oscillator test case.

Figure 5 :

 5 Figure 5: Three examples (out of 200) of the sequence of levels selected by the MR-SUR strategy for the oscillator test case.

  and Y soot = 0.027 kg • kg -1 . The environmental inputs u e are assumed random and integrated according to an environmental distribution P ue , which is a trivariate normal distribution with mean (10 • C, 100 kPa, 22.5 • C), variances equal to (20/3 • C, 2/3 kPa, 2.5 • C) 2 , and a correlation coefficient of 0.8 between the temperatures. The QoI is Q = Ue p(u e , u s , δ ref ) dP ue (u e ), where δ ref = 20 cm is the reference level and p(u e , u s , δ ref ) = P V < z crit |u e , u s , δ ref .

HFigure 6 :

 6 Figure 6: Result of four repetitions for the fire safety example. (a) Estimated probability as a function of the cost. The horizontal lines correspond to the Monte Carlo reference (dash-dotted line: mean; dotted lines: two-standard-deviation interval). (b) Square root of the measure of uncertainty (upper bound on the posterior standard deviation of the probability). The horizontal dotted line is the Monte Carlo standard deviation.

  Let Φ d ( • ; m, K) denote the cumulative distribution function of the d-variate normal distribution with mean m and covariance matrix K. The following identity is used byChevalier et al. (2014) for the computation of SUR criteria similar to ours.

  n+q (x)z crit λ(x) + k n+q (x, x) n+q (α(x)) = E n+q α(x) 2 -Φ   m n+q (x)z crit v n+q (x)

  (a n (x), a n (x); r n (x)) , (SM4) where we have applied the second part of Corollary 3 with m = m n (x), v = k n (x, x), m ′ = z crit , and v ′ = λ(x). For the second term we observe that, under P n , m n+q is a Gaussian process with mean m n and covariance function ν n ( •, • ; x). Therefore m n+q (x) ∼ N (m n (x), ν(x, x; x), and it follows that 2 (a n (x), a n (x); r n (x)) , (SM5)where we have used again the second part of Corollary 3 withm = m n (x), v = ν n (x, x; x), m ′ = z crit and v ′ = v n+q (x). Indeed, v + v ′ = ν n (x, x; x) + v n+q (x) = (k n (x, x)k n+q (x, x)) + (k n+q (x, x) + λ(x)) = v n (x), therefore mm ′

Table 1 :

 1 Levels of fidelity considered the damped harmonic oscillator example (Section 4.3). The highest level of fidelity is δ = 0.01 s. Note that the initial design only spans the first five levels.

	2 s 1/6 s 0.1 s 0.05 s 0.02 s 0.01 s

Table 2 :

 2 al. The levels of fidelity on FDS.

	Level δ	50 cm 33 cm 25 cm 20 cm
	Real cost	69 min	6 h	20 h	49 h
	Normalized cost 1/42	1/8	1/2.5	1
	Initial design	90	30	10	0

Acknowledgment. This work was supported by the French Association Nationale de la Recherche et de la Technologie (ANRT) under grant CIFRE 2014/1350.

SUPPLEMENTARY MATERIAL

Complementary information: Generalization and proof of the SUR sampling criterion (11), a short bibliography about nonsequential multi-fidelity designs, and more detailed information regarding the numerical experiments. (PDF file) (In the published version of the article, this supplementary material is provided as a separate PDF file, available from the web site of the publisher. In this authorgenerated post-print version: see below.)

Source code: Matlab programs to reproduce the numerical experiments of Sections 3.3, 4.2 and 4.3, including the modified version of the STK toolbox for Matlab/Octave [START_REF] Bect | STK: a Small (Matlab/Octave) Toolbox for Kriging[END_REF] that was used for these experiments.

SUPPLEMENTARY MATERIAL

SM-1 Introduction

This document contains supplementary material for the article "Sequential design of multifidelity computer experiments: maximizing the rate of stepwise uncertainty reduction". It is organized as follows. Section SM-2 provides a short literature review on nonsequential designs for multi-fidelity, which complements the literature review on sequential designs given in Section 3.1 of the article. Section SM-3 provides (with proof) a new SUR criterion, which covers as a special case the criterion provided (without proof) in Section 3.2 of the article. Finally, Section SM-4 provides additional information regarding the examples presented in Sections 4.2 and 4.3 of the article.

SM-2 Nonsequential designs in multi-fidelity

In this section, we provide a very brief literature review on nonsequential designs for multifidelity.

A common recommendation for multi-level designs is nesting. A multi-level design is nested when any observed point at a level δ (s) is also observed at every lower-fidelity level δ (s ′ ) , s ′ < s. Furthermore, space-filling designs are also expected to ensure observations in the whole input domain, as usual in Gaussian process regression.

A simple method to create a nested design is proposed by Forrester et al. (2007). It draws a maximin Latin Hypercube Sampling (LHS) at the lowest fidelity design, and then selects subsets of this LHS at the next levels. Le Gratiet and Garnier (2014) suggest to start with the highest-fidelity level and add points on the lower-fidelity levels to ensure better space-filling properties. The reader is also referred to Rennen et al. ( 2010) for a method which applies when there is only two levels of fidelity.

In our work, we use the method proposed by [START_REF] Qian | Nested latin hypercube designs[END_REF], which construct Nested Latin Hypercube Sampling (NLHS). An NLHS is a nested design with the property of being an LHS at each level of fidelity. We add a maximin optimization step to ensure that the design is space-filling at each level.

Note that this method was extended in several directions. He andQian (2011), Yang et al. (2014), [START_REF] Guo | Construction of latin hypercube designs with nested and sliced structures[END_REF]Xu et al. (2017) propose methods to build NLHS with particular structures, such as orthogonality.