13 Sampling Plans
For many optimization problems, function evaluations can be quite expensive.
For example, evaluating a hardware design may require a lengthy fabrication
process, an aircraft design may require a wind tunnel test, and new deep learning
hyperparameters may require a week of GPU training. A common approach for
optimizing in contexts where evaluating design points is expensive is to build asurrogate model, which is a model of the optimization problem that can be efficiently
optimized in lieu of the true objective function. Further evaluations of the true
objective function can be used to improve the model. Fitting such models requires
an initial set of points, ideally points that are space-filling; that is, points that cover
the region as well as possible. This chapter covers different sampling plans forcovering the search space when we have limited resources.1

13.1 Full Factorial
The full factorial sampling plan (algorithm 13.1) places a grid of evenly spaced
points over the search space. This approach is easy to implement, does not rely
on randomness, and covers the space, but it uses a large number of points. A gridof evenly spaced points is spread over the search space as shown in figure 13.1.
Optimization over the points in a full factorial sampling plan is referred to as gridsearch.
The sampling grid is defined by a lower-bound vector a and an upper-boundvector b such that ai ≤ xi ≤ bi for each component i. For a grid with misamples in the ith dimension, the nearest points are separated by a distance( bi − a i ) / ( m i − 1 ) .
The full factorial method requires a sample count exponential in the numberof dimensions.2 For n dimensions with m samples per dimension, we have mn

There are other references that
discuss the topics in this chapter
in greater detail. See, for example:
G. E. P. Box, W. G. Hunter, and J. S.
Hunter, Statistics for Experimenters:
An Introduction to Design, Data Analysis, and Model Building, 2nd ed. Wiley, 2005. A. Dean, D. Voss, and
D. Draguljić, Design and Analysis
of Experiments, 2nd ed. Springer,
2017. D. C. Montgomery, Design
and Analysis of Experiments. Wiley,2017.
1

The full factorial method gets its
name not from a factorial sample
count (it is exponential) but fromdesigning with two or more discrete
factors. Here the factors are the mdiscretized levels associated witheach variable.
2236 c ha p te r 13. sa mp l i n g p l a n s

Figure 13.1. Full factorial searchcovers the search space in a grid ofpoints.

b1 − a1
m1 −1

b2

x2

b2 − a2
m2 −1

a2
a1

x1

b1

total samples. This exponential growth is far too high to be of practical use when
there are more than a few variables. Even when full factorial sampling is able tobe used, the grid points are generally forced to be quite coarse and therefore caneasily miss small, local features of the optimization landscape.
function samples_full_factorial(a, b, m)
ranges = [range(a[i], stop=b[i], length=m[i])for i in 1 : length(a)]collect.(collect(Iterators.product(ranges...)))
end

Algorithm 13.1. A function for obtaining all sample locations for thefull factorial grid. Here, a is a vector of variable lower bounds, b isa vector of variable upper bounds,
and m is a vector of sample counts
for each dimension.

13.2 Random SamplingA straightforward alternative to full factorial sampling is random sampling, whichsimply draws m random samples over the design space using a pseudorandom
number generator. To generate a random sample x, we can sample each variableindependently from a distribution. If we have bounds on the variables, such asai ≤ xi ≤ bi , a common approach is to use a uniform distribution over [ ai , bi ],although other distributions may be used. For some variables, it may make sense
to use a log-uniform distribution.3 The samples of design points are uncorrelatedwith each other. The hope is that the randomness, on its own, will result in an
adequate cover of the design space.

Some parameters, such as thelearning rate for deep neural networks, are best searched in logspace.
3© 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

13. 3. u n i form p roj e c ti on p l a n s 237

13.3

Uniform Projection Plans

Suppose we have a two-dimensional optimization problem discretized into anm × m sampling grid as with the full factorial method, but, instead of taking all
m2 samples, we want to sample only m positions. We could choose the samples
at random, but not all arrangements are equally useful. We want the samples tobe spread across the space, and we want the samples to be spread across eachindividual component.
no variation in
one component

