15 Probabilistic Surrogate Models
The previous chapter discussed how to construct surrogate models from evaluateddesign points. When using surrogate models for the purpose of optimization, it
is often useful to quantify our confidence in the predictions of these models. One
way to quantify our confidence is by taking a probabilistic approach to surrogate
modeling. A common probabilistic surrogate model is the Gaussian process, which
represents a probability distribution over functions. This chapter will explain
how to use Gaussian processes to infer a distribution over the values of different
design points given the values of previously evaluated design points. We will
discuss how to incorporate gradient information as well as noisy measurementsof the objective function. Since the predictions made by a Gaussian process aregoverned by a set of parameters, we will discuss how to infer these parameters
directly from data.

15.1

Gaussian Distribution

Before introducing Gaussian processes, we will first review some relevant properties of the multivariate Gaussian distribution, often also referred to as the multivariate normal distribution.1 An n-dimensional Gaussian distribution is parameterized
by its mean µ and its covariance matrix Σ. The probability density at x is


1
⊤ −1
−n/2
−1/2
(15.1)N (x | µ, Σ) = (2π )|Σ|exp − (x − µ) Σ (x − µ)2
Figure 15.1 shows contour plots of the density functions with different covariancematrices. Covariance matrices are always positive semidefinite.
A value sampled from a Gaussian is writtenx ∼ N (µ, Σ)(15.2)The univariate Gaussian distribution is discussed in appendix C.7.

1

276

c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

"

1
Σ=
0

0
1

#

"

3
Σ=
0

0
0.5

#

"

1
Σ=
0.9

0.9
1

#

x2

2

0

−2
−2

0
x1

2

−2

0
x1

Two jointly Gaussian random variables a and b can be written
#!
" # "
" #
A C
µa
a
,
∼N
C⊤ B
µb
b

−2

2

b ∼ N ( µb , B )(15.3)(15.4)

The conditional distribution for a multivariate Gaussian also has a convenient
closed-form solution:a | b ∼ N µa | b , Σ a | b(15.5)µa|b = µa + CB−1 (b − µb )

Σa|b = A − CB−1 C⊤

2Figure 15.1. Multivariate Gaussians with different covariance matrices.

The marginal distribution2 for a vector of random variables is given by its corresponding mean and covariancea ∼ N ( µa , A )0
x1

2
The marginal distribution is the
distribution of a subset of the variables when the rest are integrated,or marginalized, out. For a distribution over two variables a and b
the marginal distribution over a is:(15.6)(15.7)

Example 15.1 illustrates how to extract the marginal and conditional distributions
from a multivariate Gaussian.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

p( a) =Zp( a, b) db

15 . 2. g a u ssi a n p roc e sse s 277

For example, consider
" # "
" #
3
0
x1
∼N
,
x2
1
1

1
2

#!

Example 15.1. Marginal and conditional distributions for a multivariate Gaussian.

The marginal distribution for x1 is N (0, 3), and the marginal distribution
for x2 is N (1, 2).
The conditional distribution for x1 given x2 = 2 is
µx1 | x2 =2 = 0 + 1 · 2−1 · (2 − 1) = 0.5

Σ x1 | x2 =2 = 3 − 1 · 2−1 · 1 = 2.5
x1 | ( x2 = 2) ∼ N (0.5, 2.5)

15.2

Gaussian Processes

In the previous chapter, we approximated the objective function f using a surrogate model function fˆ fitted to previously evaluated design points. A special
type of surrogate model known as a Gaussian process allows us not only to predict f but also to quantify our uncertainty in that prediction using a probability
distribution.3
A Gaussian process is a distribution over functions. For any finite set of points({x 1) , . . . , x(m) }, the associated function evaluations {y1 , . . . , ym } are distributed
according to:

 

 
k ( x(1) , x(1) ) · · · k ( x(1) , x( m ) )m ( x(1) )
y1

 

 . 
..
..
..
..

,
 .  ∼ N 
(15.8)
.
.
.
.

 

 . 
k ( x( m ) , x(1) ) · · · k ( x( m ) , x( m ) )m (x( m ) )ym

where m(x) is a mean function and k(x, x′ ) is the covariance function, or kernel.4
The mean function can represent prior knowledge about the function. The kernel
controls the smoothness of the functions. Methods for constructing the mean
vector and covariance matrix using mean and covariance functions are given inalgorithm 15.1.

