Optimal Sliced Latin Hypercube Designs
Shan BA, William R. MYERS, and William A. BRENNEMAN
The Procter and Gamble Company,
Mason, OH 45040(ba.s@pg.com; myers.wr@pg.com; brenneman.wa@pg.com)

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

Sliced Latin hypercube designs (SLHDs) have important applications in designing computer experimentswith continuous and categorical factors. However, a randomly generated SLHD can be poor in terms ofspace-filling, and based on the existing construction method that generates the SLHD column by column
using sliced permutation matrices, it is also difficult to search for the optimal SLHD. In this article, wedevelop a new construction approach that first generates the small Latin hypercube design in each slice andthen arranges them together to form the SLHD. The new approach is intuitive and can be easily adapted
to generate orthogonal SLHDs and orthogonal array-based SLHDs. More importantly, it enables us todevelop general algorithms that can search for the optimal SLHD efficiently.
KEY WORDS:

1.

Computer experiment; Continuous and categorical factors; Maximin distance criterion;
Space-filling design.

INTRODUCTION

Computer experiments, which simulate real-world phenomena using computational methods such as finite element analysis
or computational fluid dynamics, are becoming widely used inalmost every field of scientific research and product development. Different from conducting experiments physically, computer simulations are usually deterministic (not subject to random error) and are often characterized by a large number of
input factors whereas only few of them are important. As a result, replication is unnecessary in computer experiments and agood design should be noncollapsing, in the sense that when thedesign points are projected onto the subspace of these few important factors, replications can be avoided (Santner, Williams,
and Notz 2003; Fang, Li, and Sudjianto 2005). These requirements are perfectly met by the Latin hypercube designs (LHDs;
McKay, Beckman, and Conover 1979), which are most popularly used in designing computer experiments with continuous (quantitative) inputs and possess the desirable property
that when projecting an n-run LHD onto any factor, it always
achieves n different levels for that variable. Suppose we denote the n levels of a factor by 1, . . . , n, then an n-run LHD inp factors can be generated by using a random permutation of{1, . . . , n} for each column of its (n × p) design matrix.
The motivation of this work originates from two computer experiments from The Procter & Gamble Company, whose inputs
not only contain continuous variables but also include one ormore categorical (qualitative) variables. Examples of categorical variables are the type or shape of equipment, the presence orabsence of a part or when a continuous variable is restricted totwo or three discrete levels. The presence of such categorical inputs makes it quite challenging to design computer experiments
because the categorical factor usually cannot accommodate as
many as n different levels and consequently the LHD structure
cannot be directly applied.
For computer experiments with both continuous and categorical inputs, Qian (2012) proposed generating a sliced Latin

hypercube design (SLHD) for the p continuous factors. An nrun SLHD is a special type of LHD that can be partitioned
into t slices (blocks), each of which is also an LHD containing
n/t runs, and then each slice can be used under one of the t different level combinations of categorical factors. This structure
guarantees the maximum one-dimension uniformity of continuous factors in the whole design as well as in each slice. For
example, with two continuous variables (p = 2) and one categorical variable with three levels (t = 3), a 12-run SLHD X forthe two continuous factors is shown (in transpose) in (1), where
the solid lines divide X into three slices X 1 , X 2 , and X 3 . Obviously, X itself is an LHD with 12 different levels and for eachc = 1, . . . , 3, ⌈X c /3⌉ is a smaller LHD with four levels, where⌈ ⌉ denotes the ceiling function such that ⌈a⌉ equals the smallest
integer no less than a. After assigning each slice to one levelof the categorical variable, the full design for both the continuous and categorical inputs can be easily obtained as [X,C].
The benefit of using an SLHD is manifest. For instance, even
if the categorical factor turns out to be insignificant, the design
projected onto the remaining continuous factors still forms a
12-run LHD with no projection redundancy.
!
"
7 12 1 6 9 2 10 5 3 4 11 8T,X =4 9 3 11 1 6 12 7 10 2 5 8$#(1)CT = 1 1 1 1 2 2 2 2 3 3 3 3 .
A construction method was provided by Qian (2012) to randomly generate an n-run SLHD in t slices with each slice con-479

© 2015 American Statistical Association and
the American Society for Quality
TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4
DOI: 10.1080/00401706.2014.957867
Color versions of one or more of the figures in the article can be foundonline at www.tandfonline.com/r/tech.

480

SHAN BA, WILLIAM R. MYERS, AND WILLIAM A. BRENNEMAN

10
8
6
4
2
0

0

2

4

6

8

10

12

(b)

12

(a)

0

2

4

6

8

10

12

0

2

4

10

12

8

10

12

12
10
0

2

4

6

8

10
8
6
4
2
0

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

8

(d)

12

(c)

6

0

2

4

6

8

10

12

0

2

4

6

Figure 1. Several random SLHDs (n = 12, p = 2, m = 4, t = 3), where symbols (◦, !, △) represent points belonging to the three different
slices.

sisting of m points (n = mt). To impose the sliced structure,
Qian (2012) proposed determining the p random permutationsof {1, . . . , n} for each column of the SLHD by p-independent
sliced permutation matrices (SPMs) H 1 , . . . , H p . In this way,
the construction effort of SLHD goes into the generation of psuch SPMs, where each H j is defined to be an (m × t) matrix containing {1, . . . , n} as its elements and each column of⌈H j /t⌉ forms a permutation of {1, . . . , m}. Qian (2012) proposed randomly generating each SPM in three steps: (i) fill in theith (i = 1, . . . , m) row of SPM with {a ∈ {1, . . . , n}|⌈a/t⌉ =i}; (ii) randomly permute the entries in each row of SPM; (iii)
randomly permute the entries in each column of SPM. After
stacking the columns of each H j into a single vector, Qian
(2012) showed that it can be used as the jth column in the
SLHD design matrix.
Figure 1(a) depicts the 12-run two-dimensional SLHD matrix
X in (1), whose two columns correspond to the following two
SPMs:
⎡
⎡
⎤
⎤
7
9
3
4
1 10
⎢ 12 2
⎢9
4⎥
6
2 ⎥
⎥,
⎥ . (2)H1 = ⎢
H2 = ⎢⎣1⎣3
10 11 ⎦
12 5 ⎦
6
5
8
11 7
8

