14

Surrogate Models

The previous chapter discussed methods for producing a sampling plan. This
chapter shows how to use these samples to construct models of the objectivefunction that can be used in place of the real objective function. Such surrogate
models are designed to be smooth and inexpensive to evaluate so that they can be
efficiently optimized. The surrogate model can then be used to help direct the
search for the optimum of the real objective function.

Fitting Surrogate Models

y

14.1

A surrogate model fˆ parameterized by θ is designed to mimic the true objectivefunction f . The parameters θ can be adjusted to fit the model based on samples
collected from f . An example surrogate model is shown in figure 14.1.
Suppose we have m design pointsnoX = x(1) , x(2) , . . . , x( m )(14.1)
and associated function evaluationsnoy = y (1) , y (2) , . . . , y ( m )For a particular set of parameters, the model will predictonŷ = fˆθ (x(1) ), fˆθ (x(2) ), . . . , fˆθ (x(m) )(14.2)

θky − ŷk p

design points
surrogate model
true objective function
Figure 14.1. Surrogate models approximate the true objective function. The model is fitted to the evaluated design points but deviates
farther away from them.

(14.3)

Fitting a model to a set of points requires tuning the parameters to minimize
the difference between the true evaluations and those predicted by the model,
typically according to an L p norm:1
minimizex(14.4)

It is common to use the L2 norm.
Minimizing this equation with an
L2 norm is equivalent to minimizing the mean squared error at thosedata points.

1

254

c ha p te r 14 . su rrog a te mod e l s

Equation (14.4) penalizes the deviation of the model only at the data points.
There is no guarantee that the model will continue to fit well away from observed
data, and model accuracy typically decreases the farther we go from the sampled
points.
This form of model fitting is called regression. A large body of work exists forsolving regression problems, and it is extensively studied in machine learning.2
The rest of this chapter covers several popular surrogate models and algorithms
for fitting surrogate models to data, and concludes with methods for choosing
between types of models.

2
K. P. Murphy, Machine Learning: A
Probabilistic Perspective. MIT Press,
2012.

14.2 Linear Models
A simple surrogate model is the linear model, which has the form3fˆ = w0 + w⊤ x

θ = { w0 , w }This equation may seem familiar.
It is the equation for a hyperplane.

3(14.5)For an n-dimensional design space, the linear model has n + 1 parameters, andthus requires at least n + 1 samples to fit unambiguously.
Instead of having both w and w0 as parameters, it is common to construct asingle vector of parameters θ = [w0 , w] and prepend 1 to the vector x to getfˆ = θ⊤ x(14.6)

Finding an optimal θ requires solving a linear regression problem:
minimizeθky − ŷk22(14.7)

ky − Xθk22(14.8)

which is equivalent to solving
minimizeθ

where X is a design matrix formed from m data points( x(1) ) ⊤ (2) ⊤ 
 (x ) 

X=
..




.

(14.9)( x( m ) ) ⊤© 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

14 . 3. ba si s fu n c ti on s 255

Algorithm 14.1. A method for constructing a design matrix from a listof design points X and a method
for fitting a surrogate model using
linear regression to a list of design
points X and a vector of objectivefunction values y.

function design_matrix(X)n, m = length(X[1]), length(X)return [j==0 ? 1.0 : X[i][j] for i in 1:m, j in 0:n]end
function linear_regression(X, y)θ = pinv(design_matrix(X))*yreturn x -> θ⋅[1; x]
end

Algorithm 14.1 implements methods for computing a design matrix and forsolving a linear regression problem. Several cases for linear regression are shownin figure 14.2.
Linear regression has an analytic solution(14.10)θ = X+ ywhere X+ is the Moore-Penrose pseudoinverse of X.
If X⊤ X is invertible, the pseudoinverse can be computed as

 −1
X⊤
X+ = X⊤ X(14.11) −1
X+ = X⊤ XX⊤(14.12)

If XX⊤ is invertible, the pseudoinverse can be computed as

The function pinv computes the pseudoinverse of a given matrix.4

14.3

The function pinv uses the singular value decomposition, X =UΣV∗ , to compute the pseudoinverse:
4Basis Functions

X+ = VΣ+ U∗

The linear model is a linear combination of the components of x:nfˆ(x) = θ1 x1 + · · · + θn xn = ∑ θi xi = θ⊤ x

(14.13)i =1

which is a specific example of a more general linear combination of basis functionsqfˆ(x) = θ1 b1 (x) + · · · + θq bq (x) = ∑ θi bi (x) = θ⊤ b(x)where the pseudoinverse of the diagonal matrix Σ is obtained by taking the reciprocal of each nonzero
element of the diagonal and thentransposing the result.

(14.14)i =1© 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

256 c ha p te r 14 . su rrog a te mod e l s

m = n+1yym < n+1x

nonindependent points

m > n+1yyxFigure 14.2. Models resulting fromlinear regression, which minimizes
the square vertical distance of themodel from each point. The pseudoinverse produces a unique solution for any nonempty point configuration.
The bottom-left subfigure
shows the model obtained for
two repeated points, in this case,m = n + 1. Because the two entries
are repeated, the matrix X is
singular. Although X does nothave an inverse in this case, thepseudoinverse produces a uniquesolution that passes between thetwo points.

xx© 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

14 . 3. ba si s fu n c ti on s

257

In the case of linear regression, the basis functions simply extract each component,bi ( x ) = x i .
Any surrogate model represented as a linear combination of basis functions
can be fit using regression:
minimizeθky − Bθk22

where B is the basis matrix formed from m data points:b ( x(1) ) ⊤ b ( x(2) ) ⊤ 


B=
..



.
(
m
)
⊤
b(x )

(14.15)(14.16)

The weighting parameters can be obtained using the pseudoinverse
θ = B+ y(14.17)

Algorithm 14.2 implements this more general regression procedure.
using LinearAlgebra
function regression(X, y, bases)
B = [b(x) for x in X, b in bases]θ = pinv(B)*yreturn x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))
end

Algorithm 14.2. A method for fitting a surrogate model to a list ofdesign points X and corresponding
objective function values y using regression with basis functions contained in the bases array.

Linear models cannot capture nonlinear relations. There are a variety of otherfamilies of basis functions that can represent more expressive surrogate models.
The remainder of this section discusses a few common families.

14.3.1 Polynomial Basis Functions
Polynomial basis functions consist of a product of design vector components, eachraised to a power. Linear basis functions are a special case of polynomial basis
functions.
From the Taylor series expansion5 we know that any infinitely differentiable
function can be closely approximated by a polynomial of sufficient degree. We
can construct these bases using algorithm 14.3.

5

Covered in appendix C.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

258 c ha p te r 14 . su rrog a te mod e l s

In one dimension, a polynomial model of degree k has the formkfˆ( x ) = θ0 + θ1 x + θ2 x2 + θ3 x3 + · · · + θk x k = ∑ θi xi(14.18)i =0Hence, we have a set of basis functions bi ( x ) = xi for i ranging from 0 to k.
In two dimensions, a polynomial model of degree k has basis functions of theformjbij (x) = x1i x2 for i, j ∈ {0, . . . , k}, i + j ≤ k(14.19)Fitting a polynomial surrogate model is a regression problem, so a polynomial
model is linear in higher dimensional space (figure 14.3). Any linear combinationof basis functions can be viewed as linear regression in a higher dimensional
space.

yyf (x)fˆ( x )xx

Figure 14.3. A polynomial model
is linear in higher dimensions. The
function exists in the plane formed
from its bases, but it does not occupy the entire plane because theterms are not independent.

x2

polynomial_bases_1d(i, k) = [x->x[i]^p for p in 0:k]function polynomial_bases(n, k)bases = [polynomial_bases_1d(i, k) for i in 1 : n]terms = Function[]for ks in Iterators.product([0:k for i in 1:n]...)
if sum(ks) ≤ k
push!(terms,
x->prod(b[j+1](x) for (j,b) in zip(ks,bases)))end
end
return terms
end

Algorithm 14.3. A method for constructing an array of polynomial
basis functions up to a degree k
for the ith component of a design
point, and a method for constructing a list of n-dimensional polynomial bases for terms up to degreek.

© 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

14 . 3. ba si s fu n c ti on s 259

14.3.2 Sinusoidal Basis Functions
Any continuous function over a finite domain can be represented using an infinite set of sinusoidal basis functions.6 A Fourier series can be constructed for any
integrable univariate function f on an interval [ a, b]
f (x) =





∞
2πix
2πix
θ0
(sin)
(cos)
sin
cos
+ θi
+ ∑ θi
2
b−a
b−a
i =1

The Fourier series is also exact forfunctions defined over the entire
real line if the function is periodic.

6(14.20)

where
Z

b
2
f ( x ) dx
b−a a


Z b
2
2πix
(sin)
dx
f ( x ) sin
=
θi
b−a a
b−a


Z b
2
2πix
(cos)
θi
=
f ( x ) cos
dx
b−a a
b−a

θ0

=

(14.21)(14.22)(14.23)

Just as the first few terms of a Taylor series are used in polynomial models, so
too are the first few terms of the Fourier series used in sinusoidal models. The
bases for a single component over the domain x ∈ [ a, b] are:



b (x)
= 1/2



 0
(sin)
bi ( x ) = sin 2πix
(14.24)
b
−
a




(
cos
)

2πix
b
( x ) = cos b−ai

We can combine the terms for multidimensional sinusoidal models in the same
way we combine terms in polynomial models. Algorithm 14.4 can be used toconstruct sinusoidal basis functions. Several cases for sinusoidal regression are
shown in figure 14.4.

14.3.3 Radial Basis Functions
A radial function ψ is one which depends only on the distance of a point fromsome center point c, such that it can be written ψ(x, c) = ψ(kx − ck) = ψ(r ).
Radial functions are convenient basis functions because placing a radial function
contributes a hill or valley to the function landscape. Some common radial basis
functions are shown in figure 14.5.
© 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

260

c ha p te r 14 . su rrog a te mod e l sfunction sinusoidal_bases_1d(j, k, a, b)
T = b[j] - a[j]bases = Function[x->1/2]
for i in 1 : k
push!(bases, x->sin(2π*i*x[j]/T))
push!(bases, x->cos(2π*i*x[j]/T))
end
return bases
endfunction sinusoidal_bases(k, a, b)
n = length(a)bases = [sinusoidal_bases_1d(i, k, a, b) for i in 1 : n]terms = Function[]for ks in Iterators.product([0:2k for i in 1:n]...)
powers = [div(k+1,2) for k in ks]
if sum(powers) ≤ k
push!(terms,
x->prod(b[j+1](x) for (j,b) in zip(ks,bases)))end
end
return terms
end

Algorithm 14.4.

The method

sinusoidal_bases_1d produces alist of basis functions up to degree
k for the ith component of the design vector given lower bound a
and upper bound b. The method
sinusoidal_bases produces allbase function combinations up todegree k for lower-bound vector aand upper-bound vector b.

© 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

14 . 3. ba si s fu n c ti on s

m = n+1

Figure 14.4. Fitting sinusoidal
models to noisy points.

yym < n+1

261

x

nonindependent points

m > n+1

y

y

x

x

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

262 c ha p te r 14 . su rrog a te mod e l sthin plate spline: r2 log r

cubic: r3

linear: r

3

8

2

ψ

6
1

2

4

1

2
0

0

0

Gaussian: e−r /2σ22

multiquadratic: (r2 + σ2 ) 2

1

inverse multiquadratic: (r2 + σ2 )− 2
1

1

1
2
ψ

0.8
0.5

1.5
0.6
1

0
0

1
r

2

0

1
r

2

0.4

0

1
r

2

Figure 14.5. Several radial basis
functions.

© 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

14 . 4 . fi tti n g n oi sy obj e c ti ve fu n c ti on s 263

Radial basis functions require specifying the center points. One approach when
fitting radial basis functions to a set of data points is to use the data points as thecenters. For a set of m points, one thus constructs m radial basis functionsbi (x) = ψ(kx − x(i) k)

for i ∈ {1, . . . , m}

(14.25)The corresponding m × m basis matrix is always semidefinite. Algorithm 14.5 canbe used to construct radial basis functions with known center points. Surrogate
models with different radial basis functions are shown in figure 14.6.
radial_bases(ψ, C, p=2) = [x->ψ(norm(x - c, p)) for c in C]

0.5

x sin(5x )ψ = exp(−2r2 )ψ = exp(−5r2 )y0

Algorithm 14.5. A method for obtaining a list of basis functions
given a radial basis function ψ, alist of centers C, and an L p norm
parameter p.
Figure 14.6. Several different Gaussian radial basis functions used tofit x sin(5x ) based on 4 noise-freesamples.

ψ = exp(−10r2 )

−0.5

−1

−1.5

0

0.2

0.6

0.4

0.8

1

x

14.4

Fitting Noisy Objective Functions

Models fit using regression will pass as close as possible to every design point.
When the objective function evaluations are noisy, complex models are likely to
excessively contort themselves to pass through every point. However, smootherfits are often better predictors of the true underlying objective function.

© 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

264

c ha p te r 14 . su rrog a te mod e l s

The basis regression problem specified in equation (14.15) can be augmented
to prefer smoother solutions. A regularization term is added in addition to theprediction error in order to give preference to solutions with lower weights. The
resulting basis regression problem with L2 regularization7 is:minimizeθky − Bθk22 + λkθk22(14.26)

where λ ≥ 0 is a smoothing parameter, with λ = 0 resulting in no smoothing.
The optimal parameter vector is given by:8

 −1
B⊤ y
θ = B⊤ B + λI

(14.27)

where I is the identity matrix.
Algorithm 14.6 implements regression with L2 regularization. Surrogate models with different radial basis functions fit to noisy samples are shown in figure 14.7.
function regression(X, y, bases, λ)B = [b(x) for x in X, b in bases]θ = (B'B + λ*I)\B'y
return x -> sum(θ[i] * bases[i](x) for i in 1 : length(θ))
end

Other L p -norms, covered in appendix C.4, can be used as well.
Using the L1 norm will encouragesparse solutions with less influential component weights set to zero,
which can be useful in identifying
important basis functions.

8
The matrix B⊤ B + λI is not always invertible if λ = 0. However,
we can always produce an invertible matrix with a positive λ.
7

Algorithm 14.6. A method for regression in the presence of noise,
where λ is a smoothing term. It returns a surrogate model fitted toa list of design points X and corresponding objective function values
y using regression with basis functions bases.

14.5 Model Selection
So far, we have discussed how to fit a particular model to data. This section
explains how to select which model to use. We generally want to minimize generalization error, which is a measure of the error of the model on the full design
space, including points that may not be included in the data used to train themodel. One way to measure generalization error is to use the expected squared
error of its predictions:2 ǫgen = Ex∼X f (x) − fˆ(x)(14.28)© 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

14 . 5 . mod e l se l e c ti on

2yFigure 14.7. Several different Gaussian radial basis functions used tofit x sin(5x ) with zero mean, standard deviation 0.1 error based on
10 noisy samples and radial basis
function ψ = exp(−5r2 ).

x sin(5x )λ=0
λ = 0.1
λ = 0.5

1

265

0

−1

−2

0

0.2

0.6

0.4

0.8

1

x

Of course, we cannot calculate this generalization error exactly because it requires
knowing the function we are trying to approximate. It may be tempting to estimate
the generalization error of a model from the training error. One way to measuretraining error is to use the mean squared error (MSE) of the model evaluated on
the m samples used for training:ǫtrain =2
1 m
f (x(i) ) − fˆ(x(i) )∑m i =1(14.29)However, performing well on the training data does not necessarily correspondto low generalization error. Complex models may reduce the error on the training
set, but they may not provide good predictions in other points in the design space
as illustrated in example 14.1.9
This section discusses several methods for estimating generalization error.
These methods train and test on subsets of the data. We introduce the TrainTest
type (algorithm 14.7), which contains a list of training indices and a list of test
indices. The method fit takes in a training set and produces a model. The method
metric takes a model and a test set and produces a metric, such as the meansquared error. The method train_and_validate (algorithm 14.7) is a utility
function for training and then evaluating a model. Although we train on subsetsof the data when estimating the generalization error, once we have decided whichmodel to use, we can train on the full dataset.

A major theme in machine learning is balancing model complexity to avoid overfitting the training
data. K. P. Murphy, Machine Learning: A Probabilistic Perspective. MIT
Press, 2012.
9© 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

266 c ha p te r 14 . su rrog a te mod e l s

Consider fitting polynomials of varying degrees to evaluations of the objective functionf ( x ) = x/10 + sin( x )/4 + exp − x2

Example 14.1. A comparison oftraining and generalization error
as the degree of a polynomial surrogate model is varied.

Below we plot polynomial surrogate models of varying degrees using thesame nine evaluations evenly spaced over [−4, 4]. The training and generalization error are shown as well, where generalization is calculated over[−5, 5].

y

1

k=0

k=1

k=2

0
ǫtrain = 0.199
ǫgen = 0.187

ǫtrain = 0.089
ǫgen = 0.1

ǫtrain = 0.132ǫgen = 0.127

−1

y

1

k=3

k=5

k=4

0
ǫtrain = 0.062ǫgen = 0.09

ǫtrain = 0.029
ǫgen = 0.129ǫtrain = 0.027
ǫgen = 0.112

−1

y

1

k=6

k=8

k=7

0
ǫtrain = 0.008
ǫgen = 1.162

−1

−4 −2

0
x

2

4

ǫtrain = 0.008
ǫgen = 1.157

−4 −2

0
x

2

4

ǫtrain = 0.0ǫgen = 87.633

−4 −2

0
x

2

4

The plot shows that the generalization error is high for both very low andhigh values of k, and that training error decreases as we increase the polynomial degree. The high-degree polynomials are particularly poor predictors
for designs outside [−4, 4].

© 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

14 . 5 . mod e l se l e c ti on 267

struct TrainTest
train
test
end
function train_and_validate(X, y, tt, fit, metric)model = fit(X[tt.train], y[tt.train])return metric(model, X[tt.test], y[tt.test])
end

14.5.1 Holdout
test

train

Figure 14.8. The holdout method(left) partitions the data into train
and test sets.

generalization error estimate

test( fˆ, •)A simple approach to estimating the generalization error is the holdout method,which partitions the available data into a test set Dh with h samples and a training
set Dt consisting of all remaining m − h samples as shown in figure 14.8. The
training set is used to fit model parameters. The held out test set is not usedduring model fitting, and can thus be used to estimate the generalization error.
Different split ratios are used, typically ranging from 50% train, 50% test to 90%
train, 10% test, depending on the size and nature of the dataset. Using too few
samples for training can result in poor fits (figure 14.9), whereas using too many
will result in poor generalization estimates.
The holdout error for a model fˆ fit to the training set is2
1(14.30)ǫholdout =y − fˆ(x)∑h (x,y)∈D

ytrain(•)

Algorithm 14.7. A utility type andmethod for training a model andthen validating it on a metric. Here,train and test are lists of indicesinto the training data, X is a list ofdesign points, y is the vector of corresponding function evaluations,tt is a train-test partition, fit is amodel fitting function, and metricevaluates a model on a test set to
produce an estimate of generalization error.

xffˆ
training samples
holdout samples
Figure 14.9. Poor train-test splits
can result in poor model performance.

hEven if the partition ratio is fixed, the holdout error will depend on the particular train-test partition chosen. Choosing a partition at random (algorithm 14.8)
will only give a point estimate. In random subsampling (algorithm 14.9), we apply
the holdout method multiple times with randomly selected train-test partitions.
The estimated generalization error is the mean over all runs.10 Because the validation sets are chosen randomly, this method does not guarantee that we validateon all of the data points.

The standard deviation over all
runs can be used to estimate thestandard deviation of the estimated generalization error.

10© 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

268 c ha p te r 14 . su rrog a te mod e l s

Algorithm 14.8. A method for
randomly partitioning m data samples into training and holdout sets,
where h samples are assigned tothe holdout set.

function holdout_partition(m, h=div(m,2))p = randperm(m)train = p[(h+1):m]
holdout = p[1:h]return TrainTest(train, holdout)
endfunction random_subsampling(X, y, fit, metric;
h=div(length(X),2), k_max=10)m = length(X)mean(train_and_validate(X, y, holdout_partition(m, h),fit, metric) for k in 1 : k_max)
end

Algorithm 14.9. The random subsampling method used to obtainmean and standard deviation estimates for model generalization error using k_max runs of the holdout
method.

14.5.2 Cross Validation

D1

D2

train(••••)
train(••••)
train(••••)
train(••••)
train(••••)

D3
test( fˆ, •)test( fˆ, •)test( fˆ, •)test( fˆ, •)test( fˆ, •)D4

D5

generalization error estimate
generalization error estimate
generalization error estimate
generalization error estimate
generalization error estimate

Figure 14.10. Cross-validation partitions the data into equally sized
sets. Each set is the holdout set
once. Here we show 5-fold crossvalidation.

generalization error µ and σ
Using a train-test partition can be wasteful because our model tuning can takeadvantage only of a segment of our data. Better results are often obtained by using
k-fold cross validation.11 Here, the original dataset D is randomly partitioned into ksets D1 , D2 , . . . , Dk of equal, or approximately equal, size, as shown in figure 14.10
and implemented in algorithm 14.10. We then train k models, one on each subsetof k − 1 sets, and we use the withheld set to estimate the generalization error. The
cross-validation estimate of generalization error is the mean generalization error
over all folds:12

11

Also known as rotation estimation.

As with random subsampling,an estimate of variance can be obtained from the standard deviation
over folds.
© 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
1214 . 5 . mod e l se l e c ti on

1 k (i )ǫcross-validation
k i∑
=1

2
1
(i )
ˆ(i) (x)y−fǫcross-validation =∑(i )|Dtest | (x,y)∈D (i)269ǫcross-validation =(14.31)

test

(i )
(i )
where ǫcross-validation and Dtest are the cross-validation estimate and the withheld

test set, respectively, for the ith fold.

function k_fold_cross_validation_sets(m, k)
perm = randperm(m)sets = TrainTest[]
for i = 1:kvalidate = perm[i:k:m];train = perm[setdiff(1:m, i:k:m)]
push!(sets, TrainTest(train, validate))end
return sets
endfunction cross_validation_estimate(X, y, sets, fit, metric)mean(train_and_validate(X, y, tt, fit, metric)
for tt in sets)
end

Cross-validation also depends on the particular data partition. An exception isleave-one-out cross-validation with k = m, which has a deterministic partition. It
trains on as much data as possible, but it requires training m models.13 Averagming over all (m/k) possible partitions, known as complete cross-validation, is often
too expensive. While one can average multiple cross-validation runs, it is morecommon to average the models from a single cross-validation partition.
Cross-validation is demonstrated in example 14.2.

Algorithm 14.10.

The methodk_fold_cross_validation_sets

constructs the sets needed for
k-fold cross validation on m
samples, with k ≤ m. The method
cross_validation_estimate

computes the mean of the generalization error estimate bytraining and validating on the listof train-validate sets contained
in sets. The other variables are
the list of design points X, thecorresponding objective function
values y, a function fit that trains
a surrogate model, and a function
metric that evaluates a model ona data set.

M. Stone, “Cross-Validatory
Choice and Assessment of Statistical Predictions,” Journal of the Royal
Statistical Society, vol. 36, no. 2,
pp. 111–147, 1974.
13© 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

270 c ha p te r 14 . su rrog a te mod e l s

14.5.3 The Bootstrap
The bootstrap method14 uses multiple bootstrap samples, which consist of m indicesinto a dataset of size m independently chosen uniformly at random. The indices
are chosen with replacement, so some indices may be chosen multiple times andsome indices may not be chosen at all as shown in figure 14.11. The bootstrap
sample is used to fit a model that is then evaluated on the original training set. A
method for obtaining bootstrap samples is given in algorithm 14.11.
If b bootstrap samples are made, then the bootstrap estimate of the generalization error is the mean of the corresponding generalization error estimates

14
B. Efron, “Bootstrap Methods:
Another Look at the Jackknife,” The
Annals of Statistics, vol. 7, pp. 1–26,
1979.

ǫtest , . . . , ǫtest :
(1)

(b)

ǫboot ==

1 b (i )
ǫtest
b i∑
=1

(14.32)

2
1 m 1 b  ( j)(i ) ( j )
ˆ
y
−
f
(
x
)
m j∑
b∑
=1 i =1(14.33)where fˆ(i) is the model fit to the ith bootstrap sample. The bootstrap method isimplemented in algorithm 14.12.
The bootstrap error in equation (14.32) tests models on data points to which
they were fit. The leave-one-out bootstrap estimate removes this source of bias byonly evaluating fitted models to withheld data:2
1 m 1 b  y( j) − fˆ(i) (x( j) )ǫleave-one-out-boot =
m j∑
c ∑
=1 − j i =1 0

if jth index was not in the ith bootstrap sample
otherwise

where c− j is the number of bootstrap samples that do not contain index j. The
leave-one-out bootstrap method is implemented in algorithm 14.13.
The probability of a particular index not being in a bootstrap sample is:1 m≈ e−1 ≈ 0.368(14.35)1−m
so a bootstrap sample is expected to have on average 0.632m distinct indices fromthe original dataset.

© 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(14.34)

14 . 5 . mod e l se l e c ti on 271

f ( x ) = sin(2x ) cos(10x ) + ǫ/10
where x ∈ [0, 1] and ǫ is random noise with zero mean and unit variance,ǫ ∼ N (0, 1).

Example 14.2.
Cross validation used to fit a hyperparameter.
mean cross validated MSE

Suppose we want to fit a noisy objective function using radial basis functions
with the noise hyperparameter λ (section 14.4). We can use cross validation
to determine λ. We are given ten samples from our noisy objective function.
In practice, the objective function will be unknown, but this example uses

0.1
5 · 10−2
0

10−4

10−1
λ

102

Random.seed!(0)f = x->sin(2x)*cos(10x)X = rand(10)y = f.(X) + randn(length(X))/10

We will use three folds assigned randomly:
sets = k_fold_cross_validation_sets(length(X), 3)

Next, we implement our metric. We use the mean squared error:
metric = (f, X, y)->begin
m = length(X)return sum((f(X[i]) - y[i])^2 for i in m)/mend

We now loop through different values of λ and fit different radial basis
functions. We will use the Gaussian radial basis. Cross validation is used toobtain the MSE for each value:
λs = 10 .^ range(-4, stop=2, length=101)
es = []
basis = r->exp(-5r^2)
for λ in λs
fit = (X, y)->regression(X, y, radial_bases(basis, X), λ)push!(es,cross_validation_estimate(X, y, sets, fit, metric)[1])end

The resulting curve has a minimum at λ ≈ 0.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

272 c ha p te r 14 . su rrog a te mod e l s

train(•)

test( fˆ, •)generalization error estimate

bootstrap_sets(m, b) = [TrainTest(rand(1:m, m), 1:m) for i in 1:b]function bootstrap_estimate(X, y, sets, fit, metric)mean(train_and_validate(X, y, tt, fit, metric) for tt in sets)
endfunction leave_one_out_bootstrap_estimate(X, y, sets, fit, metric)m, b = length(X), length(sets)
ε = 0.0
models = [fit(X[tt.train], y[tt.train]) for tt in sets]
for j in 1 : m
c = 0
δ = 0.0
for i in 1 : b
if j ∉ sets[i].train
c += 1δ += metric(models[i], [X[j]], [y[j]])endend
ε += δ/cend
return ε/mend

Figure 14.11. A single bootstrap
sample consists of m indices into
the dataset sampled with replacement. A bootstrap sample is used
to train a model, which is evaluatedon the full dataset to obtain an estimate of the generalization error.

Algorithm 14.11. A method for obtaining b bootstrap samples, eachfor a data set of size m.
Algorithm 14.12. A method forcomputing the bootstrap generalization error estimate by training
and validating on the list of trainvalidate sets contained in sets. The
other variables are the list of design
points X, the corresponding objective function values y, a function
fit that trains a surrogate model,
and a function metric that evaluates a model on a data set.
Algorithm 14.13. A method forcomputing the leave-one-out bootstrap generalization error estimate
using the train-validate sets sets.
The other variables are the list ofdesign points X, the corresponding objective function values y, afunction fit that trains a surrogate
model, and a function metric that
evaluates a model on a data set.

© 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

14 . 5 . mod e l se l e c ti on 273

Unfortunately, the leave-one-out bootstrap estimate introduces a new bias
due to the varying test set sizes. The 0.632 bootstrap estimate15 (algorithm 14.14)
alleviates this bias:(14.36)

ǫ0.632-boot = 0.632ǫleave-one-out-boot + 0.368ǫboot

function bootstrap_632_estimate(X, y, sets, fit, metric)models = [fit(X[tt.train], y[tt.train]) for tt in sets]
ϵ_loob = leave_one_out_bootstrap_estimate(X,y,sets,fit,metric)
ϵ_boot = bootstrap_estimate(X,y,sets,fit,metric)
return 0.632ϵ_loob + 0.368ϵ_bootend

Several generalization estimation methods are compared in example 14.3.
Consider ten evenly spread samples of f ( x ) = x2 + ǫ/2 over x ∈ [−3, 3],
where ǫ is zero-mean, unit-variance Gaussian noise. We would like to test
several different generalization error estimation methods when fitting a linear
model to this data. Our metric is the root mean squared error, which is thesquare root of the mean squared error.
The methods used are the holdout method with eight training samples,five-fold cross validation, and the bootstrap methods each with ten bootstrap
samples. Each method was fitted 100 times and the resulting statistics areshown below.

The 0.632 bootstrap estimate wasintroduced in B. Efron, “Estimating the Error Rate of a Prediction Rule: Improvement on CrossValidation,” Journal of the American Statistical Association, vol. 78,
no. 382, pp. 316–331, 1983. A variant, the 0.632+ bootstrap estimate,
was introduced in B. Efron and
R. Tibshirani, “Improvements on
Cross-Validation: The .632+ Bootstrap Method,” Journal of the American Statistical Association, vol. 92,
no. 438, pp. 548–560, 1997.
15

Algorithm 14.14. A method for
obtaining the 0.632 bootstrap estimate for data points X, objectivefunction values y, number of bootstrap samples b, fitting function
fit, and metric function metric.
Example 14.3. A comparisonof generalization error estimation
methods. The vertical lines in thebox and whisker plots indicate theminimum, maximum, first andthird quartiles, and median of every generalization error estimation
method among 50 trials.

Holdout
5-Fold Cross Validation
Bootstrap

Leave-One-Out Bootstrap
0.632 Bootstrap

0

2

4

6

8

10

12

root mean squared error© 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

274 c ha p te r 14 . su rrog a te mod e l s

14.6 Summary• Surrogate models are function approximations that can be optimized insteadof the true, potentially expensive objective function.
• Many surrogate models can be represented using a linear combination of basis
functions.
• Model selection involves a bias-variance tradeoff between models with low
complexity that cannot capture important trends and models with high complexity that overfit to noise.
• Generalization error can be estimated using techniques such as holdout, k-foldcross validation, and the bootstrap.

14.7 Exercises
Exercise 14.1. Derive an expression satisfied by the optimum of the regression
problem equation (14.8) by setting the gradient to zero. Do not invert any matrices.
The resulting relation is called the normal equation.
Exercise 14.2. When would we use a more descriptive model, for example, withpolynomial features, versus a simpler model like linear regression?
Exercise 14.3. A linear regression problem of the form in equation (14.8) is notalways solved analytically, and optimization techniques are used instead. Why isthis the case?
Exercise 14.4. Suppose we evaluate our objective function at four points: 1, 2,
3, and 4, and we get back 0, 5, 4, and 6. We want to fit a polynomial modelf ( x ) = ∑ik=0 θi xi . Compute the leave-one-out cross validation estimate of themean squared error as k varies between 0 and 4. According to this metric, what isthe best value for k, and what are the best values for the elements of θ?

© 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