3
A more extensive introduction
to Gaussian processes is provided
by C. E. Rasmussen and C. K. I.
Williams, Gaussian Processes for Machine Learning. MIT Press, 2006.

The mean function produces the
expectation:

4m(x) = E [ f (x)]and the covariance function produces the covariance:k(x, x′ ) =E ( f (x) − m(x))( f (x′ ) − m(x′ ))© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

278

c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

Algorithm 15.1. The function μ forconstructing a mean vector given
a list of design points and a meanfunction m, and the function Σ forconstructing a covariance matrix
given one or two lists of design
points and a covariance functionk.

μ(X, m) = [m(x) for x in X]
Σ(X, k) = [k(x,x′) for x in X, x′ in X]
K(X, X′, k) = [k(x,x′) for x in X, x′ in X′]

A common kernel function is the squared exponential kernel, where( x − x ′ )2k( x, x ) = exp −
2ℓ2
′



(15.9)The parameter ℓ corresponds to what is called the characteristic length-scale, which
can be thought of as the distance we have to travel in design space until theobjective function value changes significantly.5 Hence, larger values of ℓ resultin smoother functions. Figure 15.2 shows functions sampled from a Gaussian
process with a zero-mean function and a squared exponential kernel with different
characteristic length-scales.

ℓ = 1/2

A mathematical definition of characteristic length-scale is provided
by C. E. Rasmussen and C. K. I.
Williams, Gaussian Processes for Machine Learning. MIT Press, 2006.
5

ℓ=1

ℓ=2

y

2

0

−2
0

2

6

4
x

8

10

0

2

6

4

8

10

0

2

x

There are many other kernel functions besides the squared exponential. Several
are shown in figure 15.3. Many kernel functions use r, which is the distance
between x and x′ . Usually the Euclidean distance is used. The Matérn kernel uses
the gamma function Γ, implemented by gamma from the SpecialFunctions.jl
package, and Kν ( x ) is the modified Bessel function of the second kind, implemented
by besselk(ν,x). The neural network kernel augments each design vector with a 1
for ease of notation: x̄ = [1, x1 , x2 , . . .] and x̄′ = [1, x1′ , x2′ , . . .].

6

4

8

10

x
Figure 15.2. Functions sampledfrom Gaussian processes withsquared exponential kernels.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

15 . 2. g a u ssi a n p roc e sse s 279nConstant: σ02

Linear: ∑ σd2 xd xd′

Exponential: exp(− ℓr )γ-Exponential: exp(−( ℓr )γ )Polynomial: ( x ⊤ x ′ + σ02 ) p

d =1

y

2
0

−2

Squared Exponential: exp(− 2rℓ2 )
2

y

2
0

−2

1− ν
Matérn: 2Γ(ν)

√

2ν ℓr

ν

Kν

√

2ν ℓr



Rational Quadratic:2
1 + 2αr ℓ2−αNeural Network: sin

4

−4

−1

q

2x̄⊤ Σx̄′(1+2x̄⊤ Σx̄)(1+2x̄′⊤ Σx̄′ )

!

y

2
0

−2
−4

−2

0
x

2

4

−4

−2

0
x

2

−2

0
x

2

4

Figure 15.3. Functions sampledfrom Gaussian processes with different kernel functions. Shownfunctions are for σ02 = σd2 = ℓ = 1,
p = 2, γ = ν = α = 0.5, and Σ = I.
© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

280c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

This chapter will focus on examples of Gaussian processes with single-dimensional
design spaces for ease of plotting. However, Gaussian processes can be definedover multidimensional design spaces, as illustrated in figure 15.4.

ℓ=1

x1

x1

ℓ=2

x2

ℓ = 1/2

As we will see in section 15.5, Gaussian processes can also incorporate priorindependent noise variance, denoted ν. A Gaussian process is thus defined bymean and covariance functions, prior design points and their function evaluations,and a noise variance. The associated type is given in algorithm 15.2.
mutable struct GaussianProcessm # meank # covariance function
X # design points
y # objective values
ν # noise variance
end

15.3

x1
Figure 15.4. Functions sampledfrom a Gaussian process with zeromean and squared-exponential
kernels over a two-dimensional design space.

Algorithm 15.2. A Gaussian process is defined by a mean function
m, a covariance function k, sampled
design vectors X and their corresponding values y, and a noise variance ν.

Prediction