In Figure 1(b)–1(d), we also provide three other random SLHDs of the same size that are produced by different random
SPMs. Although they also achieve maximum one-dimensional
uniformity, all of them are inferior to the design in Figure 1(a)for various reasons. In Figure 1(b), the two factors are perfectly
correlated in the whole design as well as in each slice. If thisTECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

design was used, the effects of the two factors would become
indistinguishable. In Figure 1(c), although overall the points
have been spread out in the design region, in one slice (denotedby symbol !) the two factors are still highly correlated. This
design is undesirable because if the responses across different
categories (slices) have very weak correlation, then the modelcannot borrow too much information among the data from different categories resulting in poor predictions based mainly on
this slice. When it comes to Figure 1(d), the small design in eachslice is desirable, but the design points for the whole SLHD get
clustered together, which leaves a large area of the design region
unexplored. This design would be very inefficient if the categorical variable turns out to be insignificant and the design points
need to be projected onto the subspace of continuous variables.
In general, a desirable SLHD for computer experiments shouldhave its design points well spread out for the whole design aswell as for each slice. However, for an n-run design consisting of t slices with each slice containing m points in p factors,
there can be as many as ((m!)t (t!)m )p different SLHDs. How to
avoid the undesirable ones such as in Figure 1(b)–1(d) and how
to efficiently find the optimal SLHD is the main focus of thisarticle.
Extensive work has been done in the literature on finding“good” LHDs. Due to the huge combinatorial nature of sucha problem, a popular strategy is to use an exchange algorithm
to iteratively search in the design space and optimize the LHD
based on some criteria. For example, in a simulated annealing
algorithm proposed by Morris and Mitchell (1995), two randomly chosen elements within a randomly selected column areDownloaded by [Ohio State University Libraries] at 13:49 19 November 2015

OPTIMAL SLHDS

exchanged in each step to find a new design. Other exchange
algorithms have been discussed by Jin, Chen, and Sudjianto(2005) and Joseph and Hung (2008). Extending these algorithms
to find optimal SLHDs, however, is not straightforward. This isbecause in the construction method proposed in Qian (2012),
permuting elements in any row of an SPM must precede therandom shuffling of its column elements. After elements in any
column of an SPM are permuted, the elements in its rows can nolonger be freely exchanged. This constraint on the order of permutations creates obstacles for freely generating new designs
in the neighborhood of an existing SLHD, which prohibits theimplementation of exchange algorithms for combinatorial optimization. For example, suppose we want to further improve thecurrent SLHD shown in Figure 1(a). If we choose to search in
the neighboring design region by exchanging elements “7” and“3” in the first row of H 1 in (2), it would immediately violate
the definition of SPM and thus cannot generate a valid candidateof the SLHD.

In this work, we develop a new construction method for the
SLHD that is not based on the SPMs. Unlike in Qian (2012)that generates the whole design column by column, we propose
to first construct the small LHD for each slice and then arrange
their levels to form an overall SLHD. This new construction
method is intuitive and it can be used to generate SLHDs withmany desirable properties. The rest of this article is organized
as follows. In Section 2, we present the new construction approach for the SLHD and also point out two of its interesting
variants that can generate SLHDs with orthogonal columns orwith higher-dimensional uniformity. In Section 3, a general algorithm for searching for the optimal SLHD is proposed, and
we also develop strategies to substantially improve its efficiency.
Two examples for designing computer experiments at The Procter & Gamble Company are provided in Section 4, and finally
we give some concluding remarks in Section 5.

2.

A NEW CONSTRUCTION METHOD

In this section, we present a simple and intuitive method
to construct an n-run SLHD with t slices, where each slice
contains m points in p factors (p, n, t, m are positive integers
and n = mt). This construction method consists of only two
steps:

481

Step 1. Independently generate t small LHDs X 1 , . . . , X t for thet slices, where each LHD contains m points in p factors.
Denote their factor levels by 1,+. . . , m and stack them toform an (n × p) matrix X = ti=1 X i .
Step 2. Independently in each column of the matrix X, replace
the t entries of level l (l = 1, . . . , m) with a random
permutation !l of elements {(l − 1)t + 1, . . . , lt}.
It is quite obvious that the design matrix X constructed above
is an n-run SLHD with t slices. Unlike the procedure introducedby Qian (2012) that generates the whole design column bycolumn, we first construct the small design in each slice andthen arrange them together. As a result, we will refer to thisnew method as slice-wise construction method. An example isprovided below to illustrate the proposed procedure, where we
generate three small LHDs X 1 , X 2 , X 3 and arrange their levels
to form a 12-run SLHD:

Comparing this design with Figure 1(a), it can be seen that this
SLHD is the same as the one generated by the SPMs in (2)using the procedure proposed by Qian (2012). In general, thestatistical equivalency of the two methods in generating random
SLHDs can be established by recognizing the following facts.
For any SLHD generated by the proposed slice-wise method,
we can construct p SPMs by assigning the jth column of theconstructed design matrix in slice X i to the ith column of the
SPM H j . Then permuting elements in the jth column of the ith
small LHD in our Step 1 is equivalent to permuting entries in
the ith column of the jth SPM H j in Qian (2012). In addition,
permuting the t elements in each !i in our Step 2 is equivalent
to permuting elements in each row of the (initial) SPMs in Qian
(2012). Due to this equivalency, all the sampling propertiesof the SLHD discussed in Qian (2012) are also true for the
SLHDs generated with the new method. Unlike the procedure
introduced by Qian (2012), the proposed slice-wise construction
has no constraint on the order of its two steps. After completing
Step 2, we can still go back to Step 1 to modify the design in any
slice and then use the same permutations !l (l = 1, . . . , m) in
Step 2 to generate new SLHDs in the neighborhood of theprevious SLHD. This property is very important in applying theexchange algorithms for finding optimal SLHDs, which will be
discussed in detail in Section 3.2. In addition, by envisioning theslice-wise construction method as a two-stage sampling process,we can further develop strategies to substantially improve the
TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

