Sequential design of experiments to estimate a probability of exceeding a threshold in a multi-fidelity stochastic simulator
Rémi Stroh Séverine Demeyer Nicolas Fischer Julien Bect Emmanuel Vazquez
stat.CO stat.ME stat.ML





  
  
  Sequential design of experiments to estimate a probability
    of exceeding a threshold in a multi-fidelity stochastic simulator





  Rémi STROH*

  lne, Trappes, France - remi.stroh@lne.fr

  l2s, CentraleSupélec, Univ. Paris-Sud,
    cnrs, Université Paris-Saclay, Gif-sur-Yvette, France




Séverine DEMEYER, Nicolas FISCHER

lne, Trappes, France - name.surname@lne.fr




Julien BECT, Emmanuel VAZQUEZ

l2s, CentraleSupélec, Univ. Paris-Sud, cnrs,
  Université Paris-Saclay, Gif-sur-Yvette,
  France - name.surname@centralesupelec.fr






  In this article, we consider a stochastic numerical simulator to
  assess the impact of some factors on a phenomenon. The simulator is
  seen as a black box with inputs and outputs.  The quality of a
  simulation, hereafter referred to as fidelity, is assumed to be
  tunable by means of an additional input of the simulator (e.g., a mesh
  size parameter): high-fidelity simulations provide more accurate
  results, but are time-consuming. Using a limited computation-time
  budget, we want to estimate, for any value of the physical inputs, the
  probability that a certain scalar output of the simulator will exceed
  a given critical threshold at the highest fidelity level.
  The problem is addressed in a Bayesian framework, using a Gaussian
  process model of the multi-fidelity simulator. We consider a Bayesian
  estimator of the probability, together with an associated measure of
  uncertainty, and propose a new multi-fidelity sequential design
  strategy, called msur, to select the value of physical
  inputs and the fidelity level of new simulations. The msur
  strategy is tested on an example.


  Keywords: Multi-fidelity;
  Sequential design; Bayesian analysis; Gaussian process.





§ INTRODUCTION


The objective of this article is to propose a new Bayesian algorithm
for sequential design of experiments in the context of
multi-fidelity stochastic numerical simulators. 
A numerical simulator is a computer program modeling a physical
phenomenon or a system. When the simulator is deterministic, running
the computer program twice using the same inputs yields the same
outputs. 
In this case, the simulator can be viewed formally as a function. When
the simulator is stochastic, running twice the simulator with the same
inputs does not return the same output. 
Moreover, we assume that the simulator has a particular input called
fidelity parameter—e.g., the mesh size of a finite-difference
partial differential equation solver—that controls a trade-off
between quality of simulation and computation time. 
A high-fidelity simulation provides an accurate result, but is
time-consuming.

Let t ∈T be the fidelity input of the simulator,
and x ∈X the vector of all other inputs. 
At fixed (x,  t), the output of the simulator Z follows
a probability distribution ^Z_x,  t. 
We denote by  the value of the fidelity input associated to the
highest available level of fidelity: 
this is the level of interest for the user of the simulator. 
In this article, we focus on comparing the output Z to a critical
threshold —typically, a level of output that should not be
exceeded by the physical phenomenon or the system under study. 
This comparison can be studied from two points of view: either
locally, by considering the function  that gives for each
point x the probability of exceeding the threshold at this point,

    (x) = ^Z_x, (Z >
        ),   x∈X,

or globally, by computing the global probability  of exceeding
the threshold, which may be written as

    = ∫_X(x) (),

where  is a probability distribution which models uncertainty
about the value of the input factors.

Our objective is to estimate these quantities by observing the
stochastic outcomes z_1, z_2… of simulations at points
(x_1,  t_1),  (x_2, t_2)… 
The estimators are built using a Bayesian model of the simulator. 
More precisely, we assume a prior model about ^Z_x, t and
compute the posterior distribution of the quantity of interest (either
 or ) given
χ_n = (x_i,  t_i;  z_i)_i≤ n. 
In this article, we suggest a new sequential design algorithm to
select points (x_1,  t_1),  (x_2, t_2)… in order to obtain a
fast reduction of the uncertainty about  or . 
Our methods deals both with the stochastic nature of the simulator and
the tunable simulation cost.

The paper is organized as follows. 
Section <ref> sets the Bayesian framework that we consider for
stochastic multi-fidelity simulators. 
Section <ref> describes our sequential algorithm and our
new sampling criterion called msur. 
Section <ref> illustrates the method on a test
problem, and Section <ref> concludes the paper.




§ BAYESIAN FRAMEWORK



Following our previous work on multi-fidelity stochastic simulators
<cit.>, assume that the output Z
of a simulation at (x,  t) follows a normal distribution:

    Z |ξ, λ ∼ (ξ (x,  t), λ(x,  t)) .