Gaussian processes are able to represent distributions over functions using conditional probabilities. Suppose we already have a set of points X and the corresponding y, but we wish to predict the values ŷ at points X ∗ . The joint distribution© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

15 . 3. p re d i c ti on

is

" #
"
# "
ŷ
m( X ∗ )K( X ∗ , X ∗ )∼N,ym( X )K( X , X ∗ )K( X ∗ , X )K( X , X )

#!

281

(15.10)In the equation above, we use the functions m and K, which are defined as follows:m( X ) = [m(x(1) ), . . . , m(x(n) )]
k(x(1) , x′(1) ) · · · k(x(1) , x′(m) )


..
..
..

K( X, X ′ ) = 
.
.
.


k(x(n) , x′(1) ) · · · k(x(n) , x′(m) )(15.11)(15.12)The conditional distribution is given byIn the language of Bayesian statistics, the posterior distribution isthe distribution of possible unobserved values conditioned on observed values.

6ŷ | y ∼ N m( X ∗ ) + K( X ∗ , X )K( X, X )−1 (y − m( X )), K( X ∗ , X ∗ ) − K( X ∗ , X )K( X, X )−1 K( X, X ∗ ) (15.13)
{z
}
{z
} |
|
mean

covariance

Note that the covariance does not depend on y. This distribution is often referred
to as the posterior distribution.6 A method for computing and sampling from theposterior distribution defined by a Gaussian process is given in algorithm 15.3.

Algorithm 15.3. The function
mvnrand samples from a multivariate Gaussian with an added inflation factor to prevent numerical issues. The method rand samples a
Gaussian process GP at the given
design points in matrix X.

function mvnrand(μ, Σ, inflation=1e-6)
N = MvNormal(μ, Σ + inflation*I)
return rand(N)endBase.rand(GP, X) = mvnrand(μ(X, GP.m), Σ(X, GP.k))

The predicted mean can be written as a function of x:µ̂(x) = m(x) + K(x, X )K( X, X )−1 (y − m( X ))= m(x) + θ⊤ K( X, x)(15.14)(15.15)where θ = K( X, X )−1 (y − m( X )) can be computed once and reused for different
values of x. Notice the similarity to the surrogate models in the previous chapter. The value of the Gaussian process beyond the surrogate models discussed
previously is that it also quantifies our uncertainty in our predictions.
The variance of the predicted mean can also be obtained as a function of x:ν̂(x) = K(x, x) − K(x, X )K( X, X )−1 K( X, x)(15.16)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

282 c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

In some cases, it is more convenient to formulate equations in terms of thestandard deviation, which is the square root of the variance:qσ̂(x) = ν̂(x)(15.17)The standard deviation has the same units as the mean. From the standard deviation, we can compute the 95% confidence region, which is an interval containing
95% of the probability mass associated with the distribution over y given x. For
a particular x, the 95% confidence region is given by µ̂(x) ± 1.96σ̂ (x). One may
want to use a confidence level different from 95%, but we will use 95% for theplots in this chapter. Figure 15.5 shows a plot of a confidence region associatedwith a Gaussian process fit to four function evaluations.

yconfidence intervaltrue objective function
predicted function mean
fit points

Figure 15.5. A Gaussian processusing the squared exponential kernel and its 95% confidence interval.
Uncertainty increases the farther
we are from a data point. The expected function value approaches
zero as we move far away from thedata point.

x

15.4 Gradient Measurements
Gradient observations can be incorporated into Gaussian processes in a mannerconsistent with the existing Gaussian process machinery.7 The Gaussian process
is extended to include both the function value and its gradient:
"
# "
#!
" #
mf
Kf f Kf∇y∼N,(15.18)m∇
K∇ f K∇∇∇y7
For an overview, see for example
A. O’Hagan, “Some Bayesian Numerical Analysis,” Bayesian Statistics, vol. 4, J. M. Bernardo, J. O.
Berger, A. P. Dawid, and A. F. M.
Smith, eds., pp. 345–363, 1992.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

15 . 4 . g ra d i e n t me a su re me n ts 283


where y ∼ N m f , K f f is a traditional Gaussian process, m∇ is a mean function

for the gradient,8 K f ∇ is the covariance matrix between function values andgradients, K∇ f is the covariance matrix between function gradients and values,
and K∇∇ is the covariance matrix between function gradients.
These covariance matrices are constructed using covariance functions. The
linearity of Gaussians causes these covariance functions to be related:(15.19)k f f (x, x′ ) = k(x, x′ )′(15.20)′k ∇ f (x, x ) = ∇x k(x, x )′(15.21)′k f ∇ (x, x ) = ∇x′ k(x, x )′8Like the mean of the function
value, m∇ is often zero.

