Gaussian random variables ◮ The mean of Z (also called expectation or first-order moment) is

and its second-order moment is defined as

Gaussian random variables ◮ A Gaussian variable Z can be used as a stochastic model of some uncertain real-valued quantity ◮ In other words, Z can be thought as a prior about some uncertain quantity of interest 31/224 Gaussian random variables ◮ Using a random generator, it is possible to "generate" sample values z 1 , z 2 , . . . of our model Z → possible values for our uncertain quantity of interest -6 -4 -2 0 2 4 6 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 n=100 32/224

Bayesian model

Formally, recall that a statistical model is a triplet M = (Z, F, P)

Def.

A Bayesian model is defined by the specification of ◮ a parametric statistical model M (model of the observations)

◮ a prior probability distribution Π on (Θ, Ξ) → probability model that describes uncertainty about θ before an observation is made

72/224

Generalization Exercise: Let (Z 0 , Z 1 , . . . , Z n ) be a centered Gaussian vector. Determine E (Z 0 | Z 1 , . . . , Z n ).

73/224

Application: computation of the posterior distrib. of a GP

what is the conditional distrib.-or the posterior distrib. in our Bayesian framework-of ξ(x 0 ) | ξ(x 1), . . . , ξ(x n ) ? ◮ More generally, what is the posterior distribution of the random process ξ(•) | ξ(x 1), . . . , ξ(x n ) ? 74/224

Computation of the posterior distribution of a GP Prop.

Let ξ ∼ GP (0, k). The random process ξ conditioned on

Computation of the posterior distrib. of a GP

◮ By property of the conditional expectation in Gaussian spaces, for all x ∈ X, ξ n (x ) is a linear combination of ξ(x 1 ), . . . , ξ(x n ):

◮ Moreover, the posterior mean ξ n (x ) is the orthogonal projection of ξ(x ) onto span{ξ(x i ), i = 1, . . . , n}, such that

is the best linear predictor (BLP) of ξ(x ) from ξ(x 1 ), . . . , ξ(x n ), also called the kriging predictor of ξ(x )

◮ Exactly computable for independent GP priors, 2 ≤ p 5 ◮ Implemented in STK (Matlab/Octave), GPareto (R). . . ◮ Dependent priors, larger p: Monte Carlo approx.

◮ Many other sampling criteria have been proposed

◮ Suppose that we are given a data set of n simulation results, i.e., evaluations results of an unknown function f : [0, 1] → R, at points x 1 , . . . , x n .

◮ A data set of size n = 8: A simple curve fitting problem ◮ Any approximation procedure of f consists in building a function f n = h(•; θ) where θ ∈ R l is a vector of parameters, to be estimated from D n and available prior information

◮ Fundamental example: linear model

f n (x ) = h(x , θ) = l i=1 θ i r i (x )
where functions r i : X → R are called regressors (e.g., r 1 (x ) = 1, r 2 (x ) = x , r 3 (x ) = x 2 . . . → polynomial model)

◮ Most classical method to obtain a good value of b: least squares → minimize the sum of squared errors

J(θ) = n i=1 z i -f n (x i ) 2 = n i=1 (z i -h(x i ; θ)) 2 14/224
A simple curve fitting problem ◮ Linear fit A simple curve fitting problem ◮ When the problem becomes ill-defined (as capacity increases), a classical solution for finding a good value of b is to minimize the sum of an approximation error and a regularization term:

J(θ) = n i=1 z i -θ t r (x i ) 2 + C θ 2 2 , C > 0
◮ θ 2 2 penalizes vectors θ with large elements ◮ C strikes a balance between regularization and data fidelity ◮ This approach is known as Tikhonov regularization (Tikhonov & Arsenin, 1977) → at the basis of numerous approximation methods (ridge regression, splines, RBF, SVM. . . ) Example ◮ Suppose we repeat measurements of a quantity of interest:

z 1 , z 2 , . . . ∈ R ◮ Model of observations: Z i iid ∼ N (θ 1 , θ 2 ), i = 1, . . . , n
◮ The statistical model can formally be written as the triplet

M = R n , B(R n ), N (θ 1 , θ 2 ) ⊗n θ 1 ,θ 2
◮ Moreover, if we a assume a prior distribution about θ 1 and θ 2 (e.g. θ 1 ∼ N (1, 1), θ 2 ∼ IG(3, 2)), we obtain a Bayesian model

◮ From this Bayesian model (model of observations + prior), we can compute the posterior distribution of (θ 1 , θ 2 ) given Z 1 , . . . , Z n (will be explained later) Simple curve fitting problem from a Bayesian approach ◮ Recall our simple curve fitting model

f n (x ) = θ t r (x )
with r (x ) = (1 cos(2πx ) sin(2πx ) . . . cos(2mπx ) sin(2mπx )) t ∈ R 2m+1 ◮ Bayesian model?

36/224

◮ Assume the following statistical model for the observations:

       Z i = ξ(x i ) + ε i , i = 1, . . . , n ξ(x ) = θ t r (x ), x ∈ X ε i iid ∼ N (0, σ 2 ε )
or equivalently,

Z i iid ∼ N (θ t r (x i ), σ 2 ε ), i = 1, . . . , n
◮ Moreover, choose a prior distribution for θ:

θ j indep ∼ N (0, σ 2 θ j ), j = 1, . . . ,

Random process

Def.

A random process ξ : (Ω, X) → R is a collection of random variables ξ(•, x ) : Ω → R indexed by x ∈ X ◮ Random processes can be viewed as a generalization of random vectors

◮ For a fixed ω ∈ Ω, the function ξ(ω, •) : X → R is called a sample path

38/224

◮ In our Bayesian setting, we can say that ◮ we use a random process ξ as a stochastic model of the unknown function f

◮ f is viewed as as sample paths of ξ ◮ ξ represents our knowledge about f before any evaluation has been made

◮ the distribution Π = P ξ is a prior about f

All real functions

Prior P ξ Bayesian approach

◮ The choice of a prior in a Bayesian approach reflects the user's knowledge about uncertain parameters ◮ In the case of function approximation → regularity of the function

42/224

Bayesian approach

◮ Where shall we go now?

◮ Objective: compute posterior distributions from data A "simplification" of the Bayesian model for the simple curve fitting problem ◮ ξ = θ t r (our prior about f ) is a Gaussian process ◮ Why?

44/224

Gaussian random vectors Def.

A real-valued random vector

Z = (Z 1 , . . . , Z d ) ∈ R d is said to be Gaussian iff any linear combination of its components d i=1 a i Z i , with a 1 , . . . , a d ∈ R, is a Gaussian variable ◮ A Gaussian random vector Z is characterized by its mean vector, µ = (E[Z 1 ], . . . , E[Z d ]) ∈ R d , and the covariance of the pairs of components (Z i , Z j ), i, j ∈ {1, . . . , d}, cov(Z i , Z j ) = E[(Z i -E(Z i ))(Z j -E(Z j ))] ◮ If Z ∈ R d is a Gaussian vector with mean µ ∈ R d and covariance matrix Σ ∈ R d×d , we shall write Z ∼ N (µ, Σ) 45/224 ◮ Exercise: Let Z ∼ N (µ, Σ). Determine E d i=1 a i Z i and var d i=1 a i Z i ◮ The correlation coefficient of two components Z i and Z j of Z is defined by ρ(Z i , Z j ) = cov(Z i , Z j ) var(Z i )var(Z j ) ∈ [-1, 1],
→ measures the similarity between Z i and Z j

46/224

Gaussian random vectors: correlation

ρ = 0 -4 -3 -2 -1 0 1 2 3 4 -3 -2 -1 0 1 2 3 ρ = 0.8 -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3

47/224

Gaussian random processes ◮ Recall that a random process is a set ξ = {ξ(x ), x ∈ X} of random variables indexed by the elements of X ◮ Gaussian random process → generalization of a Gaussian random vector

Def.

ξ is a Gaussian random process iff, ∀n ∈ N, ∀x 1 , . . . , x n ∈ X, and ∀a 1 , . . . , a n ∈ R, the real-valued random variable

n i=1 a i ξ(x i ) is Gaussian 48/224 Application ◮ If    ξ(x ) = θ t r (x ), x ∈ X θ j indep ∼ N (0, σ 2 θ j ), j = 1, . . . , 2m + 1 then, ∀x 1 , . . . , x n ∈ X, and ∀a 1 , . . . , a n ∈ R, n i=1 a i ξ(x i ) = i a i j θ j r j (x i ) = j i a i r j (x i ) θ j ∼ N 0, j i a i r j (x i ) 2 σ 2 θ j ◮ Thus, ξ = θ t r is a Gaussian process 49/224
Gaussian random processes ◮ A Gaussian process is characterized by

◮ its mean function m : x ∈ X → E[ξ(x )]
◮ and its covariance function

k : (x , y ) ∈ X 2 → cov(ξ(x ), ξ(y )) ◮ Notation: ξ ∼ GP (m, k) 50/224 ◮ Exercise: determine E d i=1 a i ξ(x i ) and var d i=1 a i ξ(x i ) ◮ What is the distribution of a i ξ(x i )?
→ The distribution of a linear combination of a Gaussian process GP(m, k) can be simply obtained as a function of m and k

51/224

Application

◮ If    ξ(x ) = θ t r (x ), x ∈ X θ j indep ∼ N (0, σ 2 θ j ), j = 1, . . . , 2m + 1 then, ξ ∼ GP(0, k) with k : (x , y ) → j σ 2 θ j r j (x )r j (y ) 52/224 ◮ Covariance function corresponding to      θ 1 ∼ N (0, 1) θ 2k , θ 2k+1 indep ∼ N 0, 1 1+(ω 0 k) α , k = 1, . . . , m with ω 0 = 2π 10 , α = 4 -0.5 0 0.5 h -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 k 52/224 ◮ Covariance function corresponding to      θ 1 ∼ N (0, 1) θ 2k , θ 2k+1 indep ∼ N 0, 1 1+(ω 0 k) α , k = 1, . . . , m with ω 0 = 2π 50 , α = 4 -0.5 0 0.5 h -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 3 k 53/224
The covariance of a Gaussian random process ◮ Main properties: a covariance function k is

◮ symmetric: ∀x , y ∈ X, k(x , y ) = k(y , x ) ◮ positive: ∀n ∈ N, ∀x 1 , . . . , x n ∈ X, ∀a 1 , . . . , a n ∈ R, n i,j=1 a i k(x i , x j )a j ≥ 0 ◮
In the following, we shall assume that the covariance of ξ ∼ GP(m, k) is invariant under translations, or stationary:

k(x + h, y + h) = k(x , y ), ∀x , y , h ∈ X ◮ When k is stationary, there exists a stationary covariance k sta : R d → R such that cov(ξ(x ), ξ(y )) = k(x , y ) = k sta (x -y ) 54/224
Stationary covariances

◮ When k is stationary, the variance var(ξ(x )) = cov(ξ(x ), ξ(x )) = k(0)
does not depend on x ◮ The covariance function can be written as

k(x -y ) = σ 2 ρ(x -y ) ,
with σ 2 = var(ξ(x )), and where ρ is the correlation function of ξ.

55/224

Stationary covariances

◮ The graph of the correlation function is a symmetric "bell curve" shape