too clustereduniform projection

A uniform projection plan is a sampling plan over a discrete grid where the distribution over each dimension is uniform. For example, in the rightmost sampling
plan in figure 13.2, each row has exactly one entry and each column has exactly
one entry.
A uniform projection plan with m samples on an m × m grid can be constructedusing an m-element permutation as shown in figure 13.3. There are therefore m!
possible uniform projection plans.44= 4 2 1 3 5

For m = 5 this is already 5! = 120
possible plans. For m = 10 there
are 3,628,800 plans.

4

Figure 13.3. Constructing a uniform projection plan using a permutation.

5p

Figure 13.2. Several ways to choosem samples from a two-dimensional
grid. We generally prefer sampling
plans that cover the space and varyacross each component.

3
2
1
1

2

3

4

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

238 c ha p te r 13. sa mp l i n g p l a n s

Sampling with uniform projection plans is sometimes called Latin-hypercube
sampling because of the connection to Latin squares (figure 13.4). A Latin squareis an m × m grid where each row contains each integer 1 through m and eachcolumn contains each integer 1 through m. Latin-hypercubes are a generalization
to any number of dimensions.
Uniform projection plans for n dimensions can be constructed using a permutation for each dimension (algorithm 13.2).
function uniform_projection_plan(m, n)perms = [randperm(m) for i in 1 : n][[perms[i][j] for i in 1 : n] for j in 1 : m]end

13.4 Stratified Sampling

4

1

3

2

1

4

2

3

3

2

1

4

2

3

4

1

Figure 13.4. A 4 × 4 Latin square.
A uniform projection plan can beconstructed by choosing a valuei ∈ {1, 2, 3, 4} and sampling allcells with that value.
Algorithm 13.2. A function for constructing a uniform projection plan
for an n-dimensional hypercube
with m samples per dimension. It
returns a vector of index vectors.

yMany sampling plans, including uniform projection and full factorial plans, arebased on an m × m grid. Such a grid, even if fully sampled, could miss important
information due to systematic regularities as shown in figure 13.5. One method
for providing an opportunity to hit every point is to use stratified sampling.
Stratified sampling modifies any grid-based sampling plan, including full
factorial and uniform projection plans. Cells are sampled at a point chosen uniformly at random from within the cell rather than at the cell’s center as shown infigure 13.6.

xf (x)sampling on gridstratified sampling
Figure 13.5. Using an evenlyspaced grid on a function with systematic regularities can miss important information.

Figure 13.6. Stratified sampling applied to a uniform projection plan.

© 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

13. 5 . sp a c e -fi l l i n g me tri c s 239

13.5 Space-Filling Metrics
A good sampling plan fills the design space since the ability for a surrogate model
to generalize from samples decays with the distance from those samples. Not
all plans, even uniform projection plans, are equally good at covering the search
space. For example, a grid diagonal (figure 13.7) is a uniform projection plan but
only covers a narrow strip. This section discusses different space-filling metrics formeasuring the degree to which a sampling plan X ⊆ X fills the design space.

13.5.1 Discrepancy
The ability of the sampling plan to fill a hyper-rectangular design space can bemeasured by its discrepancy.5 If X has low discrepancy, then a randomly chosen
subset of the design space should contain a fraction of samples proportional to thesubset’s volume.6 The discrepancy associated with X is the maximum difference
between the fraction of samples in a hyper-rectangular subset H and that subset’s
volume:T#( X H)(13.1)− λ(H)d( X ) = supremum#X
H
where #X and #( X H) are the number of points in X and the number of points
in X that lie in H, respectively. The value λ(H) is the n-dimensional volumeof H, the product of the side lengths of H. The term supremum is very similar
to maximization but allows a solution to exist for problems where H merely
approaches a particular rectangular subset, as seen in example 13.1.7
Computing the discrepancy of a sampling plan over the unit hyper-rectangle
is often difficult, and it is not always clear how we can compute the discrepancy
for nonrectangular feasible sets.
T