′k ∇∇ (x, x ) = ∇x ∇x′ k (x, x )(15.22)

Example 15.2 uses these relations to derive the higher-order covariance functions
for a particular kernel.
Example 15.2. Deriving covariancefunctions for a Gaussian processwith gradient observations.

Consider the squared exponential covariance function


1
′
′ 2
k f f (x, x ) = exp − x − x2
We can use equations (15.19) to (15.22) to obtain the other covariance functions necessary for using Gaussian processes with gradient information:


1
′
′
′ 2
k ∇ f (x, x )i = −(xi − xi ) exp − x − x
2




1
′ 2
′
′
′
k ∇∇ (x, x )ij = − (i = j) − (xi − xi )(x j − x j ) exp − x − x2As a reminder, Boolean expressions, such as (i = j), return 1 if true and 0 if
false.

Prediction can be accomplished in the same manner as with a traditional
Gaussian process. We first construct the joint distribution




 
ŷ
m f (X∗ )K f f (X∗ , X∗ )
 

 
 y  ∼ N  m f ( X )  ,  K f f ( X , X ∗ )∇ym∇ ( X )K∇ f ( X , X ∗ )K f f (X∗ , X )K f f (X , X )K∇ f ( X , X )
K f ∇ (X∗ , X )
K f ∇ ( X , X ) K∇∇ ( X , X )(15.23)© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

284 c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

For a Gaussian process over n-dimensional design vectors given m pairs offunction and gradient evaluations and ℓ query points, the covariance blocks havethe following dimensions:ℓ×ℓℓ×mℓ × nmm×ℓ
m×m
m × nmnm × ℓ nm × m nm × nm

(15.24)

Example 15.3 constructs such a covariance matrix.
Suppose we have evaluated a function and its gradient at two locations, x(1)and x(2) , and we wish to predict the function value at x̂. We can infer thejoint distribution over ŷ, y, and ∇y using a Gaussian process. The covariancematrix is:

k f f (x̂, x̂) k (x(1) , x̂) ff
 k (x(2) , x̂)
 ff

(1)
k
 ∇ f (x , x̂)1
(1)
k
 ∇ f (x , x̂)2(2) , x̂)