-4 -3 -2 -1 0 1 2 3 4 0 0.2 0.4 0.6 0.8 1 correlation range 56/224 ◮ We have, by C-S, ∀h ∈ X, |k(h)| = cov(ξ(x ), ξ(x + h)) = E[(ξ(x ) -m(x ))(ξ(x + h) -m(x + h))] ≤ E((ξ(x ) -m(x )) 2 ) 1/2 E((ξ(x + h) -m(x + h)) 2 ) 1/2 = k(0) 1/2 k(0) 1/2 = k(0) ◮ Recall, Bochner's spectral representation theorem Theorem A real function k(h), h ∈ R d is symmetric positive iff it is the Fourier transform of a finite positive measure, i.e. k(h) = R d e ı(u,h) dµ(u) ,
where µ is a finite positive measure on R d .

57/224

Gaussian process simulation ◮ Using a random generator, it is possible to "generate" sample paths f 1 , f 2 , . . . of a Gaussian process ξ

-2 -1 0 1 2 0 0.2 0.4 0.6 0.8 1 -2 -1 0 1 2 -4 -2 0 2 4 -2 -1 0 1 2 0 0.2 0.4 0.6 0.8 1 -2 -1 0 1 2 -4 -2 0 2

58/224

Gaussian process simulation

How to simulate sample paths of a zero-mean Gaussian random process?

◮ Choose a set of points x 1 , . . .

x n ∈ X ◮ Denote by K the n × n covariance matrix of the random vector ξ = (ξ(x 1 ), . . . , ξ(x n )) t (NB: ξ ∼ N (0, K )) ◮ Consider the Cholesky factorization of K K = CC t ,
with C a lower triangular matrix (such a factorization exists since K is a sdp matrix)

◮ Let ε = (ε 1 , . . . , ε n ) t a Gaussian vector with ε i i.i.d ∼ N (0, 1) ◮ Then C ε ∼ N (0, K ) 59/224
"Simplification" of the Bayesian model for the simple curve fitting problem

◮ Instead of choosing the model    ξ(x ) = θ t r (x ), x ∈ X θ j indep ∼ N (0, σ 2 θ j ), j = 1, . . . , 2m + 1
simply choose a covariance function k and assume ξ ∼ GP(0, k)

◮ More details about how to choose k will be given below

60/224

Posterior ◮ How to compute a posterior distribution, given the Gaussian-process prior ξ and data?

0 0.2 0.4 0.6 0.8 1 x -4 -3 -2 -1 0 1 2 3 4 z 61/224

Conditional distributions

◮ Let X be a random variable modeling an unknown quantity of interest ◮ Assume we observe a random variable T , or a random vector T = (T 1 , . . . , T n )

◮ Provided that T and X are not independent, T contains information about X ◮ Aim: define a notion of distribution of X "knowing" T

62/224

Conditional probabilities

Recall the following Def.

Let (Ω, A, P) be a probability space. Given two events A, B ∈ A such that P(B) = 0, define the probability of A given B (or conditional on B) by

P(A | B) = P(A ∩ B) P(B)

63/224

The notion of conditional density Def.

Assume that the pair (X , T ) ∈ R 2 has a density p (X ,T ) .

Define the conditional density of X given the event T = t by

p X |T (x | t) =        p (X ,T ) (x , t) p T (t) = p (X ,T ) (x , t) p X ,T (x , t)dx if p T (t) > 0 arbitrary density if p T (t) = 0.

64/224

Application ◮ Given a Gaussian random vector

Z = (Z 1 , Z 2 ) ∼ N (0, Σ), what is the distribution of Z 2 "knowing" Z 1 ? ◮ Define p Z 2 |Z 1 (z 2 |z 1 ) = p (Z 2 ,Z 1 ) (z 2 , z 1 ) p Z 1 (z 1 ) = σ 1 (2π) 1/2 (det Σ) 1/2 exp(-1/2(z t Σ -1 z-z 2 1 /σ 2 1 )) → Gaussian distribution! -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 ◮

Def.

Assume that the pair (X , T ) has a density p (X ,T ) . Define the conditional mean of X given T = t as

E(X | T = t) ∆ = R x p X |T (x |t)dx = h(t)
Def.

The random variable E(X | T ) = h(T ) is called the conditional expectation of X given T

66/224

Conditional expectation ◮ Why conditional expectation is an fundamental notion?

◮ We have the following

Theorem

Under the previous assumptions, the solution of the problem

X = argmin Y E (X -Y ) 2
where the minimum is taken over all functions of T is given by X = E(X | T )

◮ In other words, E(X | T ) is the best approximation (in the sense of the quadratic mean) of X by a function of T

67/224

Important properties of conditional expectation

(1) E(X | T ) is a random variable depending on T → there exists a function h such that E(X | T ) = h(T )

(2) The operator π H : X → E(X | T ) is a (linear) operator of orthogonal projection onto the space of all functions of T (for the inner product X

, Y → (X , Y ) = E (XY )) (3) Let X , Y , T ∈ L 2 . Then i) ∀α ∈ R, E(αX + Y | T ) = αE(X | T ) + E(Y | T ) a.s. ii) E(E[X | T ]) = E(X ) iii) If X ⊥ ⊥ T , E(X | T ) = E(X )

68/224

Conditional expectation: Gaussian random variables ◮ Recall that the space of second-order random variables L 2 (Ω, A, P) endowed with the inner product X , Y → (X , Y ) = E (XY ) is a Hilbert space ◮ Gaussian linear space

Def.

A linear subspace G of L 2 (Ω, A, P) is Gaussian iff Theorem (Projection theorem in centered Gaussian spaces)

∀X 1 , . . . , X n ∈ G
Let G be a centered Gaussian space. Let X , T 1 , . . . ,

T n ∈ G. Then E(X | T 1 , . . . , T n ) is the orthogonal projection (in L 2 ) of X on T = span{T 1 , . . . , T n }.
proof Let X ∈ G be the orthogonal projection of X on T .

◮ we have X = X + ε where ε ∈ G is orthogonal to T . In G, orthogonality ⇔ independence. Thus, ε ⊥ ⊥ T i , i = 1, . . . , n. ◮ Then, E(X | T 1 , . . . , Tn) = E( X | T 1 , . . . , Tn) + E(ε | T 1 , . . . , Tn) = X + E(ε) = X
The result follows.

70/224

Application Let Z = (Z 1 , Z 2 ) be a zero-mean Gaussian random vector, with covariance matrix

σ 2 1 σ 1,2 σ 2,1 σ 2 2 Then E(Z 1 | Z 2 ) is the orthogonal projection of Z 1 onto Z 2 . Thus E(Z 1 | Z 2 ) = λZ 2 with (Z 1 -λZ 2 , Z 2 ) = (Z 1 , Z 2 ) -λ(Z 2 , Z 2 ) = 0 . Hence, λ = σ 1,2 σ 2 71/224 Application ◮ Let Z = (Z 1 , Z 2 ) ∼ N (0, Σ) as above → recall that the cond. distrib. of Z 2 given Z 1 is a Gaussian distribution ◮ Hence Z 2 | Z 1 ∼ N (µ(Z 1 ), σ(Z 1 ) 2 ) → µ ? σ ? -4 -3 -2 -1 0 1 2 3 4 -4 -3 -2 -1 0 1 2 3 4 ◮ µ = E(Z 2 | Z 1 ) = λZ 1 with λ = σ1,2 σ 2 2
◮ Using the property of the orthogonal projection:

σ 2 = E (Z 2 -E(Z 2 | Z 1 )) 2 | Z 1 ⊥ = E (Z 2 -E (Z 2 | Z 1 )) 2 = E (Z 2 -λZ 1 ) 2 ⊥ = E ((Z 2 -λZ 1 )Z 2 ) = σ 2 2 -λσ 1,2 76/224
Computation of the posterior distrib. of a GP ◮ The posterior covariance, also called kriging covariance, is given by

k n (x , y ) := cov ξ(x ) -ξ n (x ), ξ(y ) -ξ n (y ) = k(x -y ) - i λ i (x ) k(y -x i ) .
◮ k n is the covariance function of the error of prediction ◮ The posterior variance of ξ, also called the kriging variance, is defined as

σ 2 n (x ) = var(ξ(x ) -ξ n (x )) = k n (x , x ) ◮ σ 2 n (x )
is the variance of the error of prediction

77/224

Kriging equations ◮ How to compute the weights λ i (x ) of the posterior mean/kriging predictor?

◮ Weights λ i (x ) are solutions of a system of linear equations Posterior distrib. of a GP ◮ For all x ∈ X, the random variable ξ(x ) | F n with distrib. N ξ n (x ), σ 2 n (x ) represents the residual uncertainty about ξ(x ) when ξ(x 1 ), . . . , ξ(x n ) are observed Generalization: prediction from noisy observations

K λ(x ) = k(x ) with -λ(x ) = (λ 1 (x ), . . . , λ n (x )) t -K : n × n covariance
-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 ξ n (x ) +2σ n (x )
◮ Let ξ ∼ GP(0, k) ◮ For i = 1, . . . , n, we observe Z i = ξ(x i ) + ε i at points x i ,
where the random variables ε i model an observation noise:

◮ ε i i.i.d ∼ N (0, σ 2 )
◮ independent from ξ ◮ As above, the posterior mean ξ n (x ) of ξ(x ) is obtained as the orthogonal projection of ξ(x ) on the linear subspace span{Z i , i = 1, . . . , n}:

ξ n (x ) = n i=1 λ i (x )Z i with λ i (x ) such that ∀i, ξ(x ) -ξ n (x ) ⊥ Z i 83/224
Generalization: prediction from noisy observations ◮ Thus, ∀i

E (ξ(x ) -ξ n (x ))Z i = E [ξ(x )(ξ(x i ) + ε i )] - n j=1 λ j (x )E [(ξ(x j ) + ε j )(ξ(x i ) + ε i )] = k(x , x i ) - n j=1 λ j (x ) k(x j , x i ) + σ 2 δ i,j
◮ Under matrix form (exercise):

84/224

Generalization: prediction from noisy observations

Generalization: prediction from noisy observations

We should use an approximation instead of an interpolation in three cases:

i) The observations are noisy (obviously): the computer code is stochastic (for instance, Monte Carlo is used) and running the code twice does not produce the same output

ii) The output of the computer code is very irregular → a smooth approximation is preferred

iii) The covariance matrix is ill-conditioned → adding a small observation noise will regularize the solution of the linear system (why?)

86/224

Lecture 1 : From meta-models to UQ ◮ What are plausible values for the minimum M = min x ξ(x )?

◮ Where is the minimizer x ⋆ = argmin x ξ(x )?

◮ What is the probability that ξ(x ) exceeds a given threshold? ◮ . . .

88/224

Example of a quantity of interest: the improvement ◮ Suppose that our objective is to minimize an unknown function f : X → R ◮ In our Bayesian approach, we choose a GP prior ξ for f (in other words, ξ is a model of f )

◮ Objective: construct a sequence (X 1 , X 2 , . . .) ∈ X that converges to X ⋆ = argmin x ξ(x )
◮ Given X 1 , . . . , X n+1 , how to define and choose a "good" point X n+1 in our setting?

◮ Let m n = min 1≤i≤n ξ(X i ) ◮ A "good" point x ∈ X is such that m n -ξ(x ) is large ◮ Define the excursion of ξ at x below m n , a.
k.a the improvement:

I n = 0 if ξ(x ) > m n m n -ξ(x ) if ξ(x ) ≤ m n 89/224
Example of a quantity of interest: the improvement

-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 ξ n (x ) +2 σ n (x ) -2 σ n (x ) 89/224
Example of a quantity of interest: the improvement