Assume moreover that conditional on ξ and λ, all runs of
the simulator are independent. 
Of course, the choice of a normal distribution here is a simplifying
hypothesis, that must be verified in practice for each particular
simulator
[
  Other (possibly non-parametric) families of distributions could be
  considered as well. 
  Note, however, that the convenient conjugacy property of the
  Gaussian process prior with respect to the Gaussian likelihood would
  be lost by doing so.].

Assume a Gaussian process prior for the mean function ξ:

    ξ| m, k
       ∼ (m, 
        k((x,  t), (x',  t'))),
      
           m  ∼ 𝖴(R),

where m is the (constant) mean of the Gaussian process, k its
covariance function, and 𝖴(R) the
(improper) uniform distribution on R. 
For the sake of simplicity, the functions λ and k are
assumed to be known. 
Under this prior, the posterior distribution of ξ conditionally
to the observations χ_n is Gaussian:

    ξ|χ_n  ∼ (m_n(x,  t), 
        k_n((x,  t), (x',  t'))),

where m_n and k_n are given by the kriging formulas <cit.>.

Since Z has a Gaussian distribution conditional on ξ and
λ, for all x∈X, the probability of exceeding a
critical value  can be written as

    (x) =
      Φ(ξ(x, ) - /√(λ(x, ))),

with Φ the cumulative distribution function of the normal
distribution. 
Then, denoting by _n = (  ·|χ_n)
and _n = (   ·|χ_n) the posterior
mean and variance operators, we have: 

    _n((x))
          = Φ(u_n(x) ),
    _n((x))
          = Φ_2(u_n(x), u_n(x);
        r_n(x, x))
        - Φ(u_n(x) )^2,

where Φ_2( ·, · ;  r) stands for the
bivariate normal distribution function with correlation r,

V_n(x,  t) = λ(x,  t) + k_n((x,  t),
    (x,  t)),

u_n(x) = m_n(x, ) - /√(V_n(x, )), and