13.5.2 Pairwise Distances

Figure 13.7. A uniform projection
plan that is not space-filling.
5
L. Kuipers and H. Niederreiter,
Uniform Distribution of Sequences.
Dover, 2012.
6
In arbitrary dimensions, we canuse the Lebesgue measure, which
is a generalization of volume toany subset of n-dimensional Euclidean space. It is length in onedimensional space, area in twodimensional space, and volume inthree-dimensional space.

7
The definition of discrepancy requires hyper-rectangles and typically assumes that X is a finite subset of a unit hypercube. The notionof discrepancy can be extended toallow H to include other sets, suchas convex polytopes.

An alternative method determining which of two m-point sampling plans is more
space-filling is to compare the pairwise distances between all points within eachsampling plan. Sampling plans that are more spread out will tend to have larger
pairwise distances.
The comparison is typically done by sorting each set’s pairwise distances inascending order. The plan with the first pairwise distance that exceeds the otheris considered more space-filling.
© 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

240c ha p te r 13. sa mp l i n g p l a n s

Consider the set:

 
 
 
 
 

1 1
2 1
1 3
9 3
1 1
3 4
X=
, , , ,
, ,
,
,
,
, ,
5 5
5 5
10 5
10 10
50 50
5 5
The discrepancy of X with respect to the unit square is determined by arectangular subset H that either has very small area but contains very many
points or has very large area and contains very few points.

Example 13.1. Computing the discrepancy for a sampling plan overthe unit square. The sizes of therectangles are slightly exaggerated
to clearly show which points they
contain.

x2

1

0

0

x1

1

i
i
h
h
1 2
, 5 , x2 ∈ 51 , 35 , has a volume of 0.12 andThe blue rectangle, x1 ∈ 10
contains 3 points. Its corresponding
discrepancy
measureis thusi 0.38.
i
h
h

9
1
+ ǫ, 10− ǫ , x2 ∈ 15 + ǫ, 54 − ǫ , produces
The purple rectangle, x1 ∈ 10
an even higher discrepancy. As ǫ approaches zero, the volume and the discrepancy approach 0.48 because the rectangle contains no points. Note thatthe limit was required, reflecting the need to use a supremum in the definitionof discrepancy.

© 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

13. 5 . sp a c e -fi l l i n g me tri c s 241

Algorithm 13.3 computes all pairwise distances between points in a sampling
plan. Algorithm 13.4 compares how well two sampling plans fill space using theirrespective pairwise distances.
import LinearAlgebra: norm
function pairwise_distances(X, p=2)m = length(X)[norm(X[i]-X[j], p) for i in 1:(m-1) for j in (i+1):m]
end

Algorithm 13.3. A function for
obtaining the list of pairwise distances between points in sampling
plan X using the L p norm specifiedby p.

function compare_sampling_plans(A, B, p=2)
pA = sort(pairwise_distances(A, p))pB = sort(pairwise_distances(B, p))for (dA, dB) in zip(pA, pB)
if dA < dB
return 1
elseif dA > dB
return -1
end
end
return 0
end

Algorithm 13.4. A function forcomparing the degree to which two
sampling plans A and B are spacefilling using the L p norm specified
by p. The function returns -1 if A is
more space-filling than B. It returns
1 if B is more space-filling than A. It
returns 0 if they are equivalent.

One method for generating a space-filling uniform projection plan is to generate
several candidates at random and then use the one that is most space-filling.
We can search for a space-filling uniform projection plan by repeatedly mutating a uniform projection plan in a way that preserves the uniform projection
property (algorithm 13.5). Simulated annealing, for example, could be used tosearch the space for a sampling plan with good coverage.
function mutate!(X)m, n = length(X), length(X[1])j = rand(1:n)
i = randperm(m)[1:2]X[i[1]][j], X[i[2]][j] = X[i[2]][j], X[i[1]][j]return X
end