482

SHAN BA, WILLIAM R. MYERS, AND WILLIAM A. BRENNEMAN

efficiency of generating space-filling SLHDs, and this topic will
be exploited in Section 3.3.
Before concluding this section, we briefly point out two interesting variants of the proposed slice-wise construction method,which can produce some random SLHDs with improved properties. (1) As an improvement on random LHDs that only stratify univariate margins, Owen (1992) and Tang (1993) independently proposed to use orthogonal arrays (OAs; Hedayat,
Sloane, and Stufken 1999; Wu and Hamada 2009) to generate
OA-based LHDs, which can achieve uniformity in more than
one dimensions. Because concatenating several OAs row by row
together leads to a larger OA with similar strength, by using t
OA-based LHDs (constructed as in Tang 1993) for the t slices
in Step 1, the proposed slice-wise construction method can produce an OA-based SLHD in which the whole design as well as
each of its slices can all achieve two- or higher-dimensional uniformities. (2) Another type of improvement over a random LHD
is to construct the orthogonal LHD (Ye 1998) in which every
pair of its columns is required to have zero correlation and/or
the second-order orthogonal LHD in which the element-wiseproduct of every two columns also need to be orthogonal to all
the columns in the design (Sun, Liu, and Lin 2009). By choosing some second-order orthogonal LHDs for the t slices in Step
1 and transforming their levels wisely in Step 2, the proposedslice-wise construction method can also be adapted to produce
the second-order orthogonal SLHDs in which the small LHD
in each slice as well as the whole LHD are all second-order orthogonal. Designs generated by this approach can achieve more
flexible run sizes and accommodate larger number of columns
than the existing ones in the literature (Yang et al. 2013). Since
the orthogonal SLHD may not necessarily be space-filling, weonly give a brief description of their construction in the appendix. We next study how to improve the space-filling propertyof an SLHD, which is more important in designing computer
experiments.
3.

OPTIMAL SLHD

For any given m, t, p values, the slice-wise construction
method proposed in Section 2 can generate ((m!)t (t!)m )p possible SLHDs in total. Among them, we define the optimal SLHD
as the one that optimizes a particular design optimality criterion.
In this section, we first discuss how to develop performancemeasures for evaluating the goodness of an SLHD, and then in
Section 3.2, we propose an algorithm that can generally be applied for finding the optimal SLHD with respect to any criteria.
In Section 3.3, we further develop a very efficient algorithm forfinding the optimal space-filling SLHDs.
3.1

Optimality Criteria for SLHDFor designing computer experiments, one of the most popularly used criteria is the maximin distance criteria (Johnson,
Moore, and Ylvisaker 1990). By definition, this space-filling
criterion tries to spread out the points in the design region insuch a way that the minimum distance among the design points
is maximized. Suppose we denote the distance,p between any two
design points x i and x j as d(x i , x j ) = { k=1 |xik − xj k |q }1/q ,in which q = 1 and q = 2 correspond to the rectangular and EuTECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

clidean distances, respectively. Because many different designs
may have the same minimum interpoint distance, an extensionof this criterion is to minimize the average reciprocal interpoint
distance of the design X = {x 1 , . . . , x n } (Morris and Mitchell
1995; Jin, Chen, and Sudjianto 2005), given by
⎞ 1r
⎛
/
1
2
⎠ .
φr (X) = ⎝n(n − 1) 1≤i<j ≤n d(x i , x j )r

It can be seen that in a special case when r → ∞, minimizing
the above φr is equivalent to maximizing the minimum distanceamong the design points.
When extending the maximin distance criterion for evaluating the goodness of an SLHD X = (X T1 , . . . , X Tt )T , we need toconsider both the space-filling properties of the whole design
as well as that of the small design in each slice. As a result, themaximin distance SLHD not only should minimize the φr (X) for
all the design points, but it should also minimize the φr (X i ) for
each slice (i = 1, . . . , t). To solve this multi-objective optimization problem, we propose a single objective function that takes
a weighted average of all the above individual measures asfollows:
3
2t
1
1/
φr (X i ) .
(3)
φMm (X) =φr (X) +2t i=1Here, more weights are given to φr (X) because the whole design
matrix X contains t times as many points as each X i does
and thus its space-filling property is more important. Based on
this combined measure, we can define the maximin distance
SLHD as the one that minimizes the above φMm (X). Similar
structure of the objective function has also appeared in Jonesand Nachtsheim (2011).
Other popular criteria for designing computer experiments
can be extended for selecting the SLHD in a similar way. For example, we can define a uniform
SLHD as the one that minimizes,
φunif (X) = (φCD (X) + ti=1 φCD (X i )/t)/2, where φCD represents the centered L2 -discrepancy measure proposed by Hickernell (1998). In the next section, we develop a general algorithm
to construct the optimal SLHD with respect to any criterion
defined above.
3.2 A General Algorithm
In the literature, a version of the simulated annealing algorithm proposed by Morris and Mitchell (1995) has been widelyused for constructing optimal LHDs. It uses an exchange algorithm to iteratively search in the design space in such a waythat in each step two randomly chosen elements within a randomly selected column in the design matrix are interchanged.
When it comes to the SLHDs, the construction process becomes more complex. In the slice-wise construction method
proposed in the beginning of Section 2, we can see that the
randomness in generating an SLHD comes from two sources:(i) the random permutation of elements in each column of theX i (i = 1, . . . , t) in Step 1; (ii) the random permutation of elements in each !l (l = 1, . . . , m) in Step 2. In this section, we
present a general algorithm to search for the optimal SLHDs,
which uses an exchange algorithm to randomly perturb both of

OPTIMAL SLHDS

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

these two types of permutations. Suppose we want to minimizea prescribed performance measure φ(X) for the SLHD and let
p0 = 0.5. Our algorithm starts with a randomly chosen SLHD
and then examines a sequence of new designs, each generated
as a perturbation of the preceding one:
Step (i). Denote the current SLHD as X. Draw a random variatez ∼ Uniform(0, 1).
Step (ii). If z ≤ p0 , randomly select a slice in X, and interchange two randomly chosen elements within a randomly selected column in this slice. Denote the new
SLHD as X try and go to Step (iv).
Step (iii). If z > p0 , randomly choose a column in X, and interchange two randomly chosen elements within arandomly selected !l (l ∈ {1, . . . , m}) for this column. Denote the resulting SLHD as X try and go to
Step (iv).
Step (iv). If φ(X try ) < φ(X), replace the current X with X try ;otherwise, replace X with X try with probability π =exp{−[φ(X try ) − φ(X)]/T }, where T is a preset parameter known as “temperature.”
Step (v). Repeat Steps (i) to (iv) until some convergence requirements are met.
All the parameters in the above algorithm are set to be the sameas those in a standard simulated annealing algorithm for which
the convergence is already established (Lundy and Mees 1986).
As a result, this proposed algorithm is guaranteed to ultimately
converge to the global optimum. In the end, the SLHD with thesmallest φ(X) value found by the algorithm is reported as our
optimal design.
We note that since each iteration of the above algorithm onlyinterchanges two elements within a column, its objective function φ(X) can usually be efficiently updated without recalculating all of its components at each step of the algorithm. Jin,
Chen, and Sudjianto (2005) discussed the computational efficiency of evaluating several design optimality criteria. As an
illustration, now we derive the updating formula of the mostcommonly used maximin distance measure φMm (X) for our algorithm. As defined in (3), the value of φMm (X) depends on
the n(n − 1)/2 interpoint distances d(x i , x j ) (1 ≤ i < j ≤ n)among the n design points {x 1 , . . . , x n }. Whenever two elements xwh and xvh are interchanged within the hth column ofthe design matrix, only the interpoint distances that are relatedwith the points x w and x v need to be updated:d ∗ (x w , x k ) = ([d(x w , x k )]q + s(w, v, h, k))1/q ,for 1 ≤ k ≤ n, k ̸= w, v∗
d (x v , x k ) = ([d(x v , x k )]q − s(w, v, h, k))1/q ,for 1 ≤ k ≤ n, k ̸= w, v,

(4)
(5)

where d ∗ (·) represents the new distance and d(·) represents the
distance from the previous step, s(w, v, h, k) = |xvh − xkh |q −|xwh − xkh |q and all the other interpoint distances remain unchanged. Based on this observation, it can be shown that the∗(X) can be efficiently calculated using the previous
new φMmφr (X) and φr (X i ) values as follows:
3
2
t
1
1/ ∗
∗
∗
φMm (X) =φ (X i ) ,(6)φr (X) +2
t i=1 rin whichφr∗ (X) =483!
n(n − 1)2
[φr (X)]rn(n − 1)2/
# ∗$d (x w , x k )−r − d(x w , x k )−r+41≤k≤n,k̸=w,v+/

1≤k≤n,k̸=w,v

⎤⎞1/r
# ∗
$
d (x v , x k )−r − d(x v , x k )−r ⎦⎠ , (7)and (i) when z ≤ p1 , suppose the two points x w and x v belong
to the cth slice (c ∈ {1, . . . , t}). Then we have: w, v ∈ Ic ={(c − 1)m + 1, (c − 1)m + 2, . . . , cm}, andφr∗ (X i )⎧φr (X i ),if i ̸= c,⎪
9⎪# m(m−1)⎪⎪2r⎨ m(m−1)[φ (X i )]
2 # r
$
(8)
= +,
∗
d
(x w , x k )−r − d(x w , x k )−r⎪k∈Ic ,k̸=w,v
⎪
:;
⎪
# ∗
$ 1/r
⎪
⎩ +,
−r
− d(x v , x k )−r, if i = c.
k∈Ic ,k̸=w,v d (x v , x k )(ii) When z > p2 , suppose x w belongs to slice cand x v belongs to slice c′ (c, c′ ∈ {1, . . . , t}). Then
we have w ∈ Ic = {(c − 1)m + 1, (c − 1)m + 2, . . . , cm}, v ∈Ic′ = {(c′ − 1)m + 1, (c′ − 1)m + 2, . . . , c′ m}, andφr∗ (X i )⎧φr (X i ),if i ̸= c, c′⎪
9⎪# m(m−1)⎪2r⎪⎪[φr (X i )]⎪ m(m−1)
2
⎪
⎪
⎨ ,
# ∗
$:;1/r−r−r, if i = c,
= 9+ k∈Ic ,k̸=w d (x w , x k ) − d(x w , x k )⎪# m(m−1)⎪2r⎪[φr (X i )]⎪⎪m(m−1)
2⎪⎪#$:;1/r⎪,⎩+d ∗ (x , x )−r − d(x , x )−r, if i = c′ .
k∈Ic′ ,k̸=v

v

k

v

k

(9)

By using the above updating formulas, considerable computation can be saved in evaluating the objective functions in each
iteration of our proposed algorithm.
3.3

Efficient Two-Stage Algorithm for Generating
Space-Filling SLHDs

In the algorithm presented in Section 3.2, the performancemeasure φ(X) actually is very generally defined and may notnecessarily be restricted to the space-filling measures. When n or
p is large, however, the convergence may be slow due to the enormous number of candidate designs. In this section, we show thatwhen our interest focuses on finding the optimal space-filling
SLHDs, we can modify the algorithm to significantly reduce therequired computation and improve its efficiency. The relevant
design criteria include, for example, the maximin distance criterion φMm (X) or the uniformity criterion φunif (X), which aremost commonly used in designing computer experiments.
The key to improving the efficiency of our algorithm is to perceive the proposed construction method for SLHD as a two-stage
sampling process. Particularly, the two steps in our slice-wise
construction method (described in the beginning of Section 2)
can be envisioned as follows: Step 1 partitions the p-dimensional
input region into mp cells using an (<m × ·=>· · × m?) coarser grid. To
p

TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

484

SHAN BA, WILLIAM R. MYERS, AND WILLIAM A. BRENNEMAN

generate an m-run LHD, a set of m cells are sampled from the
mp population of cells in such a way that the projections ofthese m cells onto any axis fall into m equally spaced intervals.
Such m cells are sampled repeatedly (with replacement) fromthe coarser grid for t times to form t independent m-run LHDs.
In total, a sample of n = mt cells are chosen in this step, andsome of them may be the same. In Step 2, the input region ispartitioned using a much finer <(n × ·=>· · × n)? grid, which further di-pcan be sampled more than once for different slices. In otherwords, an efficient Stage-I algorithm should quickly eliminate
all the undesirable designs that+contain replicated rows in the
combined design matrix X = ti=1 X i . In light of this, the improved two-stage algorithm can be stated as follows:Stage-I algorithm. The algorithm starts with an (n × p) design matrix, which is comprised of t randomly chosen (m ×p) LHDs. In each+ step, suppose the current combined design
matrix is X = ti=1 X i , and the interpoint distances among its
n design points are d(x i , x j ) (1 ≤ i ̸= j ≤ n). Denote S(X) =,1≤i<j ≤n 1{d(x i , x j ) = 0}, where 1{A} = 1 if A is true and1{A} = 0 otherwise. A positive S(X) indicates that some rows
in the current design matrix X are identical, and we refer to themas the duplicated rows. This algorithm proceeds as follows:(i) If S(X) = 0, directly go to Step (iv); otherwise, go to Step(ii) to reduce the number of duplicated rows in the combined design matrix X.
(ii) Randomly pick a duplicated row from the current X, andthen randomly select another row within the same sliceof this row. Interchange the two elements that belong toa randomly chosen column of these two rows. Denote the
new combined design matrix as X try , and update the corresponding interpoint distances using formulas (4) and (5).
Stage II

0

2

1

4

6

2

8

3

10

4

12

Stage I

0

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

pvides each large cell in Step 1 into t p small subcells. One small
subcell is sampled within each of the n sampled large cells from
Step 1 in such a way that the projections of these n selected small
subcells onto any axis fall into n equally spaced small intervals.
Overall, these n selected subcells eliminate the one-dimensional
projection redundancy and constitute a random SLHD. Figure 2
illustrates this two-stage sampling process using a coarser gridand a finer grid, respectively. It can be seen that since the sampling of subcells in Step 2 only takes place locally within thepreviously selected large cells, it has a much smaller impact on
the overall space-filling property of the SLHD. In addition, toobtain n space-filling subcells in Step 2, its n mother cells from
Step 1 must be well spread out. When n < mp , this means thatthese n sampled large cells in Step 1 must be space-filling anddifferent slices cannot sample the same cell from the coarser
grid.
The above observations motivate the development of an efficient two-stage algorithm that optimizes the two steps ofthe slice-wise construction method sequentially (instead ofsimultaneously) for finding the space-filling SLHDs. In theStage-I algorithm, we focus on how to construct the t small
LHDs X 1 , . . . , X t so that their combined design matrix X =+t
i=1 X i is optimal with respect to a chosen space-filling measure φ(X). After fixing these optimal small LHDs, we canrun a Stage-II algorithm that only optimizes the permutations!l (l = 1, . . . , m) for each of the p columns in X. Through thisway, the candidate SLHDs that we need to search have been
greatly reduced from ((m!)t (t!)m )p into a much smaller subset of only (m!)tp + (t!)mp designs. This can lead to substantial
savings in the computation.
A simple way to implement the above two-stage algorithm
is to run the general algorithm presented in Section 3.2 twice:
we first start with a randomly generated SLHD and run thealgorithm with parameter p0 fixed as 1; after obtaining the “optimized” SLHD, we use it as the starting design and run thealgorithm again in which we fix the parameter p0 as 0. It caneasily be seen that the first stage of this procedure optimizes thesmall LHDs in the t slices and then in the second stage it only
improves the !l permutations.
When mp > n, the above Stage-I algorithm can be further
sped up by cutting out a significant amount of undesirable designs from its candidate set. Note that the assumption mp > n isquite general since it can easily be met when the input dimension
p is moderate or large. Because computer experiments are usually characterized by a large number of inputs, it is quite a standard for us to have mp ≫ n in practice. As we have discussed
earlier, since the sample of n cells chosen by the Stage-I algorithm needs to be space-filling, no cell from the (<m × ·=>· · × m?) grid

0

1

2

3

4

0

2

4

6

8

10

12

Figure 2. Illustration of the slice-wise construction method as a two-stage sampling process.
TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

OPTIMAL SLHDS

It can be seen that after quickly eliminating the initial row
duplications in Steps (i) to (iii), the above algorithm never considers any candidate design that contains duplicated rows. Compared to running the general algorithm in Section 3.2 with p0 set
as 1, the above Stage-I algorithm cuts out all the candidate designs with positive S(X) from the search space, and thus it canbe expected to have a much faster convergence rate.
Stage-II algorithm. After obtaining the optimal combined
design matrix X from the Stage I, we can use it (with somerandomly chosen !l ’s) as the starting design and run the generalalgorithm in Section 3.2 with p0 set as 0. This Stage-II algorithm
optimizes the permutations !l (l = 1, . . . , m) for each columnof the matrix X, which arrange the t small space-filling LHDs
found by the Stage-I algorithm into an optimal space-filling
SLHD.
As we have discussed in the beginning of this section, since
the Stage-I algorithm plays a much more important role in determining the space-filling property of the SLHD, more resources
should be allocated for running the Stage-I algorithm if our computational budgets are tight. When p is very large and mp ≫ n,
we may even completely focus on the Stage-I algorithm and skip
the optimization in Stage II by directly using some randomly
generated !l ’s. In such cases, because the n sampled cells fromthe Stage-I algorithm are very sparse (far away from each other)× · · · × m) grid, the Stage-II algorithm has a very minor
in the (m<=>?
peffect in determining the space-filling property of the SLHD.
4.

10
8
6
4

OptimalSLHD

02

OptimalSLHD

Density (1000 random SLHDs)

10

15

(b)

5

Density (1000 random SLHDs)

EXAMPLES

Example 1. A new bottle design was needed for a beauty product of The Procter & Gamble Company that is shipped underchanging air pressure due to altitude changes and the solution ismixed by the consumer through shaking prior to use. The newbottle design was optimized through a computer simulation ofboth the air pressure changes and the solution mixing processwith the goal to minimize bottle deformation under pressure
while at the same time maximize mixing properties. There were
five continuous variables and three two-level categorical variables in the computer experiment. The computer experimentalong with the subsequent analysis allowed the project team todefine a new bottle design that was strong enough to handlechanges in air pressure seen through shipping, had very good
mixing properties, and was aesthetically pleasing.
Based on the available resources a 256-run computer experiment was proposed. Since the three categorical factors haveeight different level combinations, we constructed an SLHD
containing eight slices, each of which is a 32-run LHD for thefive continuous factors (n = 256, m = 32, t = 8, and p = 5).
To demonstrate the advantage of using an optimal SLHD over arandom SLHD, in this example we randomly generate SLHDs
for 1000 times and draw the density plot of their minimum interpoint distances in Figure 3(a). Then a maximin distance SLHDof the same size is constructed using the algorithm in Section 3
and its corresponding minimum interpoint distance is depictedby the red arrow in Figure 3(a). It can be seen that even the best
random SLHD out of the previous 1000 samples is substantiallyworse than the optimal SLHD. This clearly shows the necessity and benefits for optimizing the SLHD using the proposed
algorithm. Since it is well known that the space-filling property of a design has a huge impact on the predictive accuracyof the surrogate model (Johnson, Moore, and Ylvisaker 1990),
we can expect to obtain a much more accurate surrogate model
when using the proposed optimal design instead of a random
design.

(a)0

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015(iii) If S(X try ) < S(X), replace X with X try and go to Step (i);otherwise, discard X try and go back to Step (ii).
(iv) Randomly select a slice from the current X, and interchange two randomly chosen elements within a randomly
selected column in this slice. Denote the new combined
design matrix as X try , and update the corresponding interpoint distances using formulas (4) and (5).
(v) If S(X try ) > 0, discard X try and go back to Step (iv); otherwise, go to Step (vi).
(vi) If φ(X try ) < φ(X), replace the current X with X try ;otherwise, replace X with X try with probability π =exp{−[φ(X try ) − φ(X)]/T }, where T is a preset parameter known as “temperature.”(vii) Repeat Steps (iv) to (vi) until some convergence requirements are met.

485

0.05

0.15

0.25

0.35

Overall minimum interpoint distances

0.4

0.5

0.6

0.7

0.8

0.9

1.0

Average of three minimum distances

Figure 3. Density plots for 1000 random SLHDs: (a) Example 1: overall minimum interpoint distances. (b) Example 2: average of theminimum interpoint distances of the designs in three slices.
TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015

486

SHAN BA, WILLIAM R. MYERS, AND WILLIAM A. BRENNEMAN

Example 2. Engineers at the Procter & Gamble employedcomputer experiments to model the packing line process. The
following experiment involves modeling a particular part (ortransformation) of an oral care packing line. The computer
simulation consisted of nine continuous variables and one threelevel categorical variable representing three types of a particularequipment part. A 132-run SLHD containing three slices, eachof which is a 44-run LHD for the continuous factors was used(n = 132, m = 44, t = 3, and p = 9). Interestingly, after thedata were collected, it was found that the correlations of responses across different categories were quite weak. Since themodel cannot borrow too much information among the datafrom different categories, any new prediction needs to be mademainly based on the information from a single slice. This clearly
demonstrates the importance of also having a good space-fillingdesign for each slice. In Figure 3(b), we compare the average ofthe minimum interpoint distances of the designs in three slices
for an optimal SLHD with that for 1000 random SLHDs of thesame size. We can see that the optimal SLHD is also muchsuperior in this aspect.
5.

CONCLUSIONSIn this article, we developed a new approach for constructing
the SLHDs in which we first generate the small LHD in eachslice and then arrange them together. Unlike the method introduced by Qian (2012), the proposed slice-wise construction
method can be easily adapted to produce SLHDs with many
desirable properties. Not only can it construct the OA-based
SLHDs, but it can also be extended to generate the secondorder orthogonal SLHDs with more flexible run sizes and larger
number of factors than those existing designs in the literature.
More importantly, a general algorithm has been developed basedon the proposed slice-wise construction method for finding the
optimal SLHDs with respect to any criterion. When our interests focus on generating the space-filling SLHDs, a two-stage
strategy has also been presented to substantially improve theefficiency of our algorithm. An R package “SLHD” for implementing the proposed algorithm can be downloaded fromhttp://cran.r-project.org/
APPENDIX: CONSTRUCTION PROCEDURE FOR THE
SECOND-ORDER ORTHOGONAL SLHDFor simplicity, we center the n equally spaced levels in LHD tohave mean equal zero, so that the orthogonality of any two columns
only requires their inner product to be zero. The proposed procedure
is based on Sun, Liu, and Lin (2009), which can generate any secondorder orthogonal LHD in 2c+1 or 2c+1 + 1 runs (c ≥ 1) with 2c columns
as follows:
For c = 1, let
4
@
4
@
1 1
1 2
, T1 =
.
S1 =
1 −1
2 −1
For c ≥ 2, iteratively define Sc and T c as@
4
Sc−1 −S∗c−1,
Sc =∗Sc−1 Sc−1@
4−(T ∗c−1 + 2c−1 S∗c−1 )
T c−1
,
Tc =
∗
c−1
T c−1
T c−1 + 2 Sc−1
TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4

where the * operator works on any matrix with an even number ofrows by multiplying the entries in the top half of the matrix by −1 andleaving those in the bottom half unchanged. Sun, Liu, and Lin (2009)
showed that L c = (T Tc , 02c , −T Tc )T is a second-order orthogonal LHD
in 2c+1 + 1 runs with 2c columns, where 02c denotes a (2c × 1) column vector with all elements equal zero, and L c = (H Tc , −H Tc )T is
a second-order orthogonal LHD in 2c+1 runs with 2c columns, whereH c = T c − Sc /2.
Based on Sun, Liu, and Lin’s (2009) method, the slice-wise construction procedure can be modified as follows to generate the n-run secondorder orthogonal SLHD consisting of t slices, in which each slice isa small second-order orthogonal LHD in m = 2c+1 or 2c+1 + 1 runs(c ≥ 1 and n = mt):(i) For m = 2c+1 , any t ≥ 1 and p ≤ 2c (∀c ≥ 1):
Step 1. Use Sun, Liu, and Lin’s (2009) method to independently
generate t small (m × p) orthogonal LHDs X 1 , . . . , X t as
the design matrices for the t slices. Code their factor levels
as {−2c + 0.5, . . . , −1.5, −0.5, 0.5, 1.5, . . . , 2c − 0.5}, andcombinethem row by row to form an (n × p) matrix X =+t
i=1 X i .
Step 2. Within each X i (for i = 1, . . . , t), transform every positive element l to (l − 0.5)t + (i − 0.5) and transform every negative
element l to (l + 0.5)t − (i − 0.5).
′′′(ii) For m = 2c+1 + 1, t = 2c +1 or 2c +1 + 1 and p ≤ 2min(c,c ) (∀c, c′ ≥1):
Step 1. Use Sun, Liu, and Lin’s (2009) method to independently generate t small (m × p) orthogonal LHDs X 1 , . . . , X t for the tslices. Code their factor levels as {−2c , . . . , −1, 0, 1, . . . , 2c },and combine
them row by row to form an (n × p) matrix+
X = ti=1 X i .
Step 2. Within each X i (for i = 1, . . . , t), transform every positive element l to (l − 0.5)t + (i − 0.5) and transform every negative
element l to (l + 0.5)t − (i − 0.5).
Step 3. Use Sun, Liu, and Lin’s (2009) method to construct a (t × p) orthogonal LHD whose levels are taken as {1 − (1 + t)/2, 2 −(1 + t)/2, . . . , t − (1 + t)/2}. For j = 1, . . . , p, replace the tentries of level 0 in the jth column of the matrix X with thesequence in the jth column of this (t × p) orthogonal LHD.
Lemma A.1. For any constant a and b, define W c = aT c − bSc .
We have W Tc W c = wc I 2c , where wc = (a 2 2c (2c + 1)(2c+1 + 1)/6 −ab(22c + 2c ) + b2 2c ) and I 2c is a (2c × 2c ) identity matrix.
Proof. W Tc W c = (aT Tc − bSTc )(aT c − bSc ) = a 2 T Tc T c − ab(STcT c + T Tc Sc ) + b2 STc Sc . According to the Lemma 1 and theproof of Theorem 1 in Sun, Liu, and Lin (2009), we knowthat STc Sc = 2c I 2c , STc T c + T Tc Sc = (22c + 2c )I 2c and T Tc T c =Therefore,
wehave
W Tc W c =2c (2c + 1)(2c+1 + 1)/6I 2c .
2 c c
c+1
2c
c
2 c
!
(a 2 (2 + 1)(2 + 1)/6 − ab(2 + 2 ) + b 2 )I 2c .

Theorem A.1. The design matrix X generated from the above procedure is an SLHD, in which the whole design as well as each of its
slices are all second-order orthogonal LHDs.

Proof. Since the design matrix in each slice X i (i = 1, . . . , t) is constructed according to Sun, Liu, and Lin (2009), each of them is guaranteed to be a second-order orthogonal LHD. From the above construction
procedure, it is also obvious that each column of the matrix X forms
a permutation of {1 − (1 + n)/2, 2 − (1 + n)/2, . . . , n − (1 + n)/2},
and thus X is an n-run SLHD.
Now we show that the columns in X are orthogonal to each other.
First, it is easy to see that Sc is the matrix whose entries are taken
as 1 or −1 according to the signs of the elements in matrix T c (orequivalently in matrix H c ). When m = 2c+1 , the proposed procedure constructs each X i (i = 1, . . . , t) as a set of p columns from

OPTIMAL SLHDS

Downloaded by [Ohio State University Libraries] at 13:49 19 November 2015(i)(i)T T(i)
L̃ c = (M (i)Tc , −M c ) , where M c = t H c + (i − 0.5 − 0.5t)Sc =M (i)t T c + (i − 0.5 − t)Sc . According to Lemma 1, M (i)Tcc is a diag(i)onal matrix. Therefore, all the columns in each L̃ c are orthogonal
to each other (i = 1, . . . , t), and thus X = (X T1 , . . . , X Tt )T is an orthogonal LHD. When m = 2c+1 + 1, the proposed procedure con(i)structs each X i (i = 1, . . . , t) as a set of p columns from L̃ c =(i)(i)T T(i)(M (i)Tc , 0̃2c , −M c ) , where M c = t T c + (i − 0.5 − 0.5t)Sc and(i)0̃2c corresponds to the original 2c × 1 zero vector but its elements
are relabeled according to Step 3. Using the result in Lemma 1, weM (i)know that M (i)Tcc is a diagonal matrix. Since Step 3 has relabeled the t original zero vectors from X 1 , . . . , X t in such a way thatthey form a (t × p) second-order orthogonal LHD, it is obvious thatX = (X T1 , . . . , X Tt )T is an orthogonal LHD.
Now we only need to prove that the columns in X also satisfies thesecond-order orthogonality property. Observe that the matrix X canalways be arranged into blocks as X = (K T , −K T )T when n is evenor X = (K T , 0p , −K T )T when n is odd. Thus, for any three columns
in X, no matter whether they are distinct or not, they can be denoted as cj = (kTj , −kTj )T when n is even or cj = (kTj , 0, −kTj )T whenn is odd, for j = 1, 2, 3. Then, we have (c1 ⊙ c2 ⊙ c3 )T 1n = (k1 ⊙k2 ⊙ k3 )T 1m − (k1 ⊙ k2 ⊙ k3 )T 1m = 0, where u ⊙ v represents theelement-wise product of u and v, 1m represents the m × 1 vector with allelements equal one, m = n/2 when n is even and m = (n − 1)/2 whenn is odd. Now, it follows immediately that the element-wise square ofeach column and the element-wise product of every two columns in
X are orthogonal to all the columns in the design.
!

Since any 2c+1 -run second-order orthogonal LHD can at most contain 2c columns (Sun, Liu, and Lin 2009), it is easy to see that form = 2c+1 (∀c ≥ 1), the second-order orthogonal SLHDs constructedby the proposed procedure have achieved the maximum 2c number offactors (for any number of slices t ≥ 1). This is a substantial improvement compared to the existing method in Yang et al. (2013), whose
second-order orthogonal SLHDs must satisfy the following two restrictions: (i) the number of slices is a power of two: t = 2r , r = 1, . . . , c;and (ii) each slice contains 22c−r+1 rows but only 2c columns. Obviously, the slice number t in Yang et al. (2013) is much more restrictive,
and unless the design contains exactly 2c slices (r = c), the number ofcolumns these designs can accommodate are much smaller.

ACKNOWLEDGMENTS
The authors thank the editor, an associate editor, and two
referees for their valuable comments and suggestions. This work
was partially completed when the first author was an assistant

487

professor at the Fariborz Maseeh Department of Mathematics
and Statistics in the Portland State University.
[Received December 2012. Revised July 2014.]

REFERENCES
Fang, K. T., Li, R., and Sudjianto, A. (2005), Design and Modeling for Computer
Experiments, London: Chapman and Hall. [479]
Hedayat, A., Sloane, N., and Stufken, J. (1999), Orthogonal Arrays, New York:
Springer Verlag. [482]
Hickernell, F. J. (1998), “A Generalized Discrepancy and Quadrature Error Bound,” Mathematics of Computation, 67, 299–322.
[482]
Jin, R., Chen, W., and Sudjianto, A. (2005), “An Efficient Algorithm for Constructing Optimal Design of Computer Experiments,” Journal of Statistical
Planning and Inference, 134, 268–287. [481,482,483]
Johnson, M. E., Moore, L. M., and Ylvisaker, D. (1990), “Minimax and Maximin Distance Designs,” Journal of Statistical Planning and Inference, 26,
131–148. [482,485]
Jones, B., and Nachtsheim, C. J. (2011), “Efficient Designs with Minimal Aliasing,” Technometrics, 53, 62–71. [482]Joseph, V. R., and Hung, Y. (2008), “Orthogonal-Maximin Latin Hypercube
Designs,” Statistica Sinica, 18, 171–186. [481]Lundy, M., and Mees, A. (1986), “Convergence of an Annealing Algorithm,”
Mathematical Programming, 34, 111–124. [483]
McKay, M. D., Beckman, R. J., and Conover, W. J. (1979), “A Comparison of Three Methods for Selecting Values of Input Variables in
the Analysis of Output From a Computer Code,” Technometrics, 21,
239–245. [479]Morris, M. D., and Mitchell, T. J. (1995), “Exploratory Designs for Computational Experiments,” Journal of Statistical Planning and Inference, 43,
381–402. [480,482]
Owen, A. B. (1992), “Orthogonal Arrays for Computer Experiments, Integration
and Visualization,” Statistica Sinica, 2, 439–452. [482]Qian, P. Z. G. (2012), “Sliced Latin Hypercube Designs,” Journal of the American Statistical Association, 107, 393–399. [479,481,486]
Santner, T. J., Williams, B. J., and Notz, W. (2003), The Design and Analysis ofComputer Experiments, New York: Springer. [479]Sun, F. S., Liu, M. Q., and Lin, D. K. J. (2009), “Construction of Orthogonal Latin Hypercube Designs,” Biometrika, 96,
971–974. [482,486,487]
Tang, B. (1993), “Orthogonal Array-Based Latin Hypercubes,” Journal of theAmerican Statistical Association, 88, 1392–1397. [482]Wu, C. F. J., and Hamada, M. S. (2009), Experiments: Planning, Analysis, andOptimization (2nd ed.), New York: Wiley. [482]Yang, J., Lin, C. D., Qian, P. Z. G., and Lin, D. K. J. (2013), “Construction of Sliced Orthogonal Latin Hypercube Designs,” Statistica Sinica, 23,
1117–1130. [482,487]
Ye, K. Q. (1998), “Orthogonal Column Latin Hypercubes and Their Applicationin Computer Experiments,” Journal of the American Statistical Association,
93, 1430–1439. [482]

TECHNOMETRICS, NOVEMBER 2015, VOL. 57, NO. 4