r_n(x,  x') = k_n((x, ),(x', ))/√(V_n(x, ) V_n(x', ) ).




§ PROPOSED ALGORITHM





 §.§ A SUR criterion to estimate a probability of
  exceeding a threshold in the case of stochastic outputs


The goal of sequential design of experiments is to select a simulation
point (x_n+1, t_n+1), using the result of the previous
observations χ_n, to obtain a good improvement of our knowledge
about a given quantity of interest.

Our algorithm is based on the idea of sur strategies
<cit.>. 
The main idea of SUR strategies is to define a measure of
uncertainty _n about the estimator of a quantity of
interest.

Then, assuming that there is no fidelity parameter t,
the next observation point is selected in order to minimize the
expected uncertainty using this new observation, whose outcome is
random:

    x_n+1 = _x ∈X _n(_n+1|
        X_n+1 = x).


Several measures of uncertainty _n have been proposed in the
literature <cit.> for the problem of estimating the volume or
contour of an excursion set of a deterministic function. 
In our framework, however, we deal with a stochastic simulator, and we
propose a new measure of uncertainty for this case.  Let L be the
L^2 loss function

    (f, g)
      = f - g^2_
      = ∫_X(f(x) - g(x))^2 μ(),

where μ denotes a positive measure on X. 
We suggest to measure uncertainty by the L^2
loss incurred when estimating  by _n = _n():

    _n
      = _n(
      (_n, ))
      = ∫_X_n( (x) ) μ().

Note that in the case where the goal is to estimate ,
_n provides an upper bound on the posterior variance for
any  that admits a density g with respect to μ: 
_n() ≤ G _n, where
G = ∫ g^2 and g =  /.

Consequently, a sequential design strategy that aims to
reduce _n will be useful not only to estimate  itself
but also to estimate  for a large class of probability
distributions .

Let J_n(x,  t) = _n(_n+1|
  X_n+1 = x, T_n+1 = t). Since the output of the
simulator has a Gaussian distribution, J_n(x,
  t) can be written as

    J_n(x,  t) = ∫_X(
      Φ_2(u_n(y), u_n(y);
      r_n(y, y))
      - Φ_2(u_n(y), u_n(y);
      k_n((y, ),
          (x,  t))^2/V_n(y, )V_n(x,  t))
      ) μ ().





 §.§ Dealing with tunable fidelity


The sur strategy (<ref>) is relevant when there is
no fidelity parameter t. 
In a multi-fidelity context, however, the simulator becomes more and
more time-consuming as the fidelity parameters comes closer
to . 
We propose a new sequential design algorithm which takes into account
the variable simulation cost.

Let C(x,  t) > 0 be the cost of observing the simulator at
(x,  t).  The idea is to balance the benefit of an
observation against its cost.  For us, the benefit of an observation is
the reduction of uncertainty.  Hence, we propose the following
sequential strategy called msur:

    (x_n+1; t_n+1) =
      _(x, t) ∈X×T _n - J_n(x, t)/C(x, t).

This idea is similar to the one proposed by
<cit.> in the context of optimization, where the
expected improvement is divided by the cost of an observation.
Likewise, for the purpose of providing a surrogate model of the simulator,
<cit.> suggested to choose the new level
of fidelity by comparing the cost of the level and the
reward in terms of reduction of uncertainty.

Note that in many applications, the cost C depends only on the level
t: C(x,  t) = C(t).  In this case,
(<ref>) can be divided into two steps: first, select the
input x at each level that minimizes the sur criterion; second,
select the level t that maximizes the msur criterion.  Our
sequential design algorithm can thus be rewritten in this case as

    {x̃(t)     =
           _x ∈X  J_n(x, t),
    
           t_n+1    = _t ∈T _n
             - J_n(x̃(t), t)/C(t),
    
           x_n+1    = x̃(t_n+1).
    .




§ ILLUSTRATION



We illustrate our algorithm on a two-dimensional example
inspired by <cit.>.  We consider a damped harmonic
oscillator subject to  random forcing, described by the stochastic
differential equation

    Ẍ(u) + 2ζω_0Ẋ(u)
      + ω_0^2X(u) = W(u),
         X(0) = 0,
        Ẋ(0) = 0,   u ∈[0; = 30],

where W is a Gaussian white noise (with spectral density equal to
one).  We compute an approximation (X_n)_n=0,…,
  ⌊/⌋ of X by finite differences with a time
step .  We use an explicit exponential Euler scheme <cit.>. The output of our simulator is defined as

    [                           Z:      [0; 30] ×[0; 1] ×]0; 1]                            →                            R;                                                (ω_0, ζ, )                            ↦ log(max{X_n, n
          ≤⌊/⌋}). ]


The pair (ω_0, ζ) corresponds to the input vector
x and the step time  plays the role of the fidelity
parameter. The critical threshold  is set to -3. 
The highest level of fidelity is set to = 0.01. 
The computational cost is empirically linear with respect to the
fidelity level: C() = a/ + b, with a = 0.0098 and
b = 0.02 (coefficients chosen to have C( ) = 1). 
Figure <ref> represents the contour plots of the mean function,
the standard deviation and the probability function at three different
levels of fidelity.



We compute a reference value  of  on a 100× 100
regular grid at the highest level of fidelity using 10^4 simulations
of Z (see Figure <ref>).  We suppose that the
input distribution  is the uniform distribution on [0;
  30]×[0; 1] and let μ =.  The mean value =
83.3% of  on the grid serves as a reference value for
.

We use the Bayesian model of Section <ref> with fixed
hyper-parameters, and ten fixed levels of fidelity: = 1, 0.5,
0.33, 0.25, 0.2, 0.17, 0.1, 0.05, 0.02 and 0.01.  The
initial design (initial simulations before applying our msur strategy)
consists of 180× 60× 20 × 10 × 5 nested observations on
the five lowest fidelity levels (180 on the level =
1, 60 on the level = 0.5, …). The
initial design is set using the algorithm of <cit.>.

In our experiments, we compare two strategies: five sl
strategies and the msur strategy (<ref>). 
For sl strategies, all new points are sequentially selected on a
fixed level, using then sur strategy (<ref>). 
Each sl strategy corresponds to one level of fidelity. 
Each time a new observation point must be selected, we choose the
point that achieves the best value of J_n among 500 candidate points
per level drawn from . 
We allocate a simulation-time budget of 20 for each strategy. 
All experiments are repeated 12 times. 
Integrals are approximated by a sum on the 100× 100 regular
grid.

The strategies are compared based on the mean square error between
estimations and references. 
Results are presented on Figure <ref>.
Figure <ref> represents the mean square error of the
estimation _n of , and
Figure <ref> represents the mean square error of
the estimation _n of , as a function of the
computational cost.  Each curve corresponds to an average over 12
experiments with the same method of sequential design.  The plain red
curve corresponds to the msur strategy.  The dotted crossed
blue-green curves correspond to sl strategies (one curve per
level). We can see that the best level is achieved at =
0.05.  The lower fidelity levels are too biased to
estimate the probabilities correctly, and the upper fidelity levels are
too expensive to make it possible to carry out enough observations in
order to estimate the probabilities accurately.  Moreover, we can see
that the msur strategy is as good as the sl strategy at =
0.05.  Consequently, with the msur strategy, one does
not need to know which level yields the best trade-off between  accuracy
and computational cost.






§ CONCLUSION
 

This article makes two contributions. 
First, we suggest a new sur criterion to estimate a probability
of exceeding a threshold in the case of a stochastic simulator. 
Second, we construct a sequential strategy called msur as an
adaptation of the sur strategies to deal with multi-fidelity
simulators. 
Our first results are promising, because the msur strategy
succeeds to get the better of all single-level sur strategies
without knowing which level provides the best compromise between speed
and accuracy.


apalike