-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 89/224
Example of a quantity of interest: the improvement

-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 ξ n (x ) σ n (x ) 89/224
Example of a quantity of interest: the improvement

-1 -0.5 0 0.5 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Improvement ξ n (x ) m n 90/224
Example of a quantity of interest: the improvement ◮ Regions with high values of the posterior mean of I n are promising search regions for the minimum of ξ ◮ The posterior mean of I n may be written as Example of a quantity of interest: the improvement

ρ n (x ) = E (I n | ξ(X 1 ), . . . , ξ(X n )) = mn z=-inf (m n -z) p ξ(x )|ξ(X 1 ),...,ξ(Xn) (z) dz = γ(m n -ξ n (x ), σ 2 n (x )) with γ(z, s) =    √ s Φ ′ z √ s + z Φ z √ s if s > 0, max (z, 0) if s = 0. ◮ ρ n is
-1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 0 0.05 0.1 0.15 0.2 ξ n (x ) m n x ρ n (x )

92/224

Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 92/224
Insertion into an optimization algorithm ◮ The EI algorithm:

X n+1 = argmax x ρ n (x ) -1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 -6 -4 -2 0 x log 10 ρ n (x ) 93/224
Example of quantity of interest: the minimizer ◮ Assume an unknown function f : X → R and suppose we are interested in seeking its minimizer:

x ⋆ = argmin x f (x )
◮ Choose a GP prior ξ for f . Given observations ξ(x 1 ), . . . ξ(x n ), what is the posterior distrib. of X ⋆ = argmin x ξ(x )?

◮ Unlike I n above, the distrib. of X ⋆ does not possess a closed-form expression → resort to an empirical estimation using conditional sample paths

94/224

Empirical posterior density of the minimizer Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 -1 0 1 2 -1 -0.5 0 0.5 1 0 1 2 3 x p x ⋆ |Fn 94/224
Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 2 4 6 x p x ⋆ |F + n 94/224
Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 1 2 3 4 x p x ⋆ |F + n 94/224
Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 1 2 3 x p x ⋆ |F + n 94/224
Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 1 2 3 x p x ⋆ |F + n 94/224
Empirical posterior density of the minimizer

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 1 2 3 x p x ⋆ |F + n 94/224
Empirical posterior density of the minimizer Choosing a centered Gaussian random process

-1 -0.5 0 0.5 1 -2 0 2 -1 -0.5 0 0.5 1 0 1 2 3 4 x p x ⋆ |F + n 95/224
How to choose the covariance function of a GP ξ ∼ N (0, k)?

97/224

Regularity properties of a random process Def.

Given x 0 ∈ R d , a random process ξ is said to be continuous in mean-square at x 0 iff lim

x →x 0 E (ξ(x ) -ξ(x 0 )) 2 = 0
Prop.

Let ξ be a second-order random process with continuous mean function and stationary covariance function k. ξ is continuous in mean-square iff k is continuous at zero.

98/224

Regularity properties of a random process Def.

For x , h ∈ R d , define the random variable

ξ h (x ) = ξ(x 0 + h) -ξ(x 0 ) h ξ is mean-square differentiable at x 0 iff there exists a random vector ∇ξ(x 0 ) such that lim h→0 E ξ h (x 0 ) -(∇ξ(x 0 ), h) 2 = 0
Prop.

Let ξ be a second-order random process with differentiable mean function and stationary covariance function k. ξ is differentiable in mean-square iff k is two-time differentiable at zero.

99/224

Regularity properties of a random process ◮ Differentiability of the covariance function at the origin → mean-square differentiability of ξ 100/224

Influence of the regularity

mean-square continuity -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -4 -2 0 2 4 three-time mean-square differentiability -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 0.2 0.4 0.6 0.8 1 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 -3 -2 -1 0 1 2 101/224
Choice of a covariance ◮ A Gaussian process prior carries a high amount of information about f

→ it is often difficult to elicit such a prior before any evaluation is made ◮ Covariance function of ξ is usually assumed to belong to some parametric class of positive definite functions ◮ Parameter values assumed to be unknown ◮ Two approaches:

1. The parameters can be estimated from the evaluation results by maximum likelihood, and then used as if they were known (plug-in approach)

2. We can assume a prior distrib. for the parameters of the covariance and use a fully Bayesian approach