k
(
x
1
 ∇f
k ∇ f (x(2) , x̂)2

k f f (x̂, x(1) )k f f (x(1) , x(1) )k f f (x(2) , x(1) )k ∇ f (x(1) , x(1) )1
k ∇ f (x(1) , x(1) )2
k ∇ f (x(2) , x(1) )1
k ∇ f (x(2) , x(1) )2

k f f (x̂, x(2) )k f f (x(1) , x(2) )k f f (x(2) , x(2) )k ∇ f (x(1) , x(2) )1
k ∇ f (x(1) , x(2) )2k ∇ f (x(2) , x(2) )1
k ∇ f (x(2) , x(2) )2

k f ∇ (x̂, x(1) )1
k f ∇ (x(1) , x(1) )1
k f ∇ (x(2) , x(1) )1k ∇∇ (x(1) , x(1) )11
k ∇∇ (x(1) , x(1) )21
k ∇∇ (x(2) , x(1) )11
k ∇∇ (x(2) , x(1) )21

k f ∇ (x̂, x(1) )2
k f ∇ (x(1) , x(1) )2
k f ∇ (x(2) , x(1) )2k ∇∇ (x(1) , x(1) )12
k ∇∇ (x(1) , x(1) )22
k ∇∇ (x(2) , x(1) )12
k ∇∇ (x(2) , x(1) )22

k f ∇ (x̂, x(2) )1
k f ∇ (x(1) , x(2) )1
k f ∇ (x(2) , x(2) )1
k ∇∇ (x(1) , x(2) )11
k ∇∇ (x(1) , x(2) )21
k ∇∇ (x(2) , x(2) )11
k ∇∇ (x(2) , x(2) )21

Example 15.3. Constructing the covariance matrix for a Gaussian process with gradient observations.

k f ∇ (x̂, x(2) )2
k f ∇ (x(1) , x(2) )2 
k f ∇ (x(2) , x(2) )2 


(
1
)
(
2
)
k ∇∇ (x , x )12 


(
1
)
(
2
)
k ∇∇ (x , x )22 k ∇∇ (x(2) , x(2) )12 k ∇∇ (x(2) , x(2) )22

The conditional distribution follows the same Gaussian relations as in equation (15.13):(15.25)ŷ | y, ∇y ∼ N (µ∇ , Σ∇ )where:## −1 "#⊤ "y − m f (X)K f f (X , X ) K f ∇ (X , X )K f f (X , X∗ )µ∇ = m f ( X ) +∇ y − m∇ ( X )K∇ f ( X , X ) K∇∇ ( X , X )K∇ f ( X , X ∗ )"#⊤ "# −1 "#
K f f (X , X∗ )K f f (X , X ) K f ∇ (X , X )K f f (X , X∗ )∗∗
Σ∇ = K f f ( X , X ) −K∇ f ( X , X ∗ )K∇ f ( X , X ) K∇∇ ( X , X )K∇ f ( X , X ∗ )∗"The regions obtained when including gradient observations are compared tothose without gradient observations in figure 15.6.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com(15.26)(15.27)

15 . 5 . n oi sy me a su re me n ts 285

yconfidence interval without gradientconfidence interval with gradientpredicted mean without gradientpredicted mean with gradientfit points
true function

Figure 15.6. Gaussian processes
with and without gradient information using squared exponential
kernels. Incorporating gradient information can significantly reducethe confidence intervals.

x

15.5 Noisy Measurements
So far we have assumed that the objective function f is deterministic. In practice,
however, evaluations of f may include measurement noise, experimental error,or numerical roundoff.
We can model noisy evaluations as y = f (x) + z, where f is deterministic but
z is zero-mean Gaussian noise, z ∼ N (0, ν). The variance of the noise ν can be
adjusted to control the uncertainty.9
The new joint distribution is:
#!
# "
"
" #
K( X ∗ , X ∗ ) K( X ∗ , X )m( X ∗ )ŷ(15.28),∼N
K( X , X ∗ ) K( X , X ) + νIm( X )y

The techniques covered in section 14.5 can be used to tune thevariance of the noise.
9with conditional distribution:ŷ | y, ν ∼ N (µ∗ , Σ∗ )∗∗(15.29)∗µ = m( X ) + K( X , X )(K( X, X ) + νI)−1(y − m( X ))Σ∗ = K( X ∗ , X ∗ ) − K( X ∗ , X )(K( X, X ) + νI)−1 K( X, X ∗ )(15.30)(15.31)As the equations above show, accounting for Gaussian noise is straightforward
and the posterior distribution can be computed analytically. Figure 15.7 showsa noisy Gaussian process. Algorithm 15.4 implements prediction for Gaussian
processes with noisy measurements.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

286 c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

Figure 15.7. A noisy Gaussian process using a squared exponential
kernel.

yconfidence regiontrue objective function
predicted function mean
fit pointsxfunction predict(GP, X_pred)m, k, ν = GP.m, GP.k, GP.νtmp = K(X_pred, GP.X, k) / (K(GP.X, GP.X, k) + ν*I)
μₚ = μ(X_pred, m) + tmp*(GP.y - μ(GP.X, m))
S = K(X_pred, X_pred, k) - tmp*K(GP.X, X_pred, k)
νₚ = diag(S) .+ eps() # eps prevents numerical issues
return (μₚ, νₚ)
end

Algorithm 15.4. A method for obtaining the predicted means andstandard deviations in f undera Gaussian process. The method
takes a Gaussian process GP anda list of points X_pred at which toevaluate the prediction. It returns
the mean and variance at each evaluation point.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

15 . 6 . fi tti n g g a u ssi a n p roc e sse s 287

15.6

Fitting Gaussian Processes

The choice of kernel and parameters has a large effect on the form of the Gaussian
process between evaluated design points. Kernels and their parameters can bechosen using cross validation introduced in the previous chapter. Instead ofminimizing the squared error on the test data, we maximize the likelihood of thedata.10 That is, we seek the parameters θ that maximize the probability of thefunction values, p(y | X, θ). The likelihood of the data is the probability that theobserved points were drawn from the model. Equivalently, we can maximize thelog likelihood, which is generally preferable because multiplying small probabilities
in the likelihood calculation can produce extremely small values. Given a dataset
D with n entries, the log likelihood is given by

Alternatively, we could maximize the pseudolikelihood as discussed by C. E. Rasmussen and
C. K. I. Williams, Gaussian Processes
for Machine Learning. MIT Press,
2006.

10

n
1
1
log 2π − log |Kθ ( X, X ) + νI| − (y − mθ ( X ))⊤ (Kθ ( X, X ) + νI)−1 (y − mθ ( X ))2
2
2
(15.32)where the mean and covariance functions are parameterized by θ.
Let us assume a zero mean such that mθ ( X ) = 0 and θ refers only to theparameters for the Gaussian process covariance function. We can arrive at a
maximum likelihood estimate by gradient ascent. The gradient is then given by!
1 ⊤ −1 ∂K −1
1
∂
−1 ∂K
(15.33)log p(y | X, θ) = y K
K y − tr Σθ
∂θ j
2
∂θ j
2
∂θ j
log p(y | X, ν, θ) = −where Σθ = Kθ ( X, X ) + νI. Above, we use the matrix derivative relations
∂K −1
∂K−1
= − K−1
K
∂θ j
∂θ j
!
∂ log |K|
−1 ∂K
= tr K∂θ j∂θ j(15.34)(15.35)where tr(A) denotes the trace of a matrix A, defined to be the sum of the elementson the main diagonal.

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

288

c ha p te r 15 . p roba bi l i sti c su rrog a te mod e l s

15.7 Summary• Gaussian processes are probability distributions over functions.
• The choice of kernel affects the smoothness of the functions sampled from a
Gaussian process.
• The multivariate normal distribution has analytic conditional and marginal
distributions.
• We can compute the mean and standard deviation of our prediction of an
objective function at a particular design point given a set of past evaluations.
• We can incorporate gradient observations to improve our predictions of theobjective value and its gradient.
• We can incorporate measurement noise into a Gaussian process.
• We can fit the parameters of a Gaussian process using maximum likelihood.

15.8 Exercises
Exercise 15.1. Gaussian processes will grow in complexity during the optimization process as more samples accumulate. How can this be an advantage overmodels based on regression?
Exercise 15.2. How does the computational complexity of prediction with aGaussian process increase with the number of data points m?

Exercise 15.3. Consider the function f ( x ) = sin( x )/ x2 + 1 over [−5, 5]. Plotthe 95% confidence bounds for a Gaussian process with derivative information
fitted to the evaluations at {−5, −2.5, 0, 2.5, 5}. What is the maximum standard
deviation in the predicted distribution within [−5, 5]? How many function evaluations, evenly-spaced over the domain, are needed such that a Gaussian process
without derivative information achieves the same maximum predictive standard
deviation?

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com

15 . 8 . e x e rc i se s 289

Assume zero-mean functions and noise-free observations, and use the covariance functions:1
2
k f f ( x, x ′ ) = exp − x − x ′ 2
2


1
2
k ∇ f ( x, x ′ ) = ( x ′ − x ) exp − x − x ′ 2212
k ∇∇ ( x, x ′ ) = (( x − x ′ )2 − 1) exp(− x − x ′ 2 )2
Exercise 15.4. Derive the relation k f ∇ (x, x′ )i = cov f (x), ∂x∂ ′ f (x′ ) = ∂x∂ ′ k f f (x, x′ ).
iiExercise 15.5. Suppose we have a multivariate Gaussian distribution over two
variables a and b. Show that the variance of the conditional distribution over agiven b is no greater than the variance of the marginal distribution over a. Does
this make intuitive sense?

Exercise 15.6. Suppose we observe many outliers while sampling, that is, weobserve samples that do not fall within the confidence interval given by theGaussian process. This means the probabilistic model we chose is not appropriate.
What can we do?
Exercise 15.7. Consider model selection for the function evaluation pairs ( x, y):{(1, 0), (2, −1), (3, −2), (4, 1), (5, 0)}Use leave-one-out cross-validation to select the kernel that maximizes thelikelihood of predicting the withheld pair given a Gaussian process over the otherpairs in the fold. Assume zero mean with no noise. Select from the kernels:exp(−k x − x ′ k)exp(−k x − x ′ k2 )(1 + k x − x ′ k)−1

(1 + k x − x ′ k2 ) −1

(1 + k x − x ′ k)−2

© 2019 Massachusetts Institute of Technology, shared under a Creative Commons CC-BY-NC-ND license.
2022-05-22 00:25:57-07:00, revision 47fd495, comments to bugs@algorithmsbook.com