Algorithm 13.5. A function for mutating uniform projection plan X,
while maintaining its uniform projection property.

© 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

242

c ha p te r 13. sa mp l i n g p l a n s

13.5.3 Morris-Mitchell Criterion
The comparison scheme in section 13.5.2 typically results in a challenging optimization problem with many local minima. An alternative is to optimize withrespect to the Morris-Mitchell criterion (algorithm 13.6):8
Φq ( X ) =

∑
i

−q
di

!1/q

(13.2)

where di is the ith pairwise distance between points in X and q > 0 is a tunable
parameter.9 Morris and Mitchell recommend optimizing:
minimizeXmaximizeq∈{1,2,5,10,20,50,100}

Φq ( X )(13.3)

function phiq(X, q=1, p=2)
dists = pairwise_distances(X, p)
return sum(dists.^(-q))^(1/q)
end

8
M. D. Morris and T. J. Mitchell,“Exploratory Designs for Computational Experiments,” Journal of
Statistical Planning and Inference,
vol. 43, no. 3, pp. 381–402, 1995.

Larger values of q will give higher
penalties to large distances.

9

Algorithm 13.6. An implementation of the Morris-Mitchell criterion which takes a list of design points X, the criterion parameter q > 0, and a norm parameterp ≥ 1.

Figure 13.8 shows the Morris-Mitchell criterion evaluated for several randomly
generated uniform projection plans.

© 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

13. 5 . sp a c e -fi l l i n g me tri c s 243

Φ1 = 23.7

Φ1 = 24.2

Φ1 = 24.2

Φ1 = 24.3

Φ1 = 24.3

Φ1 = 24.3

Φ1 = 24.3

Φ1 = 24.5

Φ1 = 24.8

Φ1 = 24.9

Φ1 = 25.0

Φ1 = 25.0

Φ1 = 25.3

Φ1 = 25.5

Φ1 = 26.2

Φ1 = 27.2

Φ1 = 33.0

Φ1 = 36.7

Figure 13.8. Uniform projectionplans sorted from best to worst according to Φ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

244

c ha p te r 13. sa mp l i n g p l a n s

13.6 Space-Filling Subsets
In some cases, we have a set of points X and want to find a subset of points
S ⊂ X that still maximally fills X. The need for identifying space-filling subsets ofX arises in the context of multifidelity models.10 For example, suppose we used asampling plan X to identify a variety of aircraft wing designs to evaluate using
computational fluid dynamic models in simulation. We can choose only a subsetof these design points S to build and test in a wind tunnel. We still want S to be
space filling.
The degree to which S fills the design space can be quantified using the maximum distance between a point in X and the closest point in S. This metric
generalizes to any two finite sets A and B (algorithm 13.7). We can use any L p
norm, but we typically use L2 , the Euclidean distance:
dmax ( X, S) = maximize minimizeks − xk p
x∈ X

s∈ S

(13.4)

min_dist(a, B, p) = minimum(norm(a-b, p) for b in B)
d_max(A, B, p=2) = maximum(min_dist(a, B, p) for a in A)

A space-filling sampling plan is one that minimizes this metric.11 Finding aspace-filling sampling plan with m elements is an optimization problem
minimizedmax ( X, S)

subject to

S⊆X

S

10
A. I. J. Forrester, A. Sóbester, andA. J. Keane, “Multi-Fidelity Optimization via Surrogate Modelling,” Proceedings of the Royal Society of London A: Mathematical,
Physical and Engineering Sciences,vol. 463, no. 2088, pp. 3251–3269,
2007.

Algorithm 13.7. The set L p distancemetrics between two discrete sets,
where A and B are lists of design
points and p is the L p norm parameter.
We can also minimize the MorrisMitchell criterion for S.

11(13.5)#S = m
Optimizing equation (13.5) is typically computationally intractable. A brute
force approach would try all d!/m!(d − m)! size-m subsets for a dataset of d design
points. Both greedy local search (algorithm 13.8) and the exchange algorithm (algorithm 13.9) are heuristic strategies for overcoming this difficulty. They typicallyfind acceptable space-filling subsets of X.
Greedy local search starts with a point selected randomly from X and incrementally adds the next best point that minimizes the distance metric. Points areadded until the desired number of points is reached. Because the points are initialized randomly, the best results are obtained by running greedy local searchseveral times and keeping the best sampling plan (algorithm 13.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

13. 7. q u a si -ra n d om se q u e n c e s 245

Algorithm 13.8. Greedy local
search, for finding m-element sampling plans that minimize a distance metric d for discrete set X.

function greedy_local_search(X, m, d=d_max)
S = [sample(X)]for i in 2 : m
j = argmin([x ∈ S ? Inf : d(X, push!(copy(S), x))
for x in X])
push!(S, X[j])endreturn S
end

The exchange algorithm initializes S to a random subset of X and repeatedly
replaces points that are in S with a different point in X that is not already in S
to improve on the distance metric. The exchange algorithm is also typically run
multiple times.
Figure 13.9 compares space-filling subsets obtained using greedy local searchand the exchange algorithm.

13.7

Quasi-Random Sequences

Quasi-random sequences,12 also called low-discrepancy sequences, are often used in
the context of trying to approximate an integral over a multidimensional space:ZX

f (x) dx ≈v mf (x( i ) )
m i∑
=1

12
C. Lemieux, Monte Carlo andQuasi-Monte
Carlo
Sampling.
Springer, 2009.

(13.6)

where each x(i) is sampled uniformly at random over the domain X and v is thevolume of X . This approximation is known as Monte Carlo integration.
Rather than relying on random or pseudorandom numbers to generate integration points, quasi-random sequences are deterministic sequences that fill the
space in a systematic manner so that the integral converges as fast as possible in
the number of points m.13 These quasi-Monte Carlo methods have an error conver√gence of O(1/m) as opposed to O(1/ m) for typical Monte Carlo integration, as
shown in figure 13.10.

13
Pseudorandom number sequences, such as those produced
by a sequence of calls to rand, aredeterministic given a particularseed, but they appear random.
Quasi-random numbers are also
deterministic but do not appear
random.

© 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

246

c ha p te r 13. sa mp l i n g p l a n sfunction exchange_algorithm(X, m, d=d_max)
S = sample(X, m, replace=false)δ, done = d(X, S), false
while !done
best_pair = (0,0)
for i in 1 : m
s = S[i]for (j,x) in enumerate(X)
if !in(x, S)
S[i] = xδ′ = d(X, S)
if δ′ < δ
δ = δ′
best_pair = (i,j)
end
end
end
S[i] = senddone = best_pair == (0,0)
if !donei,j = best_pair
S[i] = X[j]endendreturn S
endfunction multistart_local_search(X, m, alg, k_max, d=d_max)
sets = [alg(X, m, d) for i in 1 : k_max]
return sets[argmin([d(X, S) for S in sets])]
end

Algorithm 13.9. The exchange algorithm for finding m-element sampling plans that minimize a distance metric d for discrete set X.

Algorithm 13.10. Multistart local
search runs a particular searchalgorithm multiple times and returns the best result. Here, X isthe list of points, m is the sizeof the desired sampling plan,alg is either exchange_algorithm
or greedy_local_search, k_max isthe number of iterations to run,and d is the distance metric.

© 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

13. 7. q u a si -ra n d om se q u e n c e s 247

exchange algorithm

x1

x1

Figure 13.9. Space-filling subsets obtained with both greedy local search and the exchange algorithm.

Quasi-random sequences are typically constructed for the unit n-dimensional
hypercube, [0, 1]n . Any multidimensional function with bounds on each variable
can be transformed into such a hypercube. This transformation is implemented
in algorithm 7.9.
Various methods exist for generating quasi-random sequences. Several suchmethods are compared to random sampling in figure 13.12.

100
10−2
10−4
100
102
104
number of samples
rand

13.7.1 Additive Recurrence
Simple recurrence relations of the form:x ( k +1) = x ( k ) + c

relative error

x2

greedy local search(mod 1)(13.7)

produce space-filling sets provided that c is irrational. The value of c leading to
the smallest discrepancy is√
5−1
c = 1− ϕ =≈ 0.618034(13.8)
2
where ϕ is the golden ratio.14
We can construct a space-filling set over n dimensions using an additive recurrence sequence for each coordinate, each with its own value of c. The square rootsof the primes are known to be irrational, and can thus be used to obtain different
sequences for each coordinate:
√
√
√
√
√
c1 = 2, c2 = 3, c3 = 5, c4 = 7, c5 = 11, . . .
(13.9)

Sobol

Figure 13.10. The error from
R1
estimating 0 sin(10x ) dx using
Monte Carlo integration with random numbers from U (0, 1) anda Sobol sequence. The Sobol sequence, covered in section 13.7.3,
converges faster.

C. Schretter, L. Kobbelt, and P.-O.
Dehaye, “Golden Ratio Sequencesfor Low-Discrepancy Sampling,”
Journal of Graphics Tools, vol. 16,
no. 2, pp. 95–104, 2016.

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

248 c ha p te r 13. sa mp l i n g p l a n s

Methods for additive recurrence are implemented in algorithm 13.11.
using Primes
function get_filling_set_additive_recurrence(m; c=φ-1)
X = [rand()]for i in 2 : m
push!(X, mod(X[end] + c, 1))
end
return X
endfunction get_filling_set_additive_recurrence(m, n)
ps = primes(max(ceil(Int, n*(log(n) + log(log(n)))), 13))seqs = [get_filling_set_additive_recurrence(m, c=sqrt(p))
for p in ps[1:n]]
return [collect(x) for x in zip(seqs...)]
end

Algorithm 13.11.
Additiverecurrenceforconstructing
m-element filling sequences overn-dimensional unit hypercubes.
The Primes package is used togenerate the first n prime numbers,
where the kth prime number isbounded byk (log k + log log k)
for k > 6, and primes(a) returns
all primes up to a. Note that 13 is
the sixth prime number.

13.7.2 Halton Sequence
The Halton sequence is a multidimensional quasi-random space-filling set.15 The
single-dimensional version, called van der Corput sequences, generates sequenceswhere the unit interval is divided into powers of base b. For example, b = 2
produces:
1 1 3 1 5 3 7 1
X=, , , , , , ,,...
(13.10)2 4 4 8 8 8 8 16J. H. Halton, “Algorithm 247:
Radical-Inverse Quasi-Random
Point Sequence,” Communicationsof the ACM, vol. 7, no. 12, pp. 701–
702, 1964.

15

whereas b = 5 produces:X=

1 2 3 4 1 6 11
, , , ,,,,...
5 5 5 5 25 25 25

(13.11)

Multi-dimensional space-filling sequences use one van der Corput sequence foreach dimension, each with its own base b. The bases, however, must be coprime16
in order to be uncorrelated. Methods for constructing Halton sequences are
implemented in algorithm 13.12.

Two integers are coprime if the
only positive integer that divides
them both is 1.

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

13. 7. q u a si -ra n d om se q u e n c e s 249

using Primes
function halton(i, b)
result, f = 0.0, 1.0
while i > 0
f = f / b;
result = result + f * mod(i, b)
i = floor(Int, i / b)
end
return result
end
get_filling_set_halton(m; b=2) = [halton(i,b) for i in 1: m]function get_filling_set_halton(m, n)
bs = primes(max(ceil(Int, n*(log(n) + log(log(n)))), 6))
seqs = [get_filling_set_halton(m, b=b) for b in bs[1:n]]
return [collect(x) for x in zip(seqs...)]
end

For large primes, we can get correlation in the first few numbers. Such a correlation is shown in figure 13.11. Correlation can be avoided by the leaped Haltonmethod,17 which takes every pth point, where p is a prime different from allcoordinate bases.

13.7.3 Sobol Sequences
Sobol sequences are quasi-random space-filling sequences for n-dimensional hypercubes.18 They are generated by xor-ing the previous Sobol number with a setof direction numbers:19(k)( i −1)(i )(13.12)⊻ vj
Xj = Xj
where v j is the jth bit of the kth direction number. Tables of good direction
numbers have been provided by various authors.20
A comparison of these and previous approaches is shown in figure 13.12. For
high values several methods exhibit a clear underlying structure.
(k)

Algorithm
13.12.
Haltonquasi-random m-element filling
sequences over n-dimensional unit
hypercubes, where b is the base.
The bases bs must be coprime.

Figure 13.11. The Halton sequencewith b = [19, 23] for which thefirst 18 samples are perfectly linearly correlated.
17
L. Kocis and W. J. Whiten,“Computational
Investigationsof Low-Discrepancy Sequences,”
ACM Transactions on Mathematical
Software, vol. 23, no. 2, pp. 266–294,
1997.

I. M. Sobol, “On the Distributionof Points in a Cube and the Approximate Evaluation of Integrals,”
USSR Computational Mathematics
and Mathematical Physics, vol. 7,
no. 4, pp. 86–112, 1967.
18

The symbol ⊻ denotes the xor operation, which returns true if andonly if both inputs are different.

19

The Sobol.jl package provides
an implementation for up to 1,111
dimensions.
20© 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

c ha p te r 13. sa mp l i n g p l a n s

m = 1000

m = 100

m = 10250

RandomAdditive Recurrence

HaltonUniform Projection

Sobol

Figure 13.12. A comparison of
space-filling sampling plans in two
dimensions. Samples are coloredaccording to the order in which
they are sampled. The uniform
projection plan was generated randomly and is not optimized.

© 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

13. 8 . su mma ry 251

13.8 Summary• Sampling plans are used to cover search spaces with a limited number ofpoints.
• Full factorial sampling, which involves sampling at the vertices of a uniformly
discretized grid, requires a number of points exponential in the number ofdimensions.
• Uniform projection plans, which project uniformly over each dimension, canbe efficiently generated and can be optimized to be space filling.
• Greedy local search and the exchange algorithm can be used to find a subset
of points that maximally fill a space.

13.9 Exercises
Exercise 13.1. Filling a multidimensional space requires exponentially morepoints as the number of dimensions increases. To help build this intuition, determine the side lengths of an n-dimensional hypercube such that it fills half of thevolume of the n-dimensional unit hypercube.
Exercise 13.2. Suppose that you sample randomly inside a unit sphere in ndimensions. Compute the probability that a randomly sampled point is within ǫ
distance from the surface of the sphere as n → ∞. Hint: The volume of a sphere
is C (n)r n , where r is the radius and C (n) is a function of the dimension n only.

hypercube side length ℓ• Quasi-random sequences are deterministic procedures by which space-fillingsampling plans can be generated.

1
0.8
0.6
5
10
number of dimensions, nExercise 13.3. Suppose we have a sampling plan X = {x(1) , . . . , x(10) }, wherex(i) = [cos(2πi/10), sin(2πi/10)](13.13)

Compute the Morris-Mitchell criterion for X using an L2 norm when the parameter q is set to 2. In other words, evaluate Φ2 ( X ). If we add [2, 3] to each x(i) , will
Φ2 ( X ) change? Why or why not?
Exercise 13.4. Additive recurrence requires that the multiplicative factor c inequation (13.7) be irrational. Why can c not be rational?
© 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