κ ν (h) = 1 2 ν-1 Γ(ν) 2ν 1/2 h ν K ν 2ν 1/2 h , h ∈ R (
k θ (h) = σ 2 κ ν (|h|/ρ) , h ∈ R (2) 103/224
Choice of a parametrized covariance function: the Matérn covariance

Matérn covariance in one dimension σ 2 = 1, ρ = 0.8 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 ν = 1/2 ν = 3/2 ν = 9/2 ξ is p-time mean-square differentiable iff ν > p 104/224
Choice of a parametrized covariance function: the Matérn covariance ◮ To model a function f defined over X ⊂ R d , with d > 1, we use the anisotropic form of the Matérn covariance:

k θ (x , y ) = σ 2 κ ν   d i=1 (x [i] -y [i] ) 2 ρ 2 i   , x , y ∈ R d (3)
where

x [i] , y [i]
denote the i th coordinate of x and y , and the positive scalars ρ i represent scale parameters d, in practice, we consider the vector of parameters

◮ Since σ 2 > 0, ν > 0, ρ i > 0, i = 1, . . . ,
θ = {log σ 2 , log ν, -log ρ 1 , . . . , -log ρ d } ∈ R d+2
→ makes parameter estimation easier

105/224

Parameter estimation by maximum likelihood ◮ Assume ξ is a zero-mean Gaussian process ◮ The log-likelihood of the data ξ n = (ξ(x 1 ), . . . , ξ(x n )) t can be written as

ℓ(ξ n ; θ) = - n 2 log(2π) - 1 2 log det K (θ) - 1 2 ξ n t K (θ) -1 ξ n , ( 4
)
where K (θ) is the covariance matrix of ξ n , which depends on the parameter vector θ ◮ The log-likelihood can be maximized using a gradient-based search method

106/224

Prediction of a Gaussian process with unknown mean function ◮ In the domain of computer experiments, the mean of a Gaussian process is generally written as a linear parametric function

m(•) = β t ϕ(•) , (5) with
β a vector of unknown parameters ϕ = (ϕ 1 , . . . , ϕ l ) t an l-dimensional vector of functions (in practice, polynomials)

◮ Simplest case: the mean function is an unknown constant m, in which case β = m and ϕ :

x ∈ X → 1 107/224
Prediction of a Gaussian process with unknown mean function ◮ Define the linear space of functions

P = x → l i=1 β i ϕ i (x ); β i ∈ R ,
◮ Define Λ the linear space of finite-support measures on X, i.e.

λ ∈ Λ =⇒ λ = n i=1 λ i δ xi for some n ∈ N ◮ For f : X → R, and λ = n i=1 λ i δ xi ∈ Λ, λ, f = X f dλ = n i=1 λ i f (x i )
◮ Define the linear subspace Λ P ⊥ ⊂ Λ of finite-support measures vanishing on P, i.e.

λ ∈ Λ P ⊥ =⇒ λ, f = X fdλ = n i=1 λ i f (x i ) = 0 , ∀f ∈ P 108/224
Prediction of a Gaussian process with unknown mean function ◮ Let ξ be a Gaussian random process with an unknown mean in P, and a covariance function k

◮ For x ∈ X, the (intrinsic) kriging predictor ξ n (x ) of ξ(x ) from ξ(x 1 ), . . . , ξ(x n ) is the linear projection ξ n (x ) = i λ i (x )ξ(x i )
of ξ(x ) onto span{ξ(x i ), i = 1, . . . , n} such that the variance of the error ξ(x )ξ n (x ) is minimized, under the constraint

δ x - λ i (x )δ xi ∈ Λ P ⊥ i.e., δ x - λ i (x )δ xi , ϕ j = ϕ j (x ) - λ i (x )ϕ j (x i ) = 0 , j = 1, . . . , l ◮ The requirement δ x -λ i (x )δ xi ∈ Λ P ⊥ makes the kriging predictor unbiased, even if the mean of ξ is unknown 109/224
Prediction of a Gaussian process with unknown mean function

ξ n (x ) is the linear projection of ξ(x ) onto span{ξ(x 1 ), . . . , ξ(x n )} orthogonally to P span{ξ(x i ), i ≤ n} P ξ(x ) ξ n (x ) O 110/224
Prediction of a Gaussian process with unknown mean function ◮ The weights λ i (x ; x n ) are again solutions of a system of linear equations, which can be written under a matrix form as

K ϕ t ϕ 0 λ(x ) µ(x ) = k(x ) ϕ(x ) , (6) with
ϕ an l × n matrix with entries

ϕ i (x j ), i = 1, . . . , l, j = 1, . . . , n -µ a vector of Lagrange coefficients -K , λ(x ), k(x ) as above 111/224
Prediction of a Gaussian process with unknown mean function ◮ When the mean is unknown, the kriging covariance function (the covariance of the error of prediction) is given by

k n (x , y ) := cov ξ(x ) -ξ n (x ), ξ(y ) -ξ n (y ) = k(x -y ) -λ(x ) t k(y ) -µ(x ) t ϕ(y ) .
Prop.

Let k be a covariance function and assume m ∈ P.

If ξ | m ∼ GP (m, k) m : x → β t ϕ(x ), β ∼ U(R l ) then ξ | F n ∼ GP ξ n (•), k n (•, •) with U(R l ) the (improper) uniform distribution over R l
→ justifies the use of kriging in a Bayesian framework provided that the covariance function of ξ is known

112/224

Parameter estimation with unknown mean function ◮ Objective: estimate the covariance parameters of a Gaussian process with unknown mean ◮ Restricted Maximum Likelihood (REML) approach → maximize the likelihood of the increments (or generalized increments) of the data ◮ Let ξ be a Gaussian process with an unknown mean function in P and ξ n the random vector of observations at points

x i , i = 1, . . . , n ◮ Let ϕ = (ϕ i (x j )) l,n i,j=1
be the l × n matrix of basis functions of P evaluated on {x 1 , . . . , x n }.

113/224

Parameter estimation with unknown mean function ◮ Since the dimension of P is l, the dimension of the space of the measures with support {x 1 , . . . , x n } that cancel out the functions of P is nl.

◮ Assume an n × (nl) matrix W with rank nl has been found, such that ϕW = 0 .

(The columns of W are in the kernel of ϕ.)

◮ Then Z = W T ξ n is a Gaussian random vector taking its values in R n-l , with zero mean and covariance matrix

W T K (θ)W
where K (θ) is the covariance matrix of ξ n with entries

k θ (x i -x j )
114/224 REML ◮ The random vector Z is a contrast vector ◮ The log-likelihood of the contrasts is given by

L(z | θ) = - n -l 2 log 2π- 1 2 log det(W t K (θ)W )- 1 2 z t (W t K (θ)W ) -1 z . 115/224

REML

◮ Various methods may be employed to compute the matrix W ◮ We favor the QR decomposition of ϕ T Decision-theoretic framework (cont'd)

ϕ T = (Q 1 | Q 2 ) R 0 , where (Q 1 | Q 2 ) is an n × n orthogonal matrix
◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

122/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

123/224

States of nature

Ω = f : X → R p | f such that . . . ◮ e.g., d = 1, p = 1, X = [0; 1] and Ω = C (X; R)
◮ Until further notice, we will use this simple (but important) setting to illustrate the basics of Bayesian optimization

125/224

Uncertainty quantification (reminder from Lecture #1)

◮ The true state of nature ω ⋆ ∈ Ω is unknown

◮ Example (cont'd) ◮ a function f ⋆ ∈ Ω = C (X; R) is inside the black box (ω ⋆ ≡ f ⋆ )
◮ we don't "know" f ⋆ (x ) until we run the code with input x

125/224

Uncertainty quantification (reminder from Lecture #1)

◮ The true state of nature ω ⋆ ∈ Ω is unknown

◮ Example (cont'd) ◮ a function f ⋆ ∈ Ω = C (X; R) is inside the black box (ω ⋆ ≡ f ⋆ )
◮ we don't "know" f ⋆ (x ) until we run the code with input x ◮ Bayesian approach to UQ ◮ our knowledge of ω ⋆ is encoded by a probability distribution on the set Ω of all possible ω's ◮ technically: proba on (Ω, F) for some σ-algebra F. . .

125/224

Uncertainty quantification (reminder from Lecture #1)

◮ The true state of nature ω ⋆ ∈ Ω is unknown

◮ Example (cont'd) ◮ a function f ⋆ ∈ Ω = C (X; R) is inside the black box (ω ⋆ ≡ f ⋆ )
◮ we don't "know" f ⋆ (x ) until we run the code with input x ◮ Bayesian approach to UQ ◮ our knowledge of ω ⋆ is encoded by a probability distribution on the set Ω of all possible ω's ◮ technically: proba on (Ω, F) for some σ-algebra F. . .

◮ Sequence of decisions ⇒ sequence of distributions P 0 , P 1 , . . . Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

128/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function) ◮ quasi-optimal region (sublevel set). . .

130/224

Decisions: standard setting and notations

◮ From now on we focus on the "standard" setting

= P 0 ( • | F n ) = P n-1 ( • | I n ) ◮ Various

Notations

X n (ω) = the n th evaluation point D N+1 (ω) = the "final decision" (estimate of the QoI)

D(ω) = (X 1 (ω), . . . , X N (ω), D N+1 (ω))
◮ We cannot use information that is not yet available ◮ X n (ω) depends on ω through F n-1 only ◮ D N+1 (ω) depends on ω through F N only ◮ D is a decision strategy (sequence of decision rules)

133/224

Transitions: examples Transitions: examples Transitions: examples Transitions: examples Transitions: examples ◮ availability of gradients (e.g., adjoint code) ◮ batch setting and/or multiple outputs ◮ variable run-times, "simulation failures". . .

◮ Example (cont'd) ◮ single-output, deterministic code ◮ I n = (X n , ξ(X n )) and F n = (X 1 , ξ(X 1 ), . . . , X n , ξ(X n )) ◮ ξ(X i ) = f ⋆ (X i ) is
◮ Example (cont'd) ◮ single-output, deterministic code ◮ I n = (X n , ξ(X n )) and F n = (X 1 , ξ(X 1 ), . . . , X n , ξ(X n )) ◮ ξ(X i ) = f ⋆ (X i ) is
◮ Example (cont'd) ◮ single-output, deterministic code ◮ I n = (X n , ξ(X n )) and F n = (X 1 , ξ(X 1 ), . . . , X n , ξ(X n )) ◮ ξ(X i ) = f ⋆ (X i ) is
◮ Example (cont'd) ◮ single-output, deterministic code ◮ I n = (X n , ξ(X n )) and F n = (X 1 , ξ(X 1 ), . . . , X n , ξ(X n )) ◮ ξ(X i ) = f ⋆ (X i ) is
◮ Example (cont'd) ◮ single-output, deterministic code ◮ I n = (X n , ξ(X n )) and F n = (X 1 , ξ(X 1 ), . . . , X n , ξ(X n )) ◮ ξ(X i ) = f ⋆ (X i ) is

134/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

134/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

135/224

Loss function ◮ To guide the decisions of the Bayesian agent, we need to specify a loss function L

Notation L : Ω × D → R (ω, d) → L(ω, d)
where D is the set of all possible sequences of decisions

135/224

Loss function ◮ To guide the decisions of the Bayesian agent, we need to specify a loss function L

Notation L : Ω × D → R (ω, d) → L(ω, d)
where D is the set of all possible sequences of decisions ◮ The Bayes-optimal strategy is, by definition:

D = argmin E 0 (L(D)) = argmin Ω L(ω, D(ω)) P 0 (dω)
where D ranges over all strategies ◮ Assume that we want to find the minimizer of f

L(f , d) = L(f , x ) = f ( x ) -min f
◮ d = (x1, . . . , xn, x )
◮ with x our estimate of argmin f ◮ A standard loss function for this situation is the linear loss:

L(f , d) = L(f , x ) = f ( x ) -min f (a.k.a. opportunity cost, a.k.a. instantaneous regret) ◮ Remarks ◮ L coincides with the L 1 loss of the estimator f ( x ) f ( x ) ≥ min f ⇒ L(f , x ) = |f ( x ) -min f |
◮ L is a terminal loss (does not depend on X 1 , ξ(X 1 ), . . .)

137/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

137/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function)

◮ Our BDT framework is complete, let's use it

137/224

Decision-theoretic framework (cont'd)

◮ How does this relate to optimization ?

◮ The agent is the optimization algorithm (or you, if you will)

Ingredients of a BDT problem ◮ a set of all possible "states of nature" ◮ a prior distribution over the states of nature ◮ a description of the decisions we have to make ◮ and the corresponding "transitions" ◮ a loss function (or utility function) ◮ Recall the Bayes-optimal strategy (algorithm):

D Bayes = argmin D E 0 (L(D N+1 )) = argmin D Ω L(ω, D N+1 (ω)) P 0 (dω)
where D ranges over all strategies D = (X 1 , . . . , X N , D N+1 )

140/224

Problem statement

What does this D Bayes look like ?

Can we actually build an optimal Bayesian algorithm ?

141/224

Problem statement (more precisely)

◮ Recall that ◮ X n+1 (ω) must depend on ω through F n only ◮ D N+1 (ω) must depend on ω through F N only 141/224
Problem statement (more precisely)

◮ Recall that ◮ X n+1 (ω) must depend on ω through F n only ◮ D N+1 (ω) must depend on ω through F N only
Notations

X n+1 = ϕ n (X 1 , Z 1 , . . . , X n , Z n ) = ϕ n (F n ) D N+1 = ϕ N (X 1 , Z 1 , . . . , X N , Z N ) = ϕ N (F N ) 141/224
Problem statement (more precisely)

◮ Recall that ◮ X n+1 (ω) must depend on ω through F n only ◮ D N+1 (ω) must depend on ω through F N only
Notations

X n+1 = ϕ n (X 1 , Z 1 , . . . , X n , Z n ) = ϕ n (F n ) D N+1 = ϕ N (X 1 , Z 1 , . . . , X N , Z N ) = ϕ N (F N )
◮ Goal: find the functions ϕ 0 , . . . , ϕ N

142/224

Lecture 2 : Bayesian optimization (BO) 2.2. From Bayes-optimal to myopic strategies

The optimal terminal decision Optimal choice of the last evaluation Bayes-optimal versus "practical Bayes" optimization Sampling criteria for multi-objective and/or contrained optimization

143/224

Optimal terminal decision Notation

E n = E ( • | F n ) = conditional

Optimal terminal decision

Notation

E n = E ( • | F n ) = conditional expectation with respect to F n
= expectation with respect to the probability P n ◮ Consider any incomplete strategy X 1 , . . . X N ◮ Claim: the optimal terminal decision is

D N+1 = ϕ Bayes N (X 1 , Z 1 , . . . , X N , Z N ) = argmin d E N (L(d))
where d runs over all possible values for the terminal decision

144/224

Optimal terminal decision: proof ◮ Take any strategy D = (X 1 , . . . , X N , D N+1 )

◮ Consider the modified strategy

D ′ = X 1 , . . . , X N , D ′ N+1 where D ′ N+1 = ϕ Bayes N (X 1 , Z 1 , . . . , X N , Z N ) 144/224
Optimal terminal decision: proof ◮ Take any strategy D = (X 1 , . . . , X N , D N+1 )

◮ Consider the modified strategy

D ′ = X 1 , . . . , X N , D ′ N+1 where D ′ N+1 = ϕ Bayes N (X 1 , Z 1 , . . . , X N , Z N ) ◮ Then, by definition of ϕ Bayes N , E N (L(D N+1 )) = E N (L(d)) |d=D N+1 ≥ min d E N (L(d)) = E N L(D ′ N+1 ) 144/224
Optimal terminal decision: proof ◮ Take any strategy D = (X 1 , . . . , X N , D N+1 )

◮ Consider the modified strategy

D ′ = X 1 , . . . , X N , D ′ N+1 where D ′ N+1 = ϕ Bayes N (X 1 , Z 1 , . . . , X N , Z N ) ◮ Then, by definition of ϕ Bayes N , E N (L(D N+1 )) = E N (L(d)) |d=D N+1 ≥ min d E N (L(d)) = E N L(D ′ N+1 ) ◮ and thus E 0 (L(D N+1 )) = E 0 (E N (L(D N+1 ))) ≥ E 0 E N L(D ′ N+1 ) = E 0 (L(D ′ N+1 ))

145/224

Vocabulary: posterior (Bayes) risk at time N ◮ Define the posterior risk at time N for the decision d:

R N (d) = E N (L(d))
("risk" is a synonym for "expected loss")

145/224

Vocabulary: posterior (Bayes) risk at time N ◮ Define the posterior risk at time N for the decision d:

R N (d) = E N (L(d))
("risk" is a synonym for "expected loss") Example (cont'd): linear loss ◮ Recall the setting

◮
◮ goal: minimize f ◮ d N+1 = x is an estimate of X * (f ) = argmin f ◮ L(f , x ) = f ( x ) -min f 146/224
Example (cont'd): linear loss ◮ Recall the setting

◮ goal: minimize f ◮ d N+1 = x is an estimate of X * (f ) = argmin f ◮ L(f , x ) = f ( x ) -min f
◮ Compute the posterior risk at time N for a given x ∈ X:

R N ( x ) = E N (L(ξ, x )) = E N (ξ( x ) -min ξ) = ξ N ( x ) -E N (min ξ) 146/224
Example (cont'd): linear loss ◮ Recall the setting

◮ goal: minimize f ◮ d N+1 = x is an estimate of X * (f ) = argmin f ◮ L(f , x ) = f ( x ) -min f
◮ Compute the posterior risk at time N for a given x ∈ X:

R N ( x ) = E N (L(ξ, x )) = E N (ξ( x ) -min ξ) = ξ N ( x ) -E N (min ξ)
◮ Thus the optimal terminal decision is

D Bayes N+1 = X Bayes = argmin x ∈X ξ N (x ) 147/224
Example (cont'd): linear loss

Assume that n = N = 5 (a small budget indeed). Example (cont'd): the L 1 loss, variant ◮ To summarize, we have for this example

X Bayes = argmin ξ N R Bayes N = min ξ N -E N (min ξ) ◮ Remark: in general, X Bayes ∈ {X 1 , . . . , X N }
◮ the value of the function at X Bayes is not known

148/224

Example (cont'd): the L 1 loss, variant ◮ To summarize, we have for this example

X Bayes = argmin ξ N R Bayes N = min ξ N -E N (min ξ) ◮ Remark: in general, X Bayes ∈ {X 1 , . . . , X N }
◮ the value of the function at X Bayes is not known

◮ Variant: restrict the terminal decision to {X 1 , . . . , X N } X Bayes,1 = argmin x ∈{X 1 ,...,X N } ξ(x ) R Bayes,1 N = min i≤N ξ(X i ) -E N (min ξ) 149/224
Example (cont'd): linear loss

Assume that n = N = 5 (a small budget indeed). Finding X Bayes N (last evaluation point)

◮ Let us focus now on the last evaluation point

◮ recall that D = (X 1 , . . . , X N-1 , X N , D N+1 ) 151/224
Finding X Bayes N (last evaluation point)

◮ Let us focus now on the last evaluation point

◮ recall that D = (X 1 , . . . , X N-1 , X N , D N+1 ) Notation E n,x (Y ) will mean: "compute E n (Y ), assuming that X n+1 = x " 151/224
Finding X Bayes N (last evaluation point)

◮ Let us focus now on the last evaluation point

◮ recall that D = (X 1 , . . . , X N-1 , X N , D N+1 ) Notation E n,x (Y ) will mean: "compute E n (Y ), assuming that X n+1 = x " ◮ For example, if Y = g (X 1 , Z 1 , . . . , X n , Z n , X n+1 , Z n+1 ), E n,x (Y ) = E n (g (X 1 , Z 1 , . . . , X n , Z n , x , Z x ))
where Z x denotes the result of a new evaluation at x

152/224

Finding X Bayes N (last evaluation point)

◮ Given x N ∈ X, consider the following strategy at time N -1:

1) first, evaluate at X N = x N , 2) then, optimally, i.e., use

D Bayes N+1 = ϕ Bayes N (F N ) 152/224
Finding X Bayes N (last evaluation point)

◮ Given x N ∈ X, consider the following strategy at time N -1:

1) first, evaluate at X N = x N , then, act optimally, i.e., use

D Bayes N+1 = ϕ Bayes N (F N ) ◮ The corresponding posterior risk at time N -1 is R N-1 (x N ) = E N-1,x N L(D Bayes N+1 ) = E N-1,x N R Bayes N 152/224
Finding X Bayes N (last evaluation point)

◮ Given x N ∈ X, consider the following strategy at time N -1:

1) first, evaluate at X N = x N , 2) then, act optimally, i.e., use

D Bayes N+1 = ϕ Bayes N (F N ) ◮ The corresponding posterior risk at time N -1 is R N-1 (x N ) = E N-1,x N L(D Bayes N+1 ) = E N-1,x N R Bayes N
◮ Claim: the optimal decision rule for the last evaluation is

X Bayes N = ϕ N-1 (F N-1 ) = argmin x N ∈X R N-1 (x N ) 152/224
Finding X Bayes N (last evaluation point)

◮ Given x N ∈ X, consider the following strategy at time N -1:

1) first, evaluate at X N = x N , 2) then, act optimally, i.e., use

D Bayes N+1 = ϕ Bayes N (F N ) ◮ The corresponding posterior risk at time N -1 is R N-1 (x N ) = E N-1,x N L(D Bayes N+1 ) = E N-1,x N R Bayes N
◮ Claim: the optimal decision rule for the last evaluation is

X Bayes N = ϕ N-1 (F N-1 ) = argmin x N ∈X R N-1 (x N )
◮ Remark: R N-1 is used as a "sampling criterion" (a.k.a. "infill criterion", a.k.a. "merit function". . . )

153/224

Finding X Bayes N : proof

◮ For any strategy D = (X 1 , . . . , X N-1 , X N , D N+1 ), E N-1 (L(D N+1 )) = E N-1 (R N (F N , D N+1 )) ≥ E N-1 R Bayes N (F N ) = R N-1 (F N-1 , X N-1 ) 153/224
Finding X Bayes N : proof

◮ For any strategy D = (X 1 , . . . , X N-1 , X N , D N+1 ), E N-1 (L(D N+1 )) = E N-1 (R N (F N , D N+1 )) ≥ E N-1 R Bayes N (F N ) = R N-1 (F N-1 , X N-1 ) ◮ Let D ′ = X 1 , . . . , X N-1 , X ′ N , D ′ N+1 , where X ′ N = ϕ Bayes N-1 (F N-1 ) and D ′ N+1 = ϕ Bayes N (F N-1 , X ′ N , Z ′ N ) 153/224
Finding X Bayes N : proof

◮ For any strategy D = (X 1 , . . . , X N-1 , X N , D N+1 ), E N-1 (L(D N+1 )) = E N-1 (R N (F N , D N+1 )) ≥ E N-1 R Bayes N (F N ) = R N-1 (F N-1 , X N-1 ) ◮ Let D ′ = X 1 , . . . , X N-1 , X ′ N , D ′ N+1 , where X ′ N = ϕ Bayes N-1 (F N-1 ) and D ′ N+1 = ϕ Bayes N (F N-1 , X ′ N , Z ′ N ) ◮ Then E N-1 L(D ′ N+1 ) = R Bayes N-1 (F N-1 ) ≤ R N-1 (F N-1 , X N-1 ) 153/224
Finding X Bayes N : proof

◮ For any strategy D = (X 1 , . . . , X N-1 , X N , D N+1 ), E N-1 (L(D N+1 )) = E N-1 (R N (F N , D N+1 )) ≥ E N-1 R Bayes N (F N ) = R N-1 (F N-1 , X N-1 ) ◮ Let D ′ = X 1 , . . . , X N-1 , X ′ N , D ′ N+1 , where X ′ N = ϕ Bayes N-1 (F N-1 ) and D ′ N+1 = ϕ Bayes N (F N-1 , X ′ N , Z ′ N ) ◮ Then E N-1 L(D ′ N+1 ) = R Bayes N-1 (F N-1 ) ≤ R N-1 (F N-1 , X N-1 ) ◮ Thus E 0 (L(D N+1 )) ≥ E 0 L(D ′ N+1 ) 154/224
Finding X Bayes N : example (cont'd)

◮ Recall our linear loss example

X Bayes = argmin ξ N R Bayes N = min ξ N -E N (min ξ) 154/224
Finding X Bayes N : example (cont'd)

◮ Recall our linear loss example

X Bayes = argmin ξ N R Bayes N = min ξ N -E N (min ξ) ◮ Compute the posterior risk at time N -1 R N-1 (F N-1 , x N ) = E N-1,x N R Bayes N (F N ) = E N-1,x N min ξ N -E N-1 (min ξ) 154/224
Finding X Bayes N : example (cont'd)

◮ Recall our linear loss example

X Bayes = argmin ξ N R Bayes N = min ξ N -E N (min ξ) ◮ Compute the posterior at time N -1 R N-1 (F N-1 , x N ) = E N-1,x N R Bayes N (F N ) = E N-1,x N min ξ N -E N-1 (min ξ) ◮
The optimal decision at time N -1 is

X N = argmin x N E N-1,x N min ξ N
(first appears (in english) in Mockus, Tiesis & Žilinskas, 1978) 155/224

Finding X Bayes N : example (cont'd)

◮ Equivalently,

X N = argmax x N min ξ N-1 -E N-1,x N min ξ N ρ KG N-1 (x N ) ≥ 0 155/224
Finding X Bayes N : example (cont'd)

◮ Equivalently,

X N = argmax x N min ξ N-1 -E N-1,x N min ξ N ρ KG N-1 (x N ) ≥ 0
◮ Nowadays called the Knowledge Gradient (KG) criterion (Frazier, Powell & co-authors, 2008, 2009, 2011) 155/224

Finding X Bayes N : example (cont'd)

◮ Equivalently,

X N = argmax x N min ξ N-1 -E N-1,x N min ξ N ρ KG N-1 (x N ) ≥ 0
◮ Nowadays called the Knowledge Gradient (KG) criterion (Frazier, Powell & co-authors, 2008, 2009, 2011) ◮ Remarks ◮ applicable to "noisy" observations as well Finding X Bayes N : example (cont'd)

◮ Equivalently,

X N = argmax x N min ξ N-1 -E N-1,x N min ξ N ρ KG N-1 (x N ) ≥ 0
◮ Nowadays called the Knowledge Gradient (KG) criterion (Frazier, Powell & co-authors, 2008, 2009, 2011) ◮ Remarks ◮ applicable to "noisy" observations as well ◮ a.k.a. simulation-based optimization ◮ even with a GP prior, ρ KG is not exactly computable in general ◮ idea: approx. max over a finite grid (more about that later)

156/224

Finding X Bayes N : example (cont'd)

Same example as before, n = 5, but assume now that N = 6. Finding X Bayes N : example, variant (cont'd)

◮ Recall the following variant

X Bayes,1 = argmin x ∈{X 1 ,...,X N } ξ(x ) R Bayes,1 N = min i≤N ξ(X i ) -E N (min ξ) 157/224
Finding X Bayes N : example, variant (cont'd)

◮ Recall the following variant

X Bayes,1 = argmin x ∈{X 1 ,...,X N } ξ(x ) R Bayes,1 N = min i≤N ξ(X i ) -E N (min ξ) ◮ Set M n = min i≤n ξ(X i ).
The optimal decision at time N -1 is

X N = argmax x N M N-1 -E N-1,x N (M N ) 157/224
Finding X Bayes N : example, variant (cont'd)

◮ Recall the following variant

X Bayes,1 = argmin x ∈{X 1 ,...,X N } ξ(x ) R Bayes,1 N = min i≤N ξ(X i ) -E N (min ξ) ◮ Set M n = min i≤n ξ(X i ).
The optimal decision at time N -1 is

X N = argmax x N M N-1 -E N-1,x N (M N ) = argmax x N E N-1 (M N-1 -ξ(x N )) + ρ EI n (x N ) ≥ 0
◮ This is the Expected Improvement (EI) criterion (Mockus et al 1978;Jones, Schonlau & Wlech, 1998) 157/224

Finding X Bayes N : example, variant (cont'd)

◮ Recall the following variant

X Bayes,1 = argmin x ∈{X 1 ,...,X N } ξ(x ) R Bayes,1 N = min i≤N ξ(X i ) -E N (min ξ) ◮ Set M n = min i≤n ξ(X i ).
The optimal decision at time N -1 is

X N = argmax x N M N-1 -E N-1,x N (M N ) = argmax x N E N-1 (M N-1 -ξ(x N )) + ρ EI n (x N ) ≥ 0
◮ This is the Expected Improvement (EI) criterion (Mockus et al 1978;Jones, Schonlau & Wlech, 1998) ◮ Computable analytically for GP priors ⇒ most commonly used (for deterministic numerical models)

158/224

Finding X Bayes N : example (cont'd)

Same example as before, n = 5, but assume now that N = 6. Lecture 2 : Bayesian optimization (BO) 2.2. From Bayes-optimal to myopic strategies

The optimal terminal decision Optimal choice of the last evaluation Bayes-optimal versus "practical Bayes" optimization Sampling criteria for multi-objective and/or contrained optimization

160/224

The Bayes-optimal strategy ◮ Recall the optimal terminal decision rule

ϕ Bayes N (F N ) = argmin d E N (L(d)) R Bayes N (F N ) = min d E N (L(d)) 160/224
The Bayes-optimal strategy ◮ Recall the optimal terminal decision rule

ϕ Bayes N (F N ) = argmin d E N (L(d)) R Bayes N (F N ) = min d E N (L(d))
◮ Recall the optimal rule for the last evaluation

ϕ Bayes N-1 (F N-1 ) = argmin x N E N-1,x N R Bayes N (F N ) R Bayes N-1 (F N-1 ) = min x N E N-1,x N R Bayes N (F N ) 161/224
The Bayes-optimal strategy ◮ The entire Bayes-optimal strategy can be similarly: ∀n,

ϕ Bayes n-1 (F n-1 ) = argmin xn E n-1,xn R Bayes n (F n ) R Bayes n-1 (F n-1 ) = min xn E n-1,xn R Bayes n (F n )
◮ This is called backward induction (or dynamic programming)

161/224

The Bayes-optimal strategy ◮ The entire Bayes-optimal strategy can be written similarly: ∀n,

ϕ Bayes n-1 (F n-1 ) = argmin xn E n-1,xn R Bayes n (F n ) R Bayes n-1 (F n-1 ) = min xn E n-1,xn R Bayes n (F n )
◮ This is called backward induction (or dynamic programming)

◮ So what ? Can we use this ?

162/224

The Bayes-optimal strategy ◮ More explicitely, the optimal decision for the first evaluation is

X 1 = argmin x 1 E 0,x 1 min x 2 E 1,x 2 . . . min x N E N-1,x N min d E N (L(d))
162/224

The Bayes-optimal strategy ◮ More explicitely, the optimal decision for the first evaluation is

X 1 = argmin x 1 E 0,x 1 min x 2 E 1,x 2 . . . min x N E N-1,x N min d E N (L(d))
◮ Very difficult to use in practice beyond N = 1 or 2

◮ each "min" is an optim. problem that needs to be solved. . . ◮ each "E n,x " is an integral that needs to be computed. [START_REF] Benassi | Nouvel algorithme d'optimisation bayésien utilisant une approche Monte-Carlo séquentielle[END_REF][START_REF] Feliot | A Bayesian approach to constrained single-and multi-objective optimization[END_REF] ◮ sample according to a well-chosen sequence of densities Multi-objective problems ◮ Several objective functions to be minimized:

f = (f 1 , . . . , f p ) ◮ f j : X → R, 1 ≤ j ≤ p Pareto domination relation z ≺ z ′ if (def)    z j ≤ z ′ j for all j ≤ p, z j < z ′ j
for at least one j ≤ p.

168/224

Multi-objective problems ◮ Several objective functions to be minimized:

f = (f 1 , . . . , f p ) ◮ f j : X → R, 1 ≤ j ≤ p Pareto domination relation z ≺ z ′ if (def)    z j ≤ z ′ j for all j ≤ p, z j < z ′ j
for at least one j ≤ p.

◮ The goal is to find (estimate)

◮ the Pareto set P = {x ∈ X : ∃x ′ ∈ X, f (x ′ ) ≺ f (x )} (a.k.a. set of Pareto-efficient solutions) 168/224
Multi-objective problems ◮ Several objective functions to be minimized:

f = (f 1 , . . . , f p ) ◮ f j : X → R, 1 ≤ j ≤ p Pareto domination relation z ≺ z ′ if (def)    z j ≤ z ′ j for all j ≤ p, z j < z ′ j
for at least one j ≤ p.

◮ The goal is to find (estimate) Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

◮ the Pareto set P = {x ∈ X : ∃x ′ ∈ X, f (x ′ ) ≺ f (x )} (
f 1 f 2 z 1 z 2

Noiseless evaluations

169/224

Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

f 1 f 2 z 1 z 2 z ref Noiseless evaluations B = Π p j=1 -∞; z ref j : bounding box 169/224
Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

H ⋆ f 1 f 2 z 1 z 2 z ref Noiseless evaluations B = Π p j=1 -∞; z ref j : bounding box True dominated region: H ⋆ (f ) = {z ∈ B, ∃x ∈ X, f (x ) z} 169/224
Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

H ⋆ f 1 f 2 z 1 z 2 z ref Noiseless evaluations B = Π p j=1 -∞; z ref j : bounding box True dominated region: H ⋆ (f ) = {z ∈ B, ∃x ∈ X, f (x ) z}
Loss function:

L(f , H) = |H ⋆ (f )△ H| 169/224
Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

H ⋆ f 1 f 2 z 1 z 2 z ref Noiseless evaluations B = Π p j=1 -∞; z ref j : bounding box True dominated region: H ⋆ (f ) = {z ∈ B, ∃x ∈ X, f (x ) z}
Loss function:

L(f , H) = |H ⋆ (f )△ H|
Best "safe" estimator:

Hn = {z ∈ B, ∃i ≤ n, f (Xi ) z} 169/224
Multi-objective problems ◮ EHVI: a natural extension of EI [START_REF] Emmerich | Single-and multi-objective evolutionary optimization assisted by Gaussian random field metamodels[END_REF])

H ⋆ f 1 f 2 z 1 z 2 z 3 z ref Noiseless evaluations B = Π p j=1 -∞; z ref j : bounding box True dominated region: H ⋆ (f ) = {z ∈ B, ∃x ∈ X, f (x ) z}
Loss function:

L(f , H) = |H ⋆ (f )△ H|
Best "safe" estimator:

Hn = {z ∈ B, ∃i ≤ n, f (Xi ) z} ρ EHVI n (x n+1 ) = E n,x n+1 ( |H n+1 \ H n |) 170/224
Multi-objective problems Inequality-constrained problems ◮ Single-objective, inequality-contrained problem:

◮ f = (f o , f c,1 , . . . , f c,q ), with ◮ f o : X → R, to be minimized, ◮ f c,j : X → R, 1 ≤ j ≤ q, must be ≤ 0.
◮ Consider the following loss function

L(f , x ) =    f o ( x ) -f ⋆ o if f c ( x ) ≤ 0, +∞ otherwise.
where

f ⋆ o = min x :fc(x )≤0 f o (x ) 172/224
Inequality-constrained problems ◮ Assuming ◮ noiseless evaluations, ◮ independent priors on objective and constraint functions,

◮ ∃i ≤ n, ξ c (X i ) = f c (X i ) ≤ 0,
the following myopic criterion follows [START_REF] Schonlau | The correlated knowledge gradient for simulation optimization of continuous parameters using gaussian process regression[END_REF])

ρ EIC n (x n+1 ) = ρ EI o,n (x n+1 ) • Π q j=1 P n (ξ c,j (x n+1 ) ≤ 0)
Proba of Feasibility (PF)

.

172/224

Inequality-constrained problems ◮ Assuming ◮ noiseless evaluations, ◮ independent priors on objective and constraint functions,

◮ ∃i ≤ n, ξ c (X i ) = f c (X i ) ≤ 0,
the following myopic criterion follows [START_REF] Schonlau | The correlated knowledge gradient for simulation optimization of continuous parameters using gaussian process regression[END_REF])

ρ EIC n (x n+1 ) = ρ EI o,n (x n+1 ) • Π q j=1 P n (ξ c,j (x n+1 ) ≤ 0)
Proba of Feasibility (PF)

.

◮ Implementation ◮ Easy for independent GP priors (most commonly used) ◮ Dependent priors: harder. . . (but see [START_REF] Williams | Sequential design of computer experiments for constrained optimization[END_REF] 172/224

Inequality-constrained problems ◮ Assuming ◮ noiseless evaluations, ◮ independent priors on objective and constraint functions,

◮ ∃i ≤ n, ξ c (X i ) = f c (X i ) ≤ 0,
the following myopic criterion follows [START_REF] Schonlau | The correlated knowledge gradient for simulation optimization of continuous parameters using gaussian process regression[END_REF]) (a few words on the) Worst-case approach ◮ Principle of the worst-case (minimax) approach (a few more words on the) Worst-case approach ◮ In this lecture, we will focus on the probabilistic approach 181/224

ρ EIC n (x n+1 ) = ρ EI o,n (x n+1 ) • Π q j=1 P n (ξ c,j (x n+1 ) ≤ 0)
(a few more words on the) Worst-case approach ◮ In this lecture, we will focus on the probabilistic approach ◮ See [START_REF] Marzat | Worst-case global optimization of black-box functions through kriging and relaxation[END_REF], 2016) for a "BO treatment" of the worst-case approach (using relaxation)

181/224

(a few more words on the) Worst-case approach ◮ In this lecture, we will focus on the probabilistic approach ◮ See [START_REF] Marzat | Worst-case global optimization of black-box functions through kriging and relaxation[END_REF], 2016) for a "BO treatment" of the worst-case approach (using relaxation)

◮ An issue of terminology: in the math literature,

◮ "robust optimization" refers mainly to the worst-case setting (see Ben [START_REF] Ben-Tal | Robust optimization[END_REF], [START_REF] Bertsimas | Reliability-based design optimization using efficient global reliability analysis[END_REF] and refs)

◮ the probabilistic approach is called stochastic programming ◮ while engineers use the word "robust" for both

182/224

The probabilistic approach

◮ From now, we focus on the probabilistic approach ◮ u is considered as random → U ∼ P U ◮ can be a random vector (∈ R m ), or a more complicated object

182/224

The probabilistic approach

◮ From now, we focus on the probabilistic approach ◮ u is considered as random → U ∼ P U ◮ can be a random vector (∈ R m ), or a more complicated object ◮ Numerical models: two important settings ◮ stochastic simulators ◮ environmental variables

x RNG code Z = f (x , U) x code u z = f (x , u)
where f = (f o,1 , . . . , f o,p , f c,1 , . . . , f c,q )

183/224

The "stochastic simulator" setting

x RNG code Z = f (x , U)
◮ Features of the black box ◮ U is not directly accessible ◮ only x can be chosen by the algorithm; Z = f (x , U) is observed ◮ P U is not known explicitely

183/224

The "stochastic simulator" setting

x RNG code Z = f (x , U)
◮ Features of the black box ◮ U is not directly accessible ◮ only x can be chosen by the algorithm; Z = f (x , U) is observed ◮ P U is not known explicitely

◮ "State of nature" (the things that we don't know)

◮ the family P Z x of conditional distributions ◮ the RV U 1 , U 2 . . . that will be generated when running the computer model with inputs x 1 , x 2 . . . Example 3: Three-bar truss [START_REF] Koski | Defectiveness of weighting method in multicriterion optimization of structures[END_REF][START_REF] Das | Nonlinear multicriteria optimization and robust optimality[END_REF])

P w D L 190/224
Example 3: Three-bar truss [START_REF] Koski | Defectiveness of weighting method in multicriterion optimization of structures[END_REF][START_REF] Das | Nonlinear multicriteria optimization and robust optimality[END_REF] Example 3: Three-bar truss [START_REF] Koski | Defectiveness of weighting method in multicriterion optimization of structures[END_REF][START_REF] Das | Nonlinear multicriteria optimization and robust optimality[END_REF] F 1 (wind)

F 2 (suspended load)

P w D L a 1 a 2 a 3
◮ Design variables: x = (a 1 , a 2 , a 3 , w )

◮ a j : cross-section of bar j ◮ Environmental variables: U = (F 1 , F 2 )

◮

  Poor fit!◮ Why? Model capacity is weak ◮ Now, as an example, consider the modelf n (x ) = θ t r (x ) with r (x ) = (1 cos(2πx ) sin(2πx ) . . . cos(2mπx ) sin(2mπx )) t ∈ R 2m+1(a truncated Fourier series) ◮ The increase in the number of parameters yields an ill-defined problem (l ≫ n)

◮◮◮◮◮(

  The regularization principle alone is not enough to obtain a good approximation As modeling capacity increases, overfitting may arise 20/224 A simple curve fitting problem ◮ To avoid overfitting, we should try a regularization that penalizes high frequencies more For fitting problem ◮ n = 8, m = 50, l = 101, C = 10 -8 , α = 1Using a random generator, it is possible to "generate" sample values z 1 , z 2 , . . . of our model Z → possible values for our uncertain quantity of interest Using a random generator, it is possible to "generate" sample values z 1 , z 2 , . . . of our model Z → possible values for our uncertain quantity of interest as n → ∞, the empirical distribution of the realizations tends to the normal distribution)

◮

  Here, ξ(ω, •) = θ(ω) t r (•)◮ Fixing ω (a sample path) amounts to "choosing" a value for the random vector θ k) α , k = 1, . . . , mwith ω 0 = 2π 10 , k) α , k = 1, . . . , m with ω 0 = 2π 50 , k) α , k = 1, . . . , mwith ω 0 = 2π 10 , α = 1

  Choice of a parametrized covariance function: the Matérn covariance ◮ The Matérn covariance function is a conventional covariance function in the literature of computer experiments → offers the possibility to adjust the regularity of ξ with a single parameter ◮ The Matérn function:

  Gamma function -K ν the modified Bessel function of the second kind ◮ To model a real-valued function defined over X ⊂ R, we use the Matérn covariance:

◮

  P n corresponds to the agent's beliefs after the n th decision ◮ A prior distribution P 0 needs to be specified 126/224 Uncertainty quantification: example ◮ Example (cont'd): if f is known to look more or less like this: take P 0 = GP(m, k) with m ∼ U(R) and k a (stationary) Matérn covariance 126/224 Uncertainty quantification: example ◮ Example (cont'd): if f is known to look more or less like this: take P 0 = GP(m, k) with m ∼ U(R) and k a (stationary) Matérn covariance ◮ Gaussian process priors are commonly used because ◮ they are computationally convenient ◮ while allowing a certain modeling flexibility 127/224 Uncertainty quantification: consequences ◮ For clarity, consider again the case of a deterministic model: ◮ an unknown function f ∈ Ω is in the black box 127/224 Uncertainty quantification: consequences ◮ For clarity, consider again the case of a deterministic model: ◮ an unknown function f ∈ Ω is in the black box ◮ Given a proba P on Ω, we can ◮ compute the probability of any (measurable) statement about f ◮ compute the expectation of any (measurable) function of f i.e., the unknown f can be treated as random 127/224Uncertainty quantification: consequences ◮ For clarity, consider again the case of a deterministic model:◮ an unknown function f ∈ Ω is in the black box◮ Given a proba P on Ω, we can ◮ compute the probability of any (measurable) statement about f ◮ compute the expectation of any (measurable) function of f i.e., the unknown f can be treated as randomConvenient notationξ = random function that represents the unknown f ◮ we will write, e.g., E n (ξ(x )) instead of Ω f (x ) P n (df )128/224Decision-theoretic framework (cont'd)◮ How does this relate to optimization ?◮ The agent is the optimization algorithm (or you, if you will)

◮◮

  Final decision (ex: single-objective minimization pb.)◮ an estimate of the minimizer x * = argmin x ∈X f (x ) ◮ and/or an estimate of the minimum M * = min x ∈X f (x ) with X a "known" input space Several types of decisions in an optimization procedure: ◮ intermediate decisions ◮ stopping decision ◮ final decision ◮ Final decision (ex: single-objective minimization pb.) ◮ an estimate of the minimizer x * = argmin x ∈X f (x ) ◮ and/or an estimate of the minimum M * = min x ∈X f (x ) with X a "known" input space ◮ Other settings ◮ multi-objective: Pareto set / Pareto front (see later), ◮ inequality constraints (see later), equality constraints (harder !) ,

  the true, scalar, value of the model at X i ◮ Many other (interesting) settings are possible ! ◮ stochastic simulators ◮ the output is a random draw Zn ∼ some distrib. P Zn|Xn ◮ In = (Xn, Zn) and Fn = (X1, Z1, . . . ,Xn, Zn)

  the true, scalar, value of the model at X i ◮ Many other (interesting) settings are possible ! ◮ stochastic simulators ◮ the output is a random draw Zn ∼ some distrib. P Zn|Xn ◮ In = (Xn, Zn) and Fn = (X1, Z1, . . . ,Xn, Zn)

  (a.k.a. opportunity cost, a.k.a. instantaneous regret)

  expectation with respect to F n = expectation with respect to the probability P n 143/224

  Define the posterior Bayes risk at time N:

  Same example as before, n = 5, but assume now that N = 6. Same example as before, n = 5, but assume now that N = 6.Warning: X N = argmax ξ N-1 (uncertainty is taken into account)

  Warning: X N = argmax ξ N-1 (uncertainty is taken into account) Same example as before, n = 5, but assume now that N = 6.Warning: X N = argmax ξ N-1 (uncertainty is taken into account)

  a.k.a. set of Pareto-efficient solutions) ◮ and/or the Pareto front {z ∈ R p : ∃x ∈ P, z = f (x )} (a.k.a Pareto frontier, Pareto boundary. . . )

◮

  Define an uncertainty set U ◮ Optimize by considering the worst u ∈ U ◮ For instance, assuming a single-objective problem: minimize max u∈U f o (x , u)◮ If the problem has constraints, they become:∀j ≤ q, ∀u ∈ U, f c,j (x , u) ≤ 0 179/224Example 1: Illustration of the worst-case approachExample: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 Illustration of the worst-case approach Example: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 Illustration of the worst-case approach Example: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 fo(x , u) argmin f + oRemark: very conservative, the nominal performance is ignored179/224Example 2: Illustration of the worst-case approachExample: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 Illustration of the worst-case approach Example: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5Another example: worst-case approach for a constraintExample: f c (x , u) = x 2 -Another example: worst-case approach for a constraintExample: f c (x , u) = x 2u 2 , with u ∈ U = [u 0δ; u 0 + δ]

◮◮◮◮◮

  the RV U real that defines the f (x , U real ) ultimately realized single objective f = f o , expensive to evaluate ◮ no (expensive-to-evaluate) constraints ◮ remark: cheap constraints allowed in the definition of X ⊂ R d ◮ "environmental variables" setting ◮ Consider once again the L 1 loss functionL ((f , u real ), x ) = f ( x , u real ) -min x f (x , u real ) = f ( x , u real ) -min x f (x , u real ) Compute the posterior risk at time N for an estimate x ∈ X E N (L ((ξ, U real ), x )) = E N (ξ( x , U real )) -E N (min ξ(•, U real )) Compute the posterior risk at time N for an estimate x ∈ X E N (L ((ξ, U real ), x )) = E N (ξ( x , U real )) -E N (min ξ(•, U real )) = E N ξ( x ) -E N (min ξ(•, U real ))where ξ(x ) = ξ(x , u) P U (du) Compute the posterior risk at time N for an estimate x ∈ XE N (L ((ξ, U real ), x )) = E N (ξ( x , U real )) -E N (min ξ(•, U real )) = E N ξ( x ) -E N (min ξ(•, U real ))where ξ(x ) = ξ(x , u) P U (du) ◮ Same L 1 risk (ignoring last term) as if we were dealing with the Equivalent "deterministic" problemmin x f (x ), with f (x ) = f (x , u) P U (du)(Remark: this formulation occurs very naturally in a BO framework ) Illustration of the worst-case approachExample: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 fo(x , u) argmin f + oRemark: very conservative, the nominal performance is ignored189/224Example 1: Worst-case versus probabilistic approachExample: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 f o = E (f o (•, U)), with U ∼ N (0, s 2 ), s.t. P (|U| ≤ δ) = 99.9% 189/224Example 2: Illustration of the worst-case approach Example: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 Worst-case versus probabilistic approach Example: f o (x , u) = f (x + u), with u ∈ U = [-δ; δ], δ = 5 fo(x , u) argmin f + o f o argmin f o f o = E (f o (•, U)), with U ∼ N (0, s 2 ), s.t. P (|U| ≤ δ) = 99.9% 190/224

◮

  Design variables: x = (a 1 , a 2 , a 3 , w ) ◮ a j : cross-section of bar j 190/224

  and ∀a 1 , . . . , a n ∈ R the random variable

	i a i X i is Gaussian
	◮ In what follows, assume that G is centered, i.e., each element
	in G is a zero-mean random variable
	69/224

78/224 Kriging equations proof 79/224

	matrix of the observation vector
	-k(x ): n × 1 vector with entries k(x i , x )

80/224

	Software for kriging/GP regression Software for kriging/GP regression
	◮ Matlab/GNU Octave
	◮	DACE DACE, a matlab kriging toolbox.
	◮	FERUM Finite Element Reliability using Matlab
	◮ R packages ◮ GPML Gaussian Processes for Machine Learning
	◮ ◮	BACCO Bayesian analysis of computer code software GPStuff GP Models for Bayesian analysis
	◮ ◮	fanovaGraph Building Kriging Models from FANOVA Graphs scalaGAUSS Kriging toolbox with a focus on large datasets
	◮ ◮	DiceKriging DiceOptim Matlab Stat & ML toolbox GP regression from Mathworks GPareto Dice and ReDice packages
	◮ ◮	MuFiCokriging Multi-Fidelity Cokriging models STK Small (Matlab/GNU Octave) Toolbox for Kriging
	◮ ◮	RobustInv Robust inversion based on GP (like KrigInv) SUMO ooDACE Surrogate Modeling Lab
	◮ ◮	tgp Treed Gaussian processes UQLab Uncertainty quantification framework in Matlab
	◮ Python
	◮	scikit-learn Machine learning in Python
	◮	OpenTURNS Open source lib for UQ
	◮	Spearmint Bayesian optimization
	◮	GPy Gaussian processes framework in Python
		81/224 82/224

◮ . . .

Forrester, A. Sóbester and A. Keane, Engineering

	Books
	Lecture 1 : From meta-models to UQ
	1.1 Introduction
	1.2 Black-box modeling ◮ M. Stein, Interpolation of Spatial Data: Some Theory for 1.3 Bayesian approach kriging, Springer, 1999 1.4 Posterior distribution of a quantity of interest
	1.5 Complements on Gaussian processes
	Lecture 2 : Bayesian optimization (BO)
	2.1. Decision-theoretic framework
	2.2. From Bayes-optimal to myopic strategies
	design via surrogate modelling: a practical guide, John Wiley 2.3. Design under uncertainty
	& Sons, 2008
	116/224 117/224 122/224

and R is a l × l upper triangular matrix ◮ It is trivial to check that the columns of Q 2 form a basis of the kernel of ϕ ◮ So we may chose W = Q 2 ◮ Note that W T W = I n-l . ◮ T. Santner, B. Williams and W. Notz, The Design and Analysis of Computer Experiments, Springer, 2003 ◮ C. Rasmussen and C. Williams, Gaussian processes for Machine Learning, MIT Press, 2006 ◮ A.

  BO can deal with the stopping decision too ◮ in practice, difficult to translate into a loss (see later) ◮ some "BO papers" propose heuristic stopping rule

	Decisions Decisions: intermediate decisions Decisions: intermediate decisions Decisions: stopping decision Decisions: stopping decision Decisions: final decision
	◮ Several types of decisions in an optimization procedure: ◮ Several types of decisions in an optimization procedure: ◮ Several types of decisions in an optimization procedure: ◮ Several types of decisions in an optimization procedure: ◮ Several types of decisions in an optimization procedure: ◮ Several types of decisions in an optimization procedure:
	◮ intermediate decisions ◮ intermediate decisions ◮ intermediate decisions ◮ intermediate decisions ◮ intermediate decisions ◮ intermediate decisions
	◮ stopping decision ◮ stopping decision ◮ stopping decision ◮ stopping decision ◮ stopping decision ◮ stopping decision
	◮ final decision ◮ final decision ◮ final decision ◮ final decision ◮ final decision ◮ final decision
	◮ parallel computing: batches of input values
	◮ multi-fidelity: choosing the right fidelity level
	◮ variable run-time: choosing when to stop a computation
	129/224 129/224 129/224 129/224 129/224 129/224

◮ Intermediate decisions (simple setting)

◮ running the numerical model with a given input x ∈ X ◮ getting the corresponding output (deterministic or random)

◮ Intermediate decisions (simple setting)

◮ running the numerical model with a given input x ∈ X ◮ getting the corresponding output (deterministic or random)

◮ Intermediate decisions (various extensions) ◮ . . .

◮ Stopping decision (standard setting)

◮ a budget of evaluations, or computation time, is given ◮ the stopping decision is trivial in this case ◮ Stopping decision (standard setting) ◮ a budget of evaluations, or computation time, is given ◮ the stopping decision is trivial in this case ◮ Digression: taking the cost of observations into account?

◮ in principle,

  the agent's knowledge at time n is described by P n ◮ the (n + 1) th decision induces a transition P n → P n+1 the agent's knowledge at time n is described by P n ◮ the (n + 1) th decision induces a transition P n → P n+1

	Decisions: standard setting and notations Transitions: conditioning probability measures Transitions: conditioning probability measures
	130/224 131/224 131/224

◮ intermediate decisions ≡ evaluations (known comput. cost) ◮ stopping: a budget of N evaluations is given ◮ From now on we focus on the "standard" setting ◮ intermediate decisions ≡ evaluations (known comput. cost) ◮ stopping: a budget of N evaluations is given Notations X n (ω) = the n th evaluation point D N+1 (ω) = the "final decision" (estimate of the QoI)

D(ω) = (X 1 (ω), . . . , X N (ω), D N+1 (ω)) ◮ Recall that ◮ ◮ Recall that

◮ Notations: available information I n = the information obtained as a result of the n th decision F n = the information available at time n ◮ Mathematically, P n

  the true, scalar, value of the model at X i

	133/224

  the true, scalar, value of the model at X i

◮ Many other (interesting) settings are possible ! ◮ stochastic simulators ◮ the output is a random draw Zn ∼ some distrib. P Zn|Xn ◮ In = (Xn, Zn) and Fn = (X1, Z1, . . . , Xn, Zn) ◮ availability of gradients (e.g., adjoint code)

133/224

  the true, scalar, value of the model at X i

◮ Many other (interesting) settings are possible ! ◮ stochastic simulators ◮ the output is a random draw Zn ∼ some distrib. P Zn|Xn ◮ In = (Xn, Zn) and Fn = (X1, Z1, . . . , Xn, Zn) ◮ availability of gradients (e.g., adjoint code) ◮ batch setting and/or multiple outputs 133/224

  . . ◮ none of them are tractable, even for the nicest (GP) priors

	Practical Bayesian optimization: optimization Practical Bayesian optimization: optimization Practical Bayesian optimization: optimization
	◮ Each iteration involves an auxiliary optimization problem ◮ Each iteration involves an auxiliary optimization problem ◮ Each iteration involves an auxiliary optimization problem
	◮ Various approaches to solve it ◮ Various approaches to solve it ◮ Various approaches to solve it
	◮ Fix grid or IID random search ◮ Fix grid or IID random search ◮ Fix grid or IID random search
	◮ OK for low-dimensional, simple problems ◮ OK for low-dimensional, simple problems ◮ OK for low-dimensional, simple problems
	◮ if accurate convergence is not needed ◮ if accurate convergence is not needed ◮ if accurate convergence is not needed
	◮ External solvers ◮ External solvers ◮ External solvers
	◮ ex: DiceOptim → Rgenoud (genetic + gradient) ◮ ex: DiceOptim → Rgenoud (genetic + gradient) ◮ ex: DiceOptim → Rgenoud (genetic + gradient)
	◮ ex: Janusvekis & Le Riche (2013) → CMA-ES ◮ ex: Janusvekis & Le Riche (2013) → CMA-ES ◮ ex: Janusvekis & Le Riche (2013) → CMA-ES
	◮ Sequential Monte Carlo (Benassi, 2013; Feliot et al, 2017)
	◮ sample according to a well-chosen sequence of densities
	166/224 166/224 166/224

◮ Sequential Monte Carlo

◮ n = 8, m = 50, l = 101, C = 10 -8

Assume that n = N = 5 (a small budget indeed). Remark: the two estimates are equal when ξ does not "overshoot" (e.g., for a Brownian motion prior)

150/224

Lecture 2 : Bayesian optimization (BO) 2.2. From Bayes-optimal to myopic strategies

The optimal terminal decision Optimal choice of the last evaluation Bayes-optimal versus "practical Bayes" optimization Sampling criteria for multi-objective and/or contrained optimization

The "environmental variables" setting

◮ Features of the black box ◮ the simulator remains deterministic ◮ P U is specified separately, usually explicitely ◮ the algorithm can choose (x , u) pairs to be evaluated

184/224

The "environmental variables" setting

◮ Features of the black box ◮ the simulator remains deterministic ◮ P U is specified separately, usually explicitely ◮ the algorithm can choose (x , u) pairs to be evaluated

◮ "State of nature" (the things that we don't know)

◮ the deterministic function f : (x , u) → f (x , u) ◮ the RV U real that defines the f (x , U real ) ultimately realized

184/224

The "environmental variables" setting

◮ Features of the black box ◮ the simulator remains deterministic ◮ P U is specified separately, usually explicitely ◮ the algorithm can choose (x , u) pairs to be evaluated

◮ "State of nature" (the things that we don't know)

◮ The two settings can be mixed

185/224

Problem formulations ◮ Various "robust" formulations can be considered for the design problem, depending mainly on ◮ the number of objective functions, ◮ the presence of (expensive-to-evaluate) constraints, ◮ and, of course, how we want to deal with U real .

Problem formulations

◮ Various "robust" formulations can be considered for the design problem, depending mainly on ◮ the number of objective functions, ◮ the presence of (expensive-to-evaluate) constraints, ◮ and, of course, how we want to deal with U real .

◮ In the following, we focus on ◮ single objective problems ◮ in the "environmental variables" setting ◮ and discuss two important cases:

◮ optimization of the averaged objective function ◮ reliability-based design optimization (RBDO), a.k.a. "chance constrained" optimization, and other formulations

Prior model

◮ There are two functions of interest in this setting ◮ the one that can be observed, i.e., f : (x , u) → f (x , u), ◮ and the one that want to optimize:

◮ f is a function of f ◮ priors cannot be specified independently

Gaussian process priors are, again, very convenient

Proof ◮ ξ is Gaussian by linearity of the integral ◮ Computation of the mean function: exchange and E

◮ Computation of the covariance function: idem with bilinearity

Prior model (cont'd)

◮ Actually, we can say much better:

Jointly Gaussian processes

If ξ ∼ GP(m, k), then ξ and ξ are jointly Gaussian, and

Prior model (cont'd)

◮ Actually, we can say much better:

Jointly Gaussian processes

If ξ ∼ GP(m, k), then ξ and ξ are jointly Gaussian, and

◮ Remark: m ξ , k ξ and k ξ,ξ can be computed exactly [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] for exact formulas & approximations

195/224

Prior model (cont'd)

◮ Actually, we can say much better:

Jointly Gaussian processes

If ξ ∼ GP(m, k), then ξ and ξ are jointly Gaussian, and

◮ Remark: m ξ , k ξ and k ξ,ξ can be computed exactly [START_REF] Girard | Approximate methods for propagation of uncertainty with gaussian process models[END_REF] for exact formulas & approximations

Examples 1 and 2: discretization of P U ◮ In the two examples, Prob U = N (0, s 2 ), with s = 1.52

◮ We choose to use a regular discretization with n U = 11 points ◮ points regularly spaced on [-5; 5] ◮ weights computed using the normal cdf (using mid-points) Example 3: ◮ What decision(s) do we have to make at each step ?

◮ i.e, what do we need to provide to run the numerical model ?

Sampling strategy: what ?

◮ What decision(s) do we have to make at each step ?

◮ i.e, what do we need to provide to run the numerical model ?

◮ in this case, we must simply select a point X n+1 ∈ X ◮ In the following slides we assume the general case (adaptation to the special case poses no difficulty)

199/224

Sampling strategy: how ?

◮ How do we build a sampling strategy for this problem ?

◮ in this lecture, we will apply the standard BO machinery

Sampling strategy: how ?

◮ How do we build a sampling strategy for this problem ?

◮ in this lecture, we will apply the standard BO machinery Sampling strategy: how ?

◮ How do we build a sampling strategy for this problem ?

◮ in this lecture, we will apply the standard BO machinery ◮ Assume now the L 1 loss ◮ Recall the posterior risk at time N for an estimate x ∈ X:

Sampling strategy: one step look-ahead ◮ Let L n denote the expected loss that we would get if we stopped at time n:

Sampling strategy: one step look-ahead ◮ Let L n denote the expected loss that we would get if we stopped at time n:

◮ The one-step look-ahead (myopic) strategy is

Sampling strategy: one step look-ahead ◮ Let L n denote the expected loss that we would get if we stopped at time n:

◮ The one-step look-ahead (myopic) strategy is

Sampling strategy: one step look-ahead

where ρ n denotes the corresponding "expected improvement"

Sampling strategy: one step look-ahead

where ρ n denotes the corresponding "expected improvement"

◮ Formally, looks like the KG criterion of Frazier & co, but. . .

202/224

Sampling strategy: one step look-ahead ◮ There is no real difference mathematically: in both cases 1. the function to be optimized is not observable directly, 2. the evaluation results and the function to be optimized are jointly Gaussian.

◮ Good news: we can then derive an implementable Approximate KG criterion as in Scott et al ( 2011)

203/224

Approximate KG criterion (AKG)

◮ Let X ref n ⊂ X denote some finite "reference set"

The AKG criterion is:

where the min runs over X ref n ∪ {x n+1 }.

203/224

Approximate KG criterion (AKG)

◮ Let X ref n ⊂ X denote some finite "reference set"

The AKG criterion is:

where the min runs over X ref n ∪ {x n+1 }.

203/224

Approximate KG criterion (AKG)

◮ Let X ref n ⊂ X denote some finite "reference set"

The AKG criterion is:

where the min runs over X ref n ∪ {x n+1 }.

◮ Initially proposed by Scott et al (2011)

◮ under the name KGCP ("KG for continuous parameters")

Approximate KG criterion (AKG)

◮ Let X ref n ⊂ X denote some finite "reference set"

The AKG criterion is:

where the min runs over X ref n ∪ {x n+1 }. ◮ What about simulators with truly non-Gaussian output ?

◮ "batch trick" (CLT) ◮ see also [START_REF] Browne | Stochastic simulators based optimization by gaussian process metamodels-application to maintenance investments planning issues[END_REF]. . .

209/224

Why RBDO is harder than mean-response optimization ◮ because the thresholds p tol j are usually small Bayesian RBDO algorithms ?

◮ A few GP-based algorithms have been proposed, notably: