Latin Hypercubes and Space-filling Designs

C. Devon Lin and Boxin Tang

(Feb 16, 2014)

19.1 Introduction

This chapter discusses a general design approach to planning computer experiments, which seeks design points that fill a bounded design region as uniformly as possible. Such designs are broadly referred to as space-filling designs.

The literature on the design for computer experiments has focused mainly on deterministic computer models; that is, running computer code with the same inputs always produces the same outputs (see Chapter 17). Because of this feature, the three fundamental design principles, randomization, replication and blocking, are irrelevant in computer experiments. The true relationship between the inputs and the responses is unknown and often very complicated. To explore the relationship, one could use traditional regression models. But the most popular are Gaussian process models; see Chapter 17 for details. However, before data are collected, quite often little a priori or background knowledge is available about which model would be appropriate, and designs for computer experiments should facilitate diverse modelling methods. For this purpose, a space-filling design is the best choice. The design region in which to make prediction may be unspecified at the data collection stage. Therefore, it is appropriate to use designs that represent all portions of the design region. When the primary goal of experiments is to make prediction at unsampled points, space-filling designs allow us to build a predictor with better average accuracy.

One most commonly used class of space-filling designs for computer experiments is that of Latin hypercube designs. Such designs, introduced by McKay et al. (1979), do not have repeated runs. Latin hypercube designs have one-dimensional uniformity in that, for each input variable, if its range is divided into the same number of equally-spaced intervals as the number of observations, there is exactly one observation in each interval. However, a random Latin hypercube design may not be a good choice with respect to some optimality criteria such as maximin distance and orthogonality (discussed later). The maximin distance criterion, introduced by Johnson et al. (1990), maximizes the smallest distance between any two design points so that no two design points are too close. Therefore, a maximin distance design spreads out its points evenly over the entire design region. To further enhance the space-filling property for each individual input of a maximin distance design, Morris and Mitchell (1995) proposed the use of maximin Latin hypercube designs.

Many applications involve a large number of input variables. Finding space-filling designs with a limited number of design points that provide a good coverage of the entire high dimensional input space is a very ambitious, if not hopeless, undertaking. A more reasonable approach is to construct designs that are space-filling in the low dimensional projections. Moon et al. (2011) constructed designs that are space-filling in the two-dimensional projections and demonstrated empirically that such designs also perform well in terms of the maximin distance criterion in higher dimensions. Other designs that are space-filling in the low dimensional projections are randomized orthogonal arrays Owen (1992) and orthogonal array-based Latin hypercubes Tang (1993). Another important approach is to construct orthogonal Latin hypercube designs. The basic idea of this approach is that orthogonality can be viewed as a stepping stone to constructing designs that are space-filling in low dimensional projections Bingham et al. (2009).

Originating as popular tools in numerical analysis, low-discrepancy nets, low-discrepancy sequences and uniform designs have also been well recognized as space-filling designs for computer experiments. These designs are chosen to achieve uniformity in the design space based on the discrepancy criteria such as the $L_{p}$ discrepancy (see Section 19.3.2).

As an alternative to the use of space-filling designs, one could choose designs that perform well with respect to some model-dependent criteria such as the minimum integrated mean square error and the maximum entropy Sacks et al. (1989); Shewry and Wynn (1987). One drawback of this approach is that such designs require the prior knowledge of the model. For instance, to be able to construct maximum entropy designs and integrated mean square error optimal designs, one would need the values of the parameters in the correlation function when a Gaussian process is used to model responses. One could also consider a Bayesian approach Leatherman et al. (2014). A detailed account of model-dependent designs can be found in Santner et al. (2003), Fang et al. (2006) and the references therein.

This chapter is organized as follows. Section 19.2 gives a detailed review of Latin hypercube designs, and discusses three important types of Latin hypercube designs (Latin hypercube designs based on measures of distance; orthogonal array-based Latin hypercube designs; orthogonal and nearly orthogonal Latin hypercube designs). Section 19.3 describes other space-filling designs that are not Latin hypercube designs. Concluding remarks are provided in Section 19.4.

19.2 Latin Hypercube Designs

19.2.1 Introduction and examples

A Latin hypercube of $n$ runs for $k$ factors is represented by an $n\times k$ matrix, each column of which is a permutation of $n$ equally spaced levels. For convenience, the $n$ levels are taken to be $-(n-1)/2,-(n-3)/2,\ldots,(n-3)/2,(n-1)/2$ . For example, design $L$ in Table 19.1 is a Latin hypercube of $5$ runs for $3$ factors. Given an $n\times k$ Latin hypercube $\mbox{\boldmath{$L$}}=(l_{ij})$ , a Latin hypercube design $D$ in the design space $[0,1)^{k}$ can be generated and the design matrix of $D$ is an $n\times k$ matrix with the $(i,j)$ th entry being

\displaystyle d_{ij}=\frac{l_{ij}+(n-1)/2+u_{ij}}{n},\ \ \ i=1,\ldots,n,j=1,\ldots,k,

(19.1)

where $u_{ij}$ ’s are independent random numbers from $[0,1)$ . If each $u_{ij}$ in (19.1) is taken to be 0.5, the resulting design $D$ is termed “lattice sample” due to Patterson (1954). For each factor, Latin hypercube designs have exactly one point in each of the $n$ intervals $[0,1/n),[1/n,2/n),\ldots,[(n-1)/n,1)$ . This property is referred to as one-dimensional uniformity. For instance, design $D$ in Table 19.1 is a Latin hypercube design based on the $L$ in the table, and its pairwise plot in Figure 19.1 illustrates the one-dimensional uniformity. When the five points are projected onto each axis, there is exactly one point in each of the five equally-spaced intervals.

Table 19.1: A

5\times 3

Latin hypercube

L

and a Latin hypercube design

D

based on

L

$L$			$D$
2	0	-2	0.9253	0.5117	0.1610
1	-2	-1	0.7621	0.1117	0.3081
-2	2	0	0.1241	0.9878	0.4473
0	-1	2	0.5744	0.3719	0.8270
-1	1	1	0.3181	0.7514	0.6916

Refer to caption — Figure 19.1: The pairwise plot of the Latin hypercube design $D$ in Table 19.1 for the three factors $x_{1},x_{2},x_{3}$

The popularity of Latin hypercube designs was largely attributed to their theoretical justification for the variance reduction in numerical integration. Consider a function $y=f(\mbox{\boldmath{$x$}})$ where $f$ is known, $\mbox{\boldmath{$x$}}=(x_{1},\ldots,x_{k})$ has a uniform distribution in the unit hypercube $[0,1)^{k}$ , and $y\in R$ . (More generally, when $x_{j}$ follows a continuous distribution with a cumulative distribution function $F_{j}$ , then the inputs of $x_{j}$ can be selected via the quantile transformation $F_{j}^{-1}(u_{j})$ where $u_{j}$ follows a uniform distribution in $[0,1)$ ). The expectation of $y$ ,

\displaystyle\mu=\hbox{E}(y),

(19.2)

is of interest. When the expectation $\mu$ cannot be computed explicitly or its derivation is unwieldy, one can resort to approximate methods. Let $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ be a sample of size $n$ . One estimate of $\mu$ in (19.2) is

\displaystyle\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}f(\mbox{\boldmath{$x$}}_{i}).

(19.3)

The approach that takes $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ independently from the uniform distribution in $[0,1)^{k}$ is simple random sampling. McKay et al. (1979) suggested an approach based on a Latin hypercube sample $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ . Denote the estimator $\hat{\mu}$ in (19.3) of $\mu$ under simple random sampling and Latin hypercube sampling by $\hat{\mu}_{srs}$ and $\hat{\mu}_{lhs}$ , respectively. Note that $\hat{\mu}_{srs}$ and $\hat{\mu}_{lhs}$ have the same form but $\hat{\mu}_{srs}$ uses a simple random sample and $\hat{\mu}_{lhs}$ a Latin hypercube sample. Both samples are denoted by $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ , for convenience. McKay et al. (1979) established the following theorem.

Theorem 19.2.1.

If $y=f(\mbox{\boldmath{$x$}})$ is monotonic in each of its input variables, then $\hbox{Var}(\hat{\mu}_{lhs})\leq\hbox{Var}(\hat{\mu}_{srs}).$

Theorem 19.2.1 says that when the monotonicity condition holds, Latin hypercube sampling yields a smaller variance of the sample mean than simple random sampling. Theorem 19.2.2 below Stein (1987) provides some insights into the two methods of sampling.

Theorem 19.2.2.

We have that for $\mbox{\boldmath{$x$}}\in[0,1)^{k}$ ,

\hbox{Var}(\hat{\mu}_{srs})=\frac{1}{n}\hbox{Var}[f(\mbox{\boldmath{$x$}})]

and

\hbox{Var}(\hat{\mu}_{lhs})=\frac{1}{n}\hbox{Var}[f(\mbox{\boldmath{$x$}})]-\frac{1}{n}\sum_{j=1}^{k}\hbox{Var}[f_{j}(x_{j})]+o(\frac{1}{n}),

where $x_{j}$ is the $j$ th input of $x$ , $f_{j}(x_{j})=\hbox{E}[f(\mbox{\boldmath{$x$}})|x_{j}]-\mu$ and $o(\cdot)$ is little o notation.

The term $f_{j}(x_{j})$ in Theorem 19.2.2 is the main effect of the $j$ th input variable. Theorem 19.2.2 tells us that the variance of the sample mean under Latin hypercube sampling is smaller than the counterpart under simple random sampling by an amount contributed by main effects. The extent of the variance reduction depends on the extent to which the function $f$ is additive in the inputs. Asymptotic normality and a central limit theorem of Latin hypercube sampling were established in Stein (1987) and Owen (1992), respectively. A related approach is that of quasi-Monte Carlo methods, which selects design points in a deterministic fashion (see Niederreiter 1992, and Section 19.3.2).

A randomly generated Latin hypercube design does not necessarily perform well with respect to criteria such as those of “space-filling” or “orthogonality”, alluded to in Section 19.1. For example, when projected onto two factors, design points in a random Latin hypercube design may roughly lie on the diagonal as in the plot of $x_{1}$ versus $x_{2}$ in Figure 19.1, leaving a large area in the design space unexplored. In this case, the corresponding two columns in the design matrix are highly correlated. Examples of Latin hypercube designs with desirable properties are maximin Latin hypercube designs, orthogonal-array based Latin hypercube designs, and orthogonal or nearly orthogonal Latin hypercube designs; these will be discussed throughout the chapter.

19.2.2 Latin hypercube designs based on measures of distance

To construct space-filling Latin hypercube designs, one natural approach is to make use of distance criteria. In what follows, we review several measures of distance.

Let $\mbox{\boldmath{$u$}}=(u_{1},\ldots,u_{k})$ and $\mbox{\boldmath{$v$}}=(v_{1},\ldots,v_{k})$ be two design points in the design space $\mbox{ \boldmath${\chi}$}=[0,1]^{k}$ . For $t>0$ , define the inter-point distance between $u$ and $v$ to be

d(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}})=\big{(}\sum_{j=1}^{k}|u_{j}-v_{j}|^{t}\big{)}^{1/t}.

(19.4)

When $t=1$ and $t=2$ , the measure in (19.4) becomes the rectangular and Euclidean distances, respectively. The maximin distance criterion seeks a design $D$ of $n$ points in the design space ${\chi}$ that maximizes the smallest inter-point distance; that is, it maximizes

\min_{\begin{subarray}{c}\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}}\in\mbox{\boldmath{$D$}}\\ {\mbox{\boldmath{$u$}}\neq\mbox{\boldmath{$v$}}}\end{subarray}}d(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}}),

(19.5)

where $d(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}})$ is defined as in (19.4) for any given $t$ . This criterion attempts to place the design points such that no two points are too close to each other.

A slightly different idea is to spread out the design points of a design $D$ in such a way that every point in the design space ${\chi}$ is close to some point in $D$ . This is the minimax distance criterion which seeks a design $D$ of $n$ points in ${\chi}$ that minimizes the maximum distance between an arbitrary point $\mbox{\boldmath{$x$}}\in\mbox{ \boldmath${\chi}$}$ and the design $D$ ; that is, it minimizes

\max_{\mbox{\boldmath{$x$}}\in\mbox{ \boldmath${\chi}$}}d(\mbox{\boldmath{$x$}},\mbox{\boldmath{$D$}}),

where $d(\mbox{\boldmath{$x$}},\mbox{\boldmath{$D$}})$ , representing the distance between $x$ and the closest point in $D$ , is defined as $d(\mbox{\boldmath{$x$}},\mbox{\boldmath{$D$}})=\min_{\mbox{\boldmath{$x$}}_{i}\in\mbox{\boldmath{$D$}}}d(\mbox{\boldmath{$x$}},\mbox{\boldmath{$x$}}_{i})$ and $d(\mbox{\boldmath{$x$}},\mbox{\boldmath{$x$}}_{i})$ is given in (19.4) for any given $t$ .

Audze and Eglais (1977) introduced a distance criterion similar in spirit to the maximin distance criterion by using

\displaystyle\sum_{1\leq i<j\leq n}d(\mbox{\boldmath{$x$}}_{i},\mbox{\boldmath{$x$}}_{j})^{-2},

(19.6)

where $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ are the design points. This criterion of minimizing (19.6) was used by Liefvendahl and Stocki (2006).

Moon et al. (2011) defined a two-dimensional maximin distance criterion. Let the inter-point distance between two design points $\mbox{\boldmath{$u$}}=(u_{1},\ldots,u_{k})$ and $\mbox{\boldmath{$v$}}=(v_{1},\ldots,v_{k})$ projected onto dimensions $h$ and $l$ be

d^{(2)}_{h,l}(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}})=\big{(}|u_{h}-v_{h}|^{t}+|u_{l}-v_{l}|^{t}\big{)}^{1/t},\ \ t>0.

Then the minimum inter-point distance of a design $D$ over all two-dimensional subspaces is

d_{min}^{(2)}=\min_{\begin{subarray}{c}\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}}\in\mbox{\boldmath{$D$}}\\ {\mbox{\boldmath{$u$}}\neq\mbox{\boldmath{$v$}}},h\neq l\end{subarray}}d^{(2)}_{h,l}(\mbox{\boldmath{$u$}},\mbox{\boldmath{$v$}}).

(19.7)

The two-dimensional maximin distance criterion selects a design that maximizes $d_{min}^{(2)}$ in (19.7). Moon et al. (2011) showed by examples that optimal Latin hypercube designs based on this criterion also perform well under the maximin distance criterion (19.5).

Maximin Latin hypercube designs

We now focus on maximin distance criterion. Recall the Gaussian process model in Section 17.4.1,

\displaystyle Y(x)=\mu+Z(x),

(19.8)

where $\mu$ is the unknown but constant mean function, $Z(x)$ is a stationary Gaussian process with mean 0, variance $\sigma^{2}$ , and correlation function $R(\cdot|\theta)$ . A popular choice for the correlation function is the power exponential correlation

\displaystyle R(\mbox{\boldmath{$h$}}|\theta)=\hbox{exp}\big{(}-\theta\sum_{j=1}^{k}|h_{j}|^{p}\big{)},\ \ 0<p\leq 2,

where $h_{j}$ is the $j$ th element of $h$ . Johnson et al. (1990) showed that as the correlation parameter $\theta$ goes to infinity, a maximin design maximizes the determinant of the correlation matrix, where the correlation matrix refers to that of the outputs from running the computer model at the design points. That is, a maximin design is asymptotically D-optimal under the model in (19.8) as the correlations become weak. Thus, a maximin design is also asymptotically optimal with respect to the maximum entropy criterion Shewry and Wynn (1987).

The problem of finding maximin designs is referred to as the maximum facility dispersion problem Erkut (1990) in location theory. It is closely related to the sphere packing problem in the field of discrete and computational geometry Melissen (1997); Conway et al. (1999). The two problems are, however, different as explained in Johnson et al. (1990).

An extended definition of a maximin design was given by Morris and Mitchell (1995). Define a distance list $(d_{1},\ldots,d_{m})$ and an index list $(J_{1},\ldots,J_{m})$ respectively in the following way. The distance list contains the distinct values of inter-point distances, sorted from the smallest to the largest, and $J_{i}$ in the index list is the number of pairs of design points in the design separated by the distance $d_{i}$ , $i=1,\ldots,m$ . Note that $1\leq m\leq{n\choose 2}$ . In Morris and Mitchell (1995), a design is called a maximin design if it sequentially maximizes $d_{i}$ ’s and minimizes $J_{i}$ ’s in the order $d_{1},J_{1},d_{2},J_{2},\ldots,d_{m},J_{m}$ . They further introduced a computationally efficient scalar-value criterion

\displaystyle\phi_{q}=\big{(}\sum_{i=1}^{m}\frac{J_{i}}{d_{i}^{q}}\big{)}^{1/q},

(19.9)

where $q$ is a positive integer. Minimizing $\phi_{q}$ with a large $q$ results in a maximin design. Values of $q$ are chosen depending on the size of the design searched for, ranging from 5 for small designs to 20 for moderate-sized designs to 50 for large designs.

Table 19.2: Algorithms for generating maximin Latin hypercube designs

Article	Algorithm	Criterion
Morris and Mitchell (1995)	simulated annealing	$\phi_{q}^{(a)}$
Ye et al. (2000)	columnwise-pairwise	$\phi_{q}$ and entropy
Jin et al. (2005)	enhanced stochastic	$\phi_{q}$ , entropy and $L_{2}$ discrepancy
	evolutionary algorithm
Liefvendahl and Stocki (2006)	columnwise-pairwise and	maximin and the Audze-Eglais
	genetic algorithm	function^(b)
van Dam et al. (2007)	branch-and-bound	maximin with Euclidean distance
Forrester et al. (2008)	evolutionary operation	$\phi_{q}$
Grosso et al. (2009)	iterated local search heuristics	$\phi_{q}$
Viana et al. (2010)	translational propagation	$\phi_{q}$
Zhu et al. (2011)	successive local enumeration	maximin
Moon et al. (2011)	smart swap algorithm	${d_{min}^{(2)}}^{(c)}$ , $\phi_{q}$
Chen et al. (2013)	particle swarm algorithm	$\phi_{q}$
(a): $\phi_{q}$ is given as in (19.9); (b): the Audze-Eglais function in (19.6); (c) $d_{min}^{(2)}$ is given in (19.7).

Maximin designs tend to place design points toward or on the boundary. For example, Figure 19.2 exhibits a maximin Euclidean distance design and a maximin rectangular distance design, both of seven points. Maximin designs are likely to have clumped projections onto one-dimension. Thus, such designs may not possess desirable one-dimensional uniformity which is guaranteed by Latin hypercube designs. To strike the balance, Morris and Mitchell (1995) examined maximin designs within Latin hypercube designs. Although this idea sounds simple, generating maximin Latin hypercube designs is a challenging task particularly for large designs. The primary approach for obtaining maximin Latin hypercube designs is using the algorithms summarized in Table 19.2. These algorithms search for maximin Latin hypercube designs that have $u_{ij}$ in (19.1) being a constant. For example, Moon et al. (2011) used $u_{ij}=0.5$ , which corresponds to the midpoint of the interval $[(i-1)/n,i/n]$ for $i=1,\ldots,n$ . For detailed descriptions of these algorithms, see the respective references.

Some available implementations of the algorithms in Table 19.2 include the Matlab code provided in Viana et al. (2010), the function maximinLHS in the R package lhs Carnell (2009), and the function lhsdesign in the Matlab statistics toolbox. The function in R uses random $u_{ij}$ ’s in (19.1) while the function in Matlab allows both random $u_{ij}$ ’s and $u_{ij}=0.5$ . It should be noted, however, that these designs are approximate maximin Latin hypercube designs. No theoretical method is available to construct exact maximin Latin hypercube designs of flexible run sizes except Tang (1994) and van Dam et al. (2007). These two articles provided methods for constructing exact maximin Latin hypercube designs of certain run sizes and numbers of input variables. Tang (1994) constructed a Latin hypercube based on a single replicate full factorial design (see Chapter 1 and also Wu and Hamada 2011, Chapter 4) and showed that the constructed Latin hypercube has the same rectangular distance as the single replicate full factorial design, where the rectangular distance of a design is given by (19.5) with $t=1$ . van Dam et al. (2007) constructed two-dimensional maximin Latin hypercubes with the distance measures with $t=1$ and $t=\infty$ in (19.4). For the Euclidean distance measure with $t=2$ , van Dam et al. (2007) used the branch-and-bound algorithm to find maximin Latin hypercube designs with $n\leq 70$ .

19.2.3 Orthogonal array-based Latin hypercube designs

Tang (1993) proposed orthogonal array-based Latin hypercube designs, also known as U designs, which guarantee multi-dimensional space-filling. Recall the definition of an $s$ -level orthogonal array (OA) of $n$ runs, $k$ factors and strength $r$ , denoted by ${\mbox{OA}}(n,s^{k},r)$ in Chapter 9. The $s$ levels are taken to be $1,2,\ldots,s$ in this chapter. By the definition of orthogonal arrays, a Latin hypercube of $n$ runs for $k$ factors is an ${\mbox{OA}}(n,n^{k},1)$ .

Table 19.3: An

{\mbox{OA}}(9,3^{4},2)

and a corresponding OA-based Latin hypercube

${\mbox{OA}}(9,3^{4},2)$				$L$
1	1	1	1	-2	-2	-4	-2
1	2	2	3	-4	0	1	2
1	3	3	2	-3	4	2	1
2	1	2	2	-1	-4	-1	-1
2	2	3	1	1	-1	4	-3
2	3	1	3	0	2	-3	4
3	1	3	3	3	-3	3	3
3	2	1	2	2	1	-2	0
3	3	2	1	4	3	0	-4

The construction of OA-based Latin hypercubes in Tang (1993) works as follows. Let $A$ be an ${\mbox{OA}}(n,s^{k},r)$ . For each column of $A$ and $m=1,\ldots,s$ , replace the $n/s$ positions with entry $m$ by a random permutation of $(m-1)n/s+1,(m-1)n/s+2,\ldots,mn/s$ . Denote the design after the above replacement procedure by $\mbox{\boldmath{$A$}}^{\prime}$ . In our notation, an OA-based Latin hypercube is given by $\mbox{\boldmath{$L$}}=\mbox{\boldmath{$A$}}^{\prime}-(n+1)\mbox{\boldmath{$J$}}/2$ , where $J$ is an $n\times k$ matrix of all 1’s. An OA-based Latin hypercube design in the design space $[0,1)^{k}$ can be generated via (19.1). In addition to the one-dimensional uniformity, an ${\mbox{OA}}(n,s^{k},r)$ -based Latin hypercube has the $r$ -dimensional projection property that when projected onto any $r$ columns, it has exactly $\lambda=n/s^{r}$ points in each of the $s^{r}$ cells $\mathcal{P}^{r}$ where $\mathcal{P}=\{[0,1/s],[1/s,2/s),\ldots,[1-1/s,1)\}$ . Example 19.2.1 illustrates this feature of an ${\mbox{OA}}(9,3^{4},2)$ -based Latin hypercube.

Example 19.2.1.

Table 19.3 displays an OA-based Latin hypercube $L$ based on the orthogonal array ${\mbox{OA}}(9,3^{4},2)$ in the table. Figure 19.3 depicts the pairwise plot of this Latin hypercube. In each subplot, there is exactly one point in each of nine dot-dash line boxes.

A generalization of OA-based Latin hypercubes using asymmetrical orthogonal arrays (see Chapter 9) can be readily made. For example, Figure 19.4 displays a Latin hypercube design based on an asymmetrical orthogonal array of six runs for two factors with three levels in the first factor and two levels in the second factor. Note that each of the six cells separated by dot-dash lines contains exactly one point. By contrast, in the six-point randomly generated Latin hypercube design displayed in Figure 19.4, two out of six such cells do not contain any point.

Orthogonal arrays have been used directly to carry out computer experiments; see, for example, Joseph et al. (2008). Compared with orthogonal arrays, OA-based Latin hypercubes are more favorable for computer experiments. When projected onto lower dimensions, the design points in orthogonal arrays often have replicates. This is undesirable at the early stages of experimentation when relatively few factors are likely to be important.

The construction of OA-based Latin hypercubes depends on the existence of orthogonal arrays. For the existence results of orthogonal arrays, see, for example, Hedayat et al. (1999) and Mukerjee and Wu (2006). A library of orthogonal arrays is freely available on the N.J.A. Sloane website and the MktEx macro using the software SAS Kuhfeld (2009). It should be noted that, for certain given run sizes and numbers of factors, orthogonal arrays of different numbers of levels and different strengths may be available. For instance, an ${\mbox{OA}}(16,4^{5},2)$ , an ${\mbox{OA}}(16,2^{5},4)$ and an ${\mbox{OA}}(16,2^{5},2)$ all produce OA-based Latin hypercubes of 16 runs for 5 factors. However, orthogonal arrays with more levels and/or higher strength are preferred because they lead to designs with better projection space-filling properties.

Given an orthogonal array, the construction of Tang (1993) can produce many OA-based Latin hypercubes. There arises the problem of choosing a preferable OA-based Latin hypercube. Leary et al. (2003) presented an algorithm for finding optimal OA-based Latin hypercubes that minimize (19.6) using the Euclidean distance between design points. The optimization was performed via the simulated annealing algorithm Morris and Mitchell (1995) and the columnwise-pairwise algorithm Li and Wu (1997).

Recall the problem of estimating the mean $\mu$ in (19.2) of a known function $y=f(\mbox{\boldmath{$x$}})$ using a design with $n$ points $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ in Section 19.2.1. Latin hypercube sampling stratifies all univariate margins simultaneously and thus achieves a variance reduction compared with simple random sampling, as quantified in Theorem 19.2.2. Theorem 19.2.3 below from Tang (1993) shows that a further variance reduction is achieved by OA-based Latin hypercube sampling.

Theorem 19.2.3.

Suppose that $f$ is continuous on $[0,1]^{k}$ . Let $\hat{\mu}_{oalhs}$ denote the $\hat{\mu}$ in (19.3) with a randomly selected OA-based Latin hypercube design of $n$ points. Then we have that

\hbox{Var}(\hat{\mu}_{oalhs})=\frac{1}{n}\hbox{Var}[f(\mbox{\boldmath{$x$}})]-\frac{1}{n}\sum_{j=1}^{k}\hbox{Var}[f_{j}(x_{j})]-\frac{1}{n}\sum_{i<j}^{k}\hbox{Var}[f_{ij}(x_{i},x_{j})]+o(\frac{1}{n}),

where $x_{j}$ is the $j$ th input of $x$ , $f_{j}(x_{j})=\hbox{E}[f(\mbox{\boldmath{$x$}})|x_{j}]-\mu$ , and $f_{ij}(x_{i},x_{j})=\hbox{E}[f(\mbox{\boldmath{$x$}})|x_{i},x_{j}]-\mu-f_{i}(x_{i})-f_{j}(x_{j})$ .

To better understand this result, we write

f(\mbox{\boldmath{$x$}})=\mu+\sum_{j=1}^{k}f_{j}(x_{j})+\sum_{i<j}^{k}f_{ij}(x_{i},x_{j})+r(\mbox{\boldmath{$x$}}),

where the terms on the right side of the equation are uncorrelated with each other. Thus, the variance decomposition of the function $f$ is

\hbox{Var}[f(\mbox{\boldmath{$x$}})]=\sum_{j=1}^{k}\hbox{Var}[f_{j}(x_{j})]+\sum_{i<j}^{k}\hbox{Var}[f_{ij}(x_{i},x_{j})]+\hbox{Var}[r(\mbox{\boldmath{$x$}})].

We see that Latin hypercube sampling achieves a variance reduction by removing the variances of the main effects $f_{j}(x_{j})$ from $\hbox{Var}[f(\mbox{\boldmath{$x$}})]$ , and OA-based Latin hypercube sampling further removes the variances of the interactions $f_{ij}(x_{i},x_{j})$ .

We conclude this section by mentioning that randomized orthogonal arrays proposed by Owen (1992) also enjoy good space-filling properties in the low-dimensional projections. Results similar to Theorem 19.2.3 are given in Owen (1992).

19.2.4 Orthogonal and nearly orthogonal Latin hypercube designs

This section discusses the properties and constructions of Latin hypercube designs that have zero or small column-wise correlations in all two-dimensional projections. Such designs are called orthogonal and nearly orthogonal Latin hypercube designs. Orthogonal Latin hypercube designs are directly useful in fitting data using main effects models because they allow uncorrelated estimates of linear main effects. Another rationale of obtaining orthogonal or nearly orthogonal Latin hypercube designs is that they may not be space-filling, but space-filling designs should be orthogonal or nearly orthogonal. Thus we can search for space-filling designs within the class of orthogonal or nearly orthogonal Latin hypercube designs. Other justifications are given in Iman and Conover (1982), Owen (1994), Tang (1998), Joseph and Hung (2008), among others.

Extensive research has been carried out on the construction of orthogonal or nearly orthogonal Latin hypercube designs. Ye (1998) initiated this line of research and constructed orthogonal Latin hypercubes with $n=2^{m}$ or $2^{m}+1$ runs and $k=2m-2$ factors where $m\geq 2$ . Ye’s construction was extended by Cioppa and Lucas (2007) to obtain more columns for given run sizes. Steinberg and Lin (2006) constructed orthogonal Latin hypercubes of the run sizes $n=2^{2^{m}}$ by rotating groups of factors in two-level $2^{2^{m}}$ -run regular fractional factorial designs. This idea was generalized by Pang et al. (2009) who constructed orthogonal Latin hypercubes of $p^{2^{m}}$ runs and up to $(p^{2^{m}}-1)/(p-1)$ factors by rotating groups of factors in $p$ -level $p^{2^{m}}$ -run regular factorial designs, where $p$ is a prime. Lin (2008) obtained orthogonal Latin hypercube designs of small run sizes ( $n\leq 20$ ) through an algorithm that adds columns sequentially to an existing design. More recently, various methods have been proposed to construct orthogonal Latin hypercubes of more flexible run sizes and with large factor-to-run-size ratios. Here we review constructions in Lin et al. (2009), Sun et al. (2009), and Lin et al. (2010). These methods are general, simple to implement, and cover the results in Table 19.10. For other methods, see Georgiou (2009), Sun et al. (2010), and Yang and Liu (2012).

We first give some notation and background. A vector $\mbox{\boldmath{$a$}}=(a_{1},\ldots,a_{n})$ is said to be balanced if its distinct values have equal frequency. For an $n_{1}\times k_{1}$ matrix $A$ and an $n_{2}\times k_{2}$ matrix $B$ , their Kronecker product $\mbox{\boldmath{$A$}}\otimes\mbox{\boldmath{$B$}}$ is the $(n_{1}n_{2})\times(k_{1}k_{2})$ matrix

$\mbox{\boldmath{$A$}}\otimes\mbox{\boldmath{$B$}}=\left[\begin{array}[]{cccc}a_{11}\mbox{\boldmath{$B$}}&a_{12}\mbox{\boldmath{$B$}}&\ldots&a_{1k_{1}}\mbox{\boldmath{$B$}}\\ a_{21}\mbox{\boldmath{$B$}}&a_{22}\mbox{\boldmath{$B$}}&\ldots&a_{2k_{1}}\mbox{\boldmath{$B$}}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n_{1}1}\mbox{\boldmath{$B$}}&a_{n_{1}2}\mbox{\boldmath{$B$}}&\ldots&a_{n_{1}k_{1}}\mbox{\boldmath{$B$}}\end{array}\right]$

with $a_{ij}\mbox{\boldmath{$B$}}$ itself being an $n_{2}\times k_{2}$ matrix. For an $n\times k$ design or matrix $\mbox{\boldmath{$D$}}=(d_{ij})$ , define its correlation matrix to be a $k\times k$ matrix $\mbox{\boldmath{$R$}}(\mbox{\boldmath{$D$}})=(r_{ij})$ with

r_{ij}=\frac{\sum_{m=1}^{n}(d_{mi}-\bar{d}_{i})(d_{mj}-\bar{d}_{j})}{\sqrt{\sum_{m}(d_{mi}-\bar{d}_{i})^{2}\sum_{m}(d_{mj}-\bar{d}_{j})^{2}}}

(19.10)

representing the correlation between the $i$ th and $j$ th columns of $D$ , where $\bar{d}_{i}=n^{-1}\sum_{m=1}^{n}d_{mi}$ and $\bar{d}_{j}=n^{-1}\sum_{m=1}^{n}d_{mj}$ . A design or matrix $D$ is column-orthogonal if $\mbox{\boldmath{$R$}}(\mbox{\boldmath{$D$}})$ is an identity matrix. A design or matrix $\mbox{\boldmath{$D$}}=(d_{ij})$ is orthogonal if it is balanced and column-orthogonal.

To assess near orthogonality of design $D$ , Bingham et al. (2009) introduced two measures, the maximum correlation $\rho_{M}(\mbox{\boldmath{$D$}})=\hbox{max}_{i<j}|r_{ij}|$ and the average squared correlation $\rho^{2}_{ave}(\mbox{\boldmath{$D$}})=\sum_{i<j}r^{2}_{ij}/[(k(k-1)/2]$ , where $r_{ij}$ is defined as in (19.10). Smaller values of $\rho_{M}(\mbox{\boldmath{$D$}})$ and $\rho^{2}_{ave}(\mbox{\boldmath{$D$}})$ imply near orthogonality. Obviously, if $\rho_{M}(\mbox{\boldmath{$D$}})=0$ or $\rho^{2}_{ave}(\mbox{\boldmath{$D$}})=0$ , then an orthogonal design is obtained. For a concise presentation, we use ${\mbox{OLH}}(n,k)$ to denote an orthogonal Latin hypercube of $n$ runs for $k$ factors. Lin et al. (2010) established the following theorem on the existence of orthogonal Latin hypercubes.

Theorem 19.2.4.

There exists an orthogonal Latin hypercube of $n\geq 4$ runs with more than one factor if and only if $n\neq 4m+2$ where $m$ is an integer.

Theorem 19.2.4 says that a Latin hypercube of run size 2, 3, 6, 10, 14, $\ldots$ cannot be orthogonal. For smaller run sizes, this can be readily verified by exhaustive computer search. When orthogonal Latin hypercubes of certain run sizes exist, we want to construct such designs with as many columns as possible. We review three general construction methods. To generate design points in the region $[0,1]^{k}$ from a Latin hypercube, one can use (19.1) with $u_{ij}=0.5$ , which corresponds to the midpoints of the cells.

A construction method based on an orthogonal array and a small orthogonal Latin hypercube

Lin et al. (2009) constructed a large orthogonal, or nearly orthogonal, Latin hypercube by coupling an orthogonal array of index unity with a small orthogonal, or nearly orthogonal, Latin hypercube. Let $B$ be an $n\times q$ Latin hypercube, where as in Section 19.2.1, the levels are $-(n-1)/2,-(n-3)/2,\ldots,(n-3)/2,(n-1)/2$ . Then the elements in every column of $B$ add up to zero while the sum of squares of these elements is $n(n^{2}-1)/12$ . Thus the correlation matrix whose elements are defined as in (19.10) is

\mbox{\boldmath{$R$}}(\mbox{\boldmath{$B$}})=\big{[}\frac{1}{12}n(n^{2}-1)\big{]}^{-1}\mbox{\boldmath{$B$}}^{\prime}\mbox{\boldmath{$B$}}.

(19.11)

Let $A$ be an orthogonal array ${\mbox{OA}}(n^{2},n^{2f},2)$ . The construction proposed by Lin et al. (2009) proceeds as follows.

Step I: Let $b_{ij}$ be the $(i,j)$ th element of $B$ introduced above. For $1\leq j\leq q$ , obtain an $n^{2}\times(2f)$ matrix $\mbox{\boldmath{$A$}}_{j}$ from $A$ by replacing the symbols $1,2,\ldots,n$ in the latter by $b_{1j},b_{2j},\ldots,b_{nj}$ respectively, and then partition $\mbox{\boldmath{$A$}}_{j}$ as $\mbox{\boldmath{$A$}}_{j}$ = [ $\mbox{\boldmath{$A$}}_{j1},\ldots,\mbox{\boldmath{$A$}}_{jf}$ ], where each of $\mbox{\boldmath{$A$}}_{j1},\ldots,\mbox{\boldmath{$A$}}_{jf}$ has two columns.

Step II: For $1\leq j\leq q$ , obtain the $n^{2}\times(2f)$ matrix $\mbox{\boldmath{$L$}}_{j}$ = [ $\mbox{\boldmath{$A$}}_{j1}\mbox{\boldmath{$V$}},\ldots,\mbox{\boldmath{$A$}}_{jf}\mbox{\boldmath{$V$}}$ ] , where

${\mbox{\boldmath{$V$}}=\left[\begin{array}[]{rr}1&-n\\ n&1\\ \end{array}\right].}$

Step III: Obtain the matrix $L$ = [ $\mbox{\boldmath{$L$}}_{1},\ldots,\mbox{\boldmath{$L$}}_{q}$ ], of order $N\times k$ , where $N=n^{2}$ and $k=2qf$ .

For $q=1$ , the above construction is equivalent to that in Pang et al. (2009). However, by Theorem 19.2.4, we have $q\geq 2$ when $n$ is not equal to 3 or $4m+2$ for any non-negative integer $m$ . Thus, the above method provides orthogonal or nearly orthogonal Latin hypercubes with an appreciably larger number of factors as compared to the method in Pang et al. (2009). For example, Pang et al. (2009) obtained ${\mbox{OLH}}(25,6)$ , ${\mbox{OLH}}(49,8)$ , ${\mbox{OLH}}(81,40)$ , ${\mbox{OLH}}(121,12)$ and ${\mbox{OLH}}(169,14)$ while the above construction produces ${\mbox{OLH}}(25,12)$ , ${\mbox{OLH}}(49,24)$ , ${\mbox{OLH}}(81,50)$ , ${\mbox{OLH}}(121,84)$ and ${\mbox{OLH}}(169,84)$ .

Theorem 19.2.5 below shows how the correlation matrix of the large Latin hypercube $L$ depends on that of the small Latin hypercube $B$ .

Theorem 19.2.5.

For the matrix $L$ constructed from $B$ in the above steps I, II and III, we have
(i) the matrix $L$ is a Latin hypercube, and
(ii) the correlation matrix of $L$ is given by

\mbox{\boldmath{$R$}}(\mbox{\boldmath{$L$}})=\mbox{\boldmath{$R$}}(\mbox{\boldmath{$B$}})\otimes\mbox{\boldmath{$I$}}_{2f},

where $\mbox{\boldmath{$R$}}(\mbox{\boldmath{$B$}})$ , given in (19.11), is the correlation matrix of a Latin hypercube $B$ , $\mbox{\boldmath{$I$}}_{2f}$ is the identity matrix of order $2f$ and $\otimes$ denotes Kronecker product.

Corollary 19.2.1.

If $B$ is an orthogonal Latin hypercube, then so is $L$ . In general, the maximum correlation and average squared correlation of $L$ are given by

\rho_{M}(\mbox{\boldmath{$L$}})=\rho_{M}(\mbox{\boldmath{$B$}})\ \ \hbox{ and }\ \ \rho_{ave}^{2}(\mbox{\boldmath{$L$}})=\frac{q-1}{2qf-1}\rho_{ave}^{2}(\mbox{\boldmath{$B$}}).

Corollary 19.2.1 reveals that the large Latin hypercube $L$ inherits the exact or near orthogonality of the small Latin hypercube $B$ . As a result, the effort for searching for large orthogonal or nearly orthogonal Latin hypercubes can be focused on finding small orthogonal or nearly orthogonal Latin hypercubes which are easier to obtain via some general efficient robust optimization algorithms such as simulated annealing and genetic algorithms, by minimizing $\rho_{ave}^{2}$ or $\rho_{M}$ .

Example 19.2.2 below illustrates the actual construction of some orthogonal Latin hypercubes using the method of Lin et al. (2009). Example 19.2.3 is devoted to the construction of a nearly orthogonal Latin hypercube.

Example 19.2.2.

Let $n$ be a prime or prime power for which an ${\mbox{OA}}(n^{2},n^{n+1},2)$ exists Hedayat et al. (1999). For instance, consider $n=5,7,8,9,11$ . Now if we take $B$ to be an ${\mbox{OLH}}(5,2)$ , an ${\mbox{OLH}}(7,3)$ , an ${\mbox{OLH}}(8,4)$ , an ${\mbox{OLH}}(9,5)$ , or an ${\mbox{OLH}}(11,7)$ , as displayed in Table 19.4 and take $A$ respectively to be an ${\mbox{OA}}(25,5^{6},2)$ , an ${\mbox{OA}}(49,7^{8},2)$ , an ${\mbox{OA}}(64,8^{8},2)$ , an ${\mbox{OA}}(81,9^{10},2)$ , or an ${\mbox{OA}}(121,11^{12},2)$ , then the construction described in this section provides an ${\mbox{OLH}}(25,12)$ , an ${\mbox{OLH}}(49,24)$ , an ${\mbox{OLH}}(64,32)$ , an ${\mbox{OLH}}(81,50)$ , or an ${\mbox{OLH}}(121,84)$ , respectively.

Table 19.4: Orthogonal Latin hypercubes OLH(

5,2

), OLH(

7,3

), OLH(

8,4

), OLH(

9,5

) and OLH(

11,7

)

${\mbox{OLH}}(5,2)$					${\mbox{OLH}}(7,3)$						${\mbox{OLH}}(8,4)$
1	-2				-3	3	2				0.5	-1.5	3.5	2.5
2	1				-2	0	-3				1.5	0.5	2.5	-3.5
0	0				-1	-2	-1				2.5	-3.5	-1.5	-0.5
-1	2				0	-3	1				3.5	2.5	-0.5	1.5
-2	-1				1	-1	3				-3.5	-2.5	0.5	-1.5
					2	1	-2				-2.5	3.5	1.5	0.5
					3	2	0				-1.5	-0.5	-2.5	3.5
											-0.5	1.5	-3.5	-2.5
${\mbox{OLH}}(9,5)$								${\mbox{OLH}}(11,7)$
-4	-2	0	-3	3				-5	-4	-5	-5	-3	0	0
-3	4	2	1	-2				-4	2	-1	3	4	5	4
-2	-3	-4	-1	-3				-3	-2	4	5	-4	-2	-1
-1	3	-2	3	4				-2	3	-3	4	1	-4	-2
0	-4	4	4	0				-1	4	2	-4	3	2	-4
1	2	-1	0	-4				0	-5	5	-2	5	-3	2
2	0	3	-2	-1				1	5	3	-3	-5	-1	5
3	1	1	-4	2				2	-1	1	1	-2	3	-5
4	-1	-3	2	1				3	0	0	-1	0	1	-3
								4	1	-4	0	2	-5	1
								5	-3	-2	2	-1	4	3

Example 19.2.3.

Through a computer search, Lin et al. (2009) found a nearly orthogonal Latin hypercube with 13 rows and 12 columns, as given in Table 19.5. This Latin hypercube has $\rho_{ave}=0.0222$ and $\rho_{M}=0.0495$ . Together with an ${\mbox{OA}}(13^{2},13^{14},2)$ , the above procedure provides a nearly orthogonal Latin hypercube of 169 runs and 168 factors with $\rho_{ave}=0.0057$ and $\rho_{M}=0.0495$ , by Corollary 19.2.1.

Table 19.5: A nearly orthogonal Latin hypercube with 13 rows and 12 columns

\begin{array}[]{rrrrrrrrrrrr}-6&-6&-5&-4&-5&-2&2&1&-3&-2&-1&-2\\ -5&5&3&-5&3&4&-6&0&-4&1&-3&-1\\ -4&2&-4&1&2&6&5&-5&6&0&1&1\\ -3&1&2&4&-6&1&-2&6&2&3&2&6\\ -2&-2&6&-3&6&-5&3&4&4&-3&3&0\\ -1&-5&4&6&1&-1&0&-4&0&6&-5&-3\\ 0&6&0&3&-4&-6&-3&-3&3&-5&0&-4\\ 1&0&-3&5&5&0&1&2&-5&-6&-4&5\\ 2&-1&-6&0&4&-4&-5&-2&-1&5&6&2\\ 3&4&1&2&-1&2&6&3&-6&2&5&-6\\ 4&-4&5&-2&-3&3&-1&-6&-2&-4&4&3\\ 5&3&-1&-6&-2&-3&4&-1&1&4&-6&4\\ 6&-3&-2&-1&0&5&-4&5&5&-1&-2&-5\end{array}

Before ending this section, we comment on the projection space-filling property of Latin hypercubes built above using a Latin hypercube $B$ and an orthogonal array $A$ . Any pair of columns obtained using different columns of $A$ retains the two-dimensional projection property of $A$ . When projected to those pairs of columns associated with the same column of $A$ , the design points form clusters and these clusters are spread out as the corresponding two columns of $B$ .

A recursive construction method

Orthogonal Latin hypercubes allow uncorrelated estimates of main effects in a main-effect regression model. Sun et al. (2009) extended the concept of orthogonal Latin hypercubes for second-order polynomial models.

For a design $D$ with columns $\mbox{\boldmath{$d$}}_{1},\ldots,\mbox{\boldmath{$d$}}_{k}$ , let $\tilde{\mbox{\boldmath{$D$}}}$ be the $n\times[k(k+1)/2]$ matrix whose columns consist of all possible products $\mbox{\boldmath{$d$}}_{i}\odot\mbox{\boldmath{$d$}}_{j}$ , where $\odot$ denotes the element-wise product of vectors $\mbox{\boldmath{$d$}}_{i}$ and $\mbox{\boldmath{$d$}}_{j}$ , $i=1,\ldots,k$ , $j=1,\ldots,k$ and $i\leq j$ . Define the correlation matrix between $D$ and $\tilde{\mbox{\boldmath{$D$}}}$ to be

\displaystyle\mbox{\boldmath{$R$}}(\mbox{\boldmath{$D$}},\tilde{\mbox{\boldmath{$D$}}})=\left(\begin{array}[]{rrrr}r_{11}&r_{12}&\ldots&r_{1q}\\ r_{21}&r_{22}&\ldots&r_{2q}\\ \vdots&\vdots&\ddots&\vdots\\ r_{k1}&r_{k2}&\ldots&r_{kq}\\ \end{array}\right),

(19.16)

where $q=k(k+1)/2$ and $r_{ij}$ is the correlation between the $i$ th column of $D$ and the $j$ th column of $\tilde{\mbox{\boldmath{$D$}}}$ . A second-order orthogonal Latin hypercube $D$ has the properties: (a) the correlation matrix $\mbox{\boldmath{$R$}}(\mbox{\boldmath{$D$}})$ is an identity matrix, and (b) $\mbox{\boldmath{$R$}}(\mbox{\boldmath{$D$}},\tilde{\mbox{\boldmath{$D$}}})$ in (19.16) is a zero matrix.

Sun et al. (2009) proposed the following procedure for constructing second-order orthogonal Latin hypercubes of $2^{c+1}+1$ runs in $2^{c}$ factors for any integer $c\geq 1$ . Throughout this section, let $\mbox{\boldmath{$X$}}^{*}$ represent the matrix obtained by switching the signs in the top half of the matrix $X$ with an even number of rows.

Step I: For $c=1$ , let

\mbox{\boldmath{$S$}}_{1}=\left(\begin{array}[]{rr}1&1\\ 1&-1\end{array}\right)\ \mbox{\ and\ }\ \ \mbox{\boldmath{$T$}}_{1}=\left(\begin{array}[]{rr}1&2\\ 2&-1\end{array}\right).

Step II: For an integer $c\geq 2$ , define

\mbox{\boldmath{$S$}}_{c}=\left(\begin{array}[]{cc}\mbox{\boldmath{$S$}}_{c-1}&-\mbox{\boldmath{$S$}}^{*}_{c-1}\\ \mbox{\boldmath{$S$}}_{c-1}&\mbox{\boldmath{$S$}}^{*}_{c-1}\end{array}\right)\ \mbox{\ and\ }\ \mbox{\boldmath{$T$}}_{c}=\left(\begin{array}[]{cc}\mbox{\boldmath{$T$}}_{c-1}&-(\mbox{\boldmath{$T$}}^{*}_{c-1}+2^{c-1}\mbox{\boldmath{$S$}}^{*}_{c-1})\\ \mbox{\boldmath{$T$}}_{c-1}+2^{c-1}\mbox{\boldmath{$S$}}_{c-1}&\mbox{\boldmath{$T$}}^{*}_{c-1}\end{array}\right).

Step III: Obtain a $(2^{c+1}+1)\times 2^{c}$ Latin hypercube $\mbox{\boldmath{$L$}}_{c}$ as

\mbox{\boldmath{$L$}}_{c}=\left(\begin{array}[]{r}\mbox{\boldmath{$T$}}_{c}\\ \textbf{0}_{2^{c}}\\ -\mbox{\boldmath{$T$}}_{c}\\ \end{array}\right),

where $\textbf{0}_{2^{c}}$ denotes a zero row vector of length $2^{c}$ .

Example 19.2.4.

A second-order orthogonal Latin hypercube of 17 runs for 8 factors constructed using the above procedure is given by

{\left(\begin{array}[]{rrrrrrrr}1&2&3&4&5&6&7&8\\[-6.0pt] 2&-1&-4&3&6&-5&-8&7\\[-6.0pt] 3&4&-1&-2&-7&-8&5&6\\[-6.0pt] 4&-3&2&-1&-8&7&-6&5\\[-6.0pt] 5&6&7&8&-1&-2&-3&-4\\[-6.0pt] 6&-5&-8&7&-2&1&4&-3\\[-6.0pt] 7&8&-5&-6&3&4&-1&-2\\[-6.0pt] 8&-7&6&-5&4&-3&2&-1\\[-6.0pt] 0&0&0&0&0&0&0&0\\[-6.0pt] -1&-2&-3&-4&-5&-6&-7&-8\\[-6.0pt] -2&1&4&-3&-6&5&8&-7\\[-6.0pt] -3&-4&1&2&7&8&-5&-6\\[-6.0pt] -4&3&-2&1&8&-7&6&-5\\[-6.0pt] -5&-6&-7&-8&1&2&3&4\\[-6.0pt] -6&5&8&-7&2&-1&-4&3\\[-6.0pt] -7&-8&5&6&-3&-4&1&2\\[-6.0pt] -8&7&-6&5&-4&3&-2&1\\ \end{array}\right).}

Sun et al. (2010) further constructed second-order orthogonal Latin hypercubes of $2^{c+1}$ runs in $2^{c}$ factors by modifying Step III in the above procedure. In Step III, let $\mbox{\boldmath{$H$}}_{c}=\mbox{\boldmath{$T$}}_{c}-\mbox{\boldmath{$S$}}_{c}/2$ and obtain $\mbox{\boldmath{$L$}}_{c}$ as

\mbox{\boldmath{$L$}}_{c}=\left(\begin{array}[]{r}\mbox{\boldmath{$H$}}_{c}\\ -\mbox{\boldmath{$H$}}_{c}\\ \end{array}\right).

Then $\mbox{\boldmath{$L$}}_{c}$ is a second-order orthogonal Latin hypercube of $2^{c+1}$ runs in $2^{c}$ factors.

A construction method based on small orthogonal designs and small orthogonal Latin hypercubes

This section reviews the construction from Lin et al. (2010) for constructing orthogonal and nearly orthogonal Latin hypercubes. All the proofs can be found in Lin et al. (2010). Let $\mbox{\boldmath{$A$}}=(a_{ij})$ be an $n_{1}\times k_{1}$ matrix with entries $a_{ij}=\pm 1$ , $\mbox{\boldmath{$B$}}=(b_{ij})$ be an ${n_{2}\times k_{2}}$ Latin hypercube, $\mbox{\boldmath{$E$}}=(e_{ij})$ be an ${n_{1}\times k_{1}}$ Latin hypercube, and $\mbox{\boldmath{$F$}}=(f_{ij})$ be an ${n_{2}\times k_{2}}$ matrix with entries $d_{ij}=\pm 1$ . Lin et al. (2010) construct designs via

\mbox{\boldmath{$L$}}=\mbox{\boldmath{$A$}}\otimes\mbox{\boldmath{$B$}}+n_{2}\mbox{\boldmath{$E$}}\otimes\mbox{\boldmath{$F$}}.

(19.17)

The resulting design $L$ in (19.17) has $n=n_{1}n_{2}$ runs and $k=k_{1}k_{2}$ factors, and becomes an orthogonal Latin hypercube, if certain conditions on $\mbox{\boldmath{$A$}},\mbox{\boldmath{$B$}},\mbox{\boldmath{$E$}},\mbox{\boldmath{$F$}}$ are met.

Theorem 19.2.6.

Design $L$ in (19.17) is an orthogonal Latin hypercube if
(i) $A$ and $F$ are column-orthogonal matrices of $\pm 1$ ,
(ii) $B$ and $E$ are orthogonal Latin hypercubes,
(iii) at least one of the two, $\mbox{\boldmath{$A$}}^{\prime}\mbox{\boldmath{$E$}}=\textbf{0}$ and $\mbox{\boldmath{$B$}}^{\prime}\mbox{\boldmath{$F$}}=\textbf{0}$ , is true, and
(iv) at least one of the following two conditions is true:
(a) $A$ and $E$ satisfy that for any $i$ , if $p_{1}$ and $p_{2}$ are such that $e_{p_{1}i}=-e_{p_{2}i}$ , then $a_{p_{1}i}=a_{p_{2}i}$ ;
(b) $B$ and $F$ satisfy that for any $j$ , if $q_{1}$ and $q_{2}$ are such that $b_{q_{1}j}=-b_{q_{2}j}$ , then $f_{q_{1}j}=f_{q_{2}j}$ .

Condition (iv) in Theorem 19.2.6 is needed for $L$ to be a Latin hypercube. To make $L$ orthogonal, conditions (i), (ii) and (iii) are necessary. Choices for $A$ and $F$ include Hadamard matrices and orthogonal arrays with levels $\pm 1$ (see Chapter 9). Because of the orthogonality of $A$ and $F$ , $n_{1}$ and $n_{2}$ must be equal to two or multiples of four (Dey and Mukerjee 1999, p. 33). Theorem 19.2.6 requires designs $B$ and $E$ to be orthogonal Latin hypercubes. All known orthogonal Latin hypercubes of run sizes that are two or multiples of four can be used. As a result, Theorem 19.2.6 can be used to construct a vast number of orthogonal Latin hypercubes of $n=8k$ runs. Example 19.2.5 illustrates the use of Theorem 19.2.6.

Example 19.2.5.

Consider constructing an orthogonal Latin hypercube of 32 runs. Let $\mbox{\boldmath{$A$}}=(1,1)^{\prime}$ , $B$ be the $16\times 12$ orthogonal Latin hypercube in Table 19.6, $\mbox{\boldmath{$E$}}=(1/2,-1/2)^{\prime}$ , and $F$ be a matrix obtained by taking any 12 columns from a Hadamard matrix of order 16. By Theorem 19.2.6, $L$ in (19.17) with the chosen $\mbox{\boldmath{$A$}},\mbox{\boldmath{$B$}},\mbox{\boldmath{$E$}},\mbox{\boldmath{$F$}}$ constitutes a $32\times 12$ orthogonal Latin hypercube.

Table 19.6: A

16\times 12

orthogonal Latin hypercube from Steinberg and Lin (2006)

\mbox{\boldmath{$B$}}=\frac{1}{2}\left(\begin{array}[]{rrrrrrrrrrrrrrrr}-15&5&9&-3&7&11&-11&7&-9&3&-15&5\\ -13&1&1&13&-7&-11&11&-7&-1&-13&-13&1\\ -11&7&-7&-11&13&-1&-1&-13&9&-3&15&-5\\ -9&3&-15&5&-13&1&1&13&1&13&13&-1\\ -7&-11&11&-7&11&-7&7&11&5&15&-3&-9\\ -5&-15&3&9&-11&7&-7&-11&13&-1&-1&-13\\ -3&-9&-5&-15&1&13&13&-1&-5&-15&3&9\\ -1&-13&-13&1&-1&-13&-13&1&-13&1&1&13\\ 1&13&13&-1&-9&3&-15&5&11&-7&7&11\\ 3&9&5&15&9&-3&15&-5&3&9&5&15\\ 5&15&-3&-9&-3&-9&-5&-15&-11&7&-7&-11\\ 7&11&-11&7&3&9&5&15&-3&-9&-5&-15\\ 9&-3&15&-5&-5&-15&3&9&-7&-11&11&-7\\ 11&-7&7&11&5&15&-3&-9&-15&5&9&-3\\ 13&-1&-1&-13&-15&5&9&-3&7&11&-11&7\\ 15&-5&-9&3&15&-5&-9&3&15&-5&-9&3\\ \end{array}\right)

When $n_{1}=n_{2}$ , a stronger result than Theorem 19.2.6, as given in Theorem 19.2.7, can be established. It provides orthogonal Latin hypercubes with more columns than those in Theorem 19.2.6.

Theorem 19.2.7.

If $n_{1}=n_{2}$ and $A$ , $B$ , $E$ , and $F$ are chosen according to Theorem 19.2.6, then design ( $\mbox{\boldmath{$L$}},\mbox{\boldmath{$U$}}$ ) is an orthogonal Latin hypercube with $2k_{1}k_{2}$ factors, where $L$ is as in Theorem 19.2.6 and $\mbox{\boldmath{$U$}}=-n_{1}\mbox{\boldmath{$A$}}\otimes\mbox{\boldmath{$B$}}+\mbox{\boldmath{$E$}}\otimes\mbox{\boldmath{$F$}}$ .

Example 19.2.6.

To construct orthogonal Latin hypercubes of 64 runs, let $n_{1}=n_{2}=8$ and take

{\small\mbox{\boldmath{$A$}}=\mbox{\boldmath{$F$}}=\left(\begin{array}[]{rrrr}1&1&1&1\\[-6.0pt] 1&1&-1&-1\\[-6.0pt] 1&-1&1&-1\\[-6.0pt] 1&-1&-1&1\\[-6.0pt] 1&1&1&1\\[-6.0pt] 1&1&-1&-1\\[-6.0pt] 1&-1&1&-1\\[-6.0pt] 1&-1&-1&1\\ \end{array}\right)},\hbox{ and }\ {\small\mbox{\boldmath{$B$}}=\mbox{\boldmath{$E$}}=\frac{1}{2}\left(\begin{array}[]{rrrr}1&-3&7&5\\[-6.0pt] 3&1&5&-7\\[-6.0pt] 5&-7&-3&-1\\[-6.0pt] 7&5&-1&3\\[-6.0pt] -1&3&-7&-5\\[-6.0pt] -3&-1&-5&7\\[-6.0pt] -5&7&3&1\\[-6.0pt] -7&-5&1&-3\\ \end{array}\right)}.

Then design $(\mbox{\boldmath{$L$}},\mbox{\boldmath{$U$}})$ in Theorem 19.2.7 is a $64\times 32$ orthogonal Latin hypercube.

Theorem 19.2.8.

Suppose that an ${\mbox{OLH}}(n,k)$ is available, where $n$ is a multiple of 4 such that a Hadamard matrix of order $n$ exists. Then the following orthogonal Latin hypercubes, an OLH( $2n,k$ ), an OLH( $4n,2k$ ), an OLH( $8n,4k$ ) and an OLH( $16n,8k$ ), can all be constructed.

We give a sketch of the proof for Theorem 19.2.8. The proof provides the actual construction of these orthogonal Latin hypercubes. The theorem results from an application of Theorem 19.2.6 and the use of orthogonal designs in Table 19.7. Note that each of the four matrices in Table 19.7 can be written as $(\mbox{\boldmath{$X$}}^{\prime},-\mbox{\boldmath{$X$}}^{\prime})^{\prime}$ . From such an $X$ , define $S$ to be the matrix obtained by choosing $x_{i}=1$ for all $i$ ’s. Let $\mbox{\boldmath{$A$}}=(\mbox{\boldmath{$S$}}^{\prime},\mbox{\boldmath{$S$}}^{\prime})^{\prime}$ . Further let $E$ be an orthogonal Latin hypercube derived from a matrix in Table 19.7 by letting $x_{i}=(2i-1)/2$ for $i=1,\ldots,n/2$ . Now we choose $B$ to be a given ${\mbox{OLH}}(n,k)$ and $F$ be the matrix obtained by taking any $k$ columns from a Hadamard matrix order $n$ . Such matrices $A$ , $B$ , $E$ , and $F$ meet conditions (i), (ii), (iii) and (iv) in Theorem 19.2.6, from which Theorem 19.2.8 follows.

Table 19.7: Four orthogonal designs

$n$
2	$4$	$8$	$16$
$\begin{array}[]{r}x_{1}\\ -x_{1}\\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \\ \end{array}$	$\begin{array}[]{rr}x_{1}&x_{2}\\ x_{2}&-x_{1}\\ -x_{1}&-x_{2}\\ -x_{2}&x_{1}\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ &\\ \end{array}$	$\begin{array}[]{rrrr}x_{1}&-x_{2}&x_{4}&x_{3}\\ x_{2}&x_{1}&x_{3}&-x_{4}\\ x_{3}&-x_{4}&-x_{2}&-x_{1}\\ x_{4}&x_{3}&-x_{1}&x_{2}\\ -x_{1}&x_{2}&-x_{4}&-x_{3}\\ -x_{2}&-x_{1}&-x_{3}&x_{4}\\ -x_{3}&x_{4}&x_{2}&x_{1}\\ -x_{4}&-x_{3}&x_{1}&-x_{2}\\ &&&\\ &&&\\ &&&\\ &&&\\ &&&\\ &&&\\ &&&\\ &&&\\ \end{array}$	$\begin{array}[]{rrrrrrrr}x_{1}&-x_{2}&-x_{4}&-x_{3}&-x_{8}&x_{7}&x_{5}&x_{6}\\ x_{2}&x_{1}&-x_{3}&x_{4}&-x_{7}&-x_{8}&-x_{6}&x_{5}\\ x_{3}&-x_{4}&x_{2}&x_{1}&-x_{6}&-x_{5}&x_{7}&-x_{8}\\ x_{4}&x_{3}&x_{1}&-x_{2}&-x_{5}&x_{6}&-x_{8}&-x_{7}\\ x_{5}&-x_{6}&-x_{8}&x_{7}&x_{4}&x_{3}&-x_{1}&-x_{2}\\ x_{6}&x_{5}&-x_{7}&-x_{8}&x_{3}&-x_{4}&x_{2}&-x_{1}\\ x_{7}&-x_{8}&x_{6}&-x_{5}&x_{2}&-x_{1}&-x_{3}&x_{4}\\ x_{8}&x_{7}&x_{5}&x_{6}&x_{1}&x_{2}&x_{4}&x_{3}\\ -x_{1}&x_{2}&x_{4}&x_{3}&x_{8}&-x_{7}&-x_{5}&-x_{6}\\ -x_{2}&-x_{1}&x_{3}&-x_{4}&x_{7}&x_{8}&x_{6}&-x_{5}\\ -x_{3}&x_{4}&-x_{2}&-x_{1}&x_{6}&x_{5}&-x_{7}&x_{8}\\ -x_{4}&-x_{3}&-x_{1}&x_{2}&x_{5}&-x_{6}&x_{8}&x_{7}\\ -x_{5}&x_{6}&x_{8}&-x_{7}&-x_{4}&-x_{3}&x_{1}&x_{2}\\ -x_{6}&-x_{5}&x_{7}&x_{8}&-x_{3}&x_{4}&-x_{2}&x_{1}\\ -x_{7}&x_{8}&-x_{6}&x_{5}&-x_{2}&x_{1}&x_{3}&-x_{4}\\ -x_{8}&-x_{7}&-x_{5}&-x_{6}&-x_{1}&-x_{2}&-x_{4}&-x_{3}\\ \end{array}$

Theorem 19.2.8 is a very powerful result. By repeatedly applying Theorem 19.2.8, one can obtain many infinite series of orthogonal Latin hypercubes. For example, starting with an ${\mbox{OLH}}(12,6)$ , we can obtain an ${\mbox{OLH}}(192,48)$ , which can be used in turn to construct an ${\mbox{OLH}}(768,96)$ and so on. For another example, an ${\mbox{OLH}}(256,248)$ in Steinberg and Lin (2006) can be used to construct an ${\mbox{OLH}}(1024,496)$ , an ${\mbox{OLH}}(4096,1984)$ and so on.

Another result from Lin et al. (2010) shows how the method in (19.17) can be used to construct nearly orthogonal Latin hypercubes.

Theorem 19.2.9.

Suppose that condition (iv) in Theorem 19.2.6 is satisfied so that design $L$ in (19.17) is a Latin hypercube. If conditions (i) and (iii) in Theorem 19.2.6 hold, we then have that
(i) $\rho_{ave}^{2}(\mbox{\boldmath{$L$}})=w_{1}\rho_{ave}^{2}(\mbox{\boldmath{$B$}})+w_{2}\rho_{ave}^{2}(\mbox{\boldmath{$E$}})$ , and
(ii) $\rho_{M}(\mbox{\boldmath{$L$}})=\hbox{max}\{w_{3}\rho_{M}(\mbox{\boldmath{$B$}}),w_{4}\rho_{M}(\mbox{\boldmath{$E$}})\}$ ,
where $w_{1}$ , $w_{2}$ , $w_{3}$ and $w_{4}$ are given by $w_{1}=(k_{2}-1)(n_{2}^{2}-1)^{2}/[(k_{1}k_{2}-1)(n^{2}-1)^{2}]$ , $w_{2}=n_{2}^{4}(k_{1}-1)(n_{1}^{2}-1)^{2}/[(k_{1}k_{2}-1)(n^{2}-1)^{2}]$ , $w_{3}=(n_{2}^{2}-1)/(n^{2}-1)$ and $w_{4}=n_{2}^{2}(n_{1}^{2}-1)/(n^{2}-1)$ .

Theorem 19.2.9 says that if $B$ and $E$ are nearly orthogonal Latin hypercubes, the resulting Latin hypercube $L$ is also nearly orthogonal. An example, illustrating the use of this result, is given below.

Table 19.8: Design matrix of

\mbox{\boldmath{$B$}}_{0}

in Example 19.2.7

\left(\begin{array}[]{rrrrrrrrrrrrrrr}-15&15&-13&13&-5&-13&5&3&-1&5&-7&5&-9&-9&5\\ -13&-15&-3&3&7&3&15&-11&13&-5&7&-13&-7&-3&-3\\ -11&-9&-5&-11&-15&13&-5&11&-9&9&9&3&-5&-1&-11\\ -9&-1&9&-15&-11&1&-1&-13&5&-1&-15&7&1&3&15\\ -7&1&-7&7&15&15&-13&9&-5&-13&-3&-1&-1&7&13\\ -5&13&11&-5&9&-7&-3&-9&-13&11&13&-9&-3&13&1\\ -3&-5&13&15&-9&-9&-11&1&7&-9&15&11&9&1&-1\\ -1&-11&3&-7&11&-15&13&15&-7&-3&-9&9&7&9&-5\\ 1&3&-9&-3&-1&-5&-15&-1&11&3&-11&-15&15&5&-15\\ 3&-3&15&11&3&9&1&-7&-15&1&-13&-3&3&-15&-9\\ 5&9&7&-1&5&11&9&13&15&15&5&1&11&-7&9\\ 7&7&-1&-13&13&-1&-7&-5&9&-7&3&15&-13&-11&-13\\ 9&5&-11&-9&-7&-3&7&-3&-11&-15&11&-7&13&-13&7\\ 11&11&5&5&-13&7&11&5&3&-11&-5&-5&-11&15&-7\\ 13&-7&-15&9&1&5&3&-15&-3&13&1&13&5&11&3\\ 15&-13&1&1&-3&-11&-9&7&1&7&-1&-11&-15&-5&11\end{array}\right)

Example 19.2.7.

Let $\mbox{\boldmath{$A$}}=(1,1)^{\prime}$ and $\mbox{\boldmath{$E$}}=(1/2,-1/2)^{\prime}$ . Choose a $16\times 15$ nearly orthogonal Latin hypercube $\mbox{\boldmath{$B$}}=\mbox{\boldmath{$B$}}_{0}/2$ where $\mbox{\boldmath{$B$}}_{0}$ is displayed in Table 19.8, and $B$ has $\rho_{ave}^{2}(\mbox{\boldmath{$B$}})=0.0003$ and $\rho_{M}(\mbox{\boldmath{$B$}})=0.0765$ . Taking any 15 columns of a Hadamard matrix of order 16 to be $F$ and then applying (19.17), we obtain a Latin hypercube $L$ of 32 runs and 15 factors. As $\rho_{ave}^{2}(\mbox{\boldmath{$E$}})=\rho_{M}(\mbox{\boldmath{$E$}})=0$ , we have $\rho_{ave}^{2}(\mbox{\boldmath{$L$}})=(n_{2}^{2}-1)^{2}\rho_{ave}^{2}(\mbox{\boldmath{$B$}})/(n^{2}-1)^{2}=0.00002$ and $\rho_{M}(\mbox{\boldmath{$L$}})=(n_{2}^{2}-1)\rho_{M}(\mbox{\boldmath{$B$}})/(n^{2}-1)=0.0191$ .

Existence of orthogonal Latin hypercubes

A problem, of at least theoretical importance, in the study of orthogonal Latin hypercubes is to determine the maximum number $k^{*}$ of columns in an orthogonal Latin hypercube of a given run size $n$ . Theorem 19.2.4 tells us that $k^{*}=1$ if $n$ is 3 or $n=4m+2$ for any non-negative integer $m$ and $k^{*}\geq 2$ otherwise. Lin et al. (2010) obtained a stronger result.

Theorem 19.2.10.

The maximum number $k^{*}$ of factors for an orthogonal Latin hypercube of $n=16m+j$ runs has a lower bound given below:
(i) $k^{*}\geq 6$ for all $n=16m+j$ where $m\geq 1$ and $j\neq 2,6,10,14$ ;
(ii) $k^{*}\geq 7$ for $n=16m+11$ where $m\geq 0$ ;
(iii) $k^{*}\geq 12$ for $n=16m,16m+1$ where $m\geq 2$ ;
(iv) $k^{*}\geq 24$ for $n=32m,32m+1$ where $m\geq 2$ ;
(v) $k^{*}\geq 48$ for $n=64m,64m+1$ where $m\geq 2$ .

The above theorem provides a general lower bound on the maximum number $k^{*}$ of factors for an orthogonal Latin hypercube of $n$ runs. We now summarize the results on the best lower bound on the maximum number $k^{*}$ in an ${\mbox{OLH}}(n,k^{*})$ from all existing approaches for $n\leq 256$ . Table 19.9 lists the best lower bound on the maximum number $k^{*}$ in an ${\mbox{OLH}}(n,k^{*})$ for $n\leq 24$ . These values except the case $n=16$ were obtained by Lin (2008) through an algorithm. For $n=16$ , Steinberg and Lin (2006) obtained an orthogonal Latin hypercube with 12 columns. Table 19.10 reports the best lower bound on the maximum number $k^{*}$ in an ${\mbox{OLH}}(n,k^{*})$ for $24<n\leq 256$ as well as the corresponding approach for achieving the best lower bound.

Table 19.9: The best lower bound

k

on the maximum number

k^{*}

of factors in

{\mbox{OLH}}(n,k^{*})

for

n\leq 24

$n$	4	5	7	8	9	11	12	13	15	16	17	19	20	21	23	24
$k$	2	2	3	4	5	7	6	6	6	12	6	6	6	6	6	6

Table 19.10: The best lower bound

k

on the maximum number

k^{*}

of factors in

{\mbox{OLH}}(n,k^{*})

for

n>24

$n$	$k$	Reference	$n$	$k$	Reference
25	12	Lin et al. (2009)	144	24	Lin et al. (2010)
32	16	Sun et al. (2009)	145	12	Lin et al. (2010)
33	16	Sun et al. (2009)	160	24	Lin et al. (2010)
48	12	Lin et al. (2010)	161	24	Lin et al. (2010)
49	24	Lin et al. (2009)	169	84	Lin et al. (2009)
64	32	Sun et al. (2009)	176	12	Lin et al. (2010)
65	32	Sun et al. (2009)	177	12	Lin et al. (2010)
80	12	Lin et al. (2010)	192	48	Lin et al. (2010)
81	50	Lin et al. (2009)	193	48	Lin et al. (2010)
96	24	Lin et al. (2010)	208	12	Lin et al. (2010)
97	24	Lin et al. (2010)	209	12	Lin et al. (2010)
112	12	Lin et al. (2010)	224	24	Lin et al. (2010)
113	12	Lin et al. (2010)	225	24	Lin et al. (2010)
121	84	Lin et al. (2009)	240	12	Lin et al. (2010)
128	64	Sun et al. (2009)	241	12	Lin et al. (2010)
129	64	Sun et al. (2009)	256	248	Steinberg and Lin (2006)

19.3 Other Space-filling Designs

Section 19.2 discussed various Latin hypercube designs that are suitable for computer experiments. A Latin hypercube design does not have repeated runs and each of its factors has as many levels as the run size. Bingham et al. (2009) argued that it is absolutely unnecessary to have the same number of levels as the run size in many practical applications. They proposed the use of orthogonal and nearly orthogonal designs for computer experiments, where each factor is allowed to have repeated levels. This is a rich class of orthogonal designs, including two-level orthogonal designs and orthogonal Latin hypercubes as special cases. This section reviews the concept and constructions of orthogonal designs. We also review another class of orthogonal designs provided by Moon et al. (2011). Other classes of space-filling designs that do not fall under Latin hypercube designs are low-discrepancy sequences and uniform designs. Both types of designs originate from the field of numerical analysis and give rise to attractive space-filling designs. We provide a brief account of low-discrepancy sequences and review various measures of uniformity in uniform designs.

19.3.1 Orthogonal designs with many levels

Consider designs of $n$ runs for $k$ factors each of $s$ levels, where $2\leq s\leq n$ . For convenience, the $s$ levels are chosen to be centered and equally spaced; one such choice is $-(s-1)/2,-(s-3)/2,\ldots,(s-3)/2,(s-1)/2$ . Such a design is denoted by ${\mbox{D}}(n,s^{k})$ and can be represented by an $n\times k$ design matrix $\mbox{\boldmath{$D$}}=(d_{ij})$ with entries from the above set of $s$ levels. A Latin hypercube of $n$ runs for $k$ factors is a ${\mbox{D}}(n,s^{k})$ with $n=s$ .

Let $\mbox{\boldmath{$A$}}=(a_{ij})$ be an $n_{1}\times k_{1}$ matrix with entries $a_{ij}=\pm 1$ and $\mbox{\boldmath{$D$}}_{0}$ be a ${\mbox{D}}(n_{2},s^{k_{2}})$ . Bingham et al. (2009) constructed the $(n_{1}n_{2})\times(k_{1}k_{2})$ design

\mbox{\boldmath{$D$}}=\mbox{\boldmath{$A$}}\otimes\mbox{\boldmath{$D$}}_{0}.

(19.18)

If $A$ is column-orthogonal, then design $D$ in (19.18) is orthogonal if and only if $\mbox{\boldmath{$D$}}_{0}$ is orthogonal. This provides a powerful way to construct a rich class of orthogonal designs for computer experiments, as illustrated by Example 19.3.1.

Example 19.3.1.

Let $\mbox{\boldmath{$D$}}_{0}$ be the orthogonal Latin hypercube ${\mbox{OLH}}(16,12)$ constructed by Steinberg and Lin (2006). The construction method in (19.18) gives a series of orthogonal designs of $16m$ runs for $12m$ factors by letting $A$ be a Hadamard matrix of order $m$ , where $m$ is an integer such that a Hadamard matrix of order $m$ exists.

Higher order orthogonality and near orthogonality of $D$ in (19.18) were also discussed in Bingham et al. (2009). They considered two generalizations of the method (19.18). Let $\mbox{\boldmath{$D$}}_{j}$ be a ${\mbox{D}}(n_{2},s^{k_{2}})$ , for each $j=1,\ldots,k_{1}$ . One generalization is

\mbox{\boldmath{$D$}}=(a_{ij}\mbox{\boldmath{$D$}}_{j})=\left[\begin{array}[]{cccc}a_{11}\mbox{\boldmath{$D$}}_{1}&a_{12}\mbox{\boldmath{$D$}}_{2}&\ldots&a_{1k_{1}}\mbox{\boldmath{$D$}}_{k_{1}}\\ a_{21}\mbox{\boldmath{$D$}}_{1}&a_{22}\mbox{\boldmath{$D$}}_{2}&\ldots&a_{2k_{1}}\mbox{\boldmath{$D$}}_{k_{1}}\\ \vdots&\vdots&\ddots&\vdots\\ a_{n_{1}1}\mbox{\boldmath{$D$}}_{1}&a_{n_{1}2}\mbox{\boldmath{$D$}}_{2}&\ldots&a_{n_{1}k_{1}}\mbox{\boldmath{$D$}}_{k_{1}}\\ \end{array}\right].

(19.19)

The following results study the orthogonality of design $D$ in (19.19).

Theorem 19.3.1.

Let $A$ be column-orthogonal. We have that
(i) $\rho_{M}(\mbox{\boldmath{$D$}})=\hbox{max}\{\rho_{M}(\mbox{\boldmath{$D$}}_{1}),\ldots,\rho_{M}(\mbox{\boldmath{$D$}}_{k_{1}})\}$ ,
(ii) $\rho^{2}_{ave}(\mbox{\boldmath{$D$}})=w[\rho^{2}_{ave}(\mbox{\boldmath{$D$}}_{1})+\ldots+\rho^{2}_{ave}(\mbox{\boldmath{$D$}}_{k_{1}})]/k_{1}$ , where $w=(k_{2}-1)/(k_{1}k_{2}-1)$ , and
(iii) $D$ in (19.19) is orthogonal if and only if $\mbox{\boldmath{$D$}}_{1},\ldots,\mbox{\boldmath{$D$}}_{k_{1}}$ are all orthogonal.

The generalization (19.19) constructs designs with improved projection properties Bingham et al. (2009). The research on orthogonal designs was further pursued by Georgiou (2011) who proposed an alternative construction method and obtained many new designs.

Another class of orthogonal designs is Gram-Schmidt designs constructed by Moon et al. (2011). A Gram-Schmidt design for $n$ observations and $k$ inputs is generated from an $n\times k$ Latin hypercube design $\mbox{\boldmath{$D$}}=(d_{ij})=(\mbox{\boldmath{$d$}}_{1},\ldots,\mbox{\boldmath{$d$}}_{k})$ as follows.

Step 1:

Center the $j$ th column of $D$ to have mean zero:

$\mbox{\boldmath{$v$}}_{j}=\mbox{\boldmath{$d$}}_{j}-\sum_{i=1}^{n}d_{ij}/n,\hbox{for}\ j=1,\ldots,k.$

Step 2:

Apply the Gram-Schimidt algorithm to $\mbox{\boldmath{$v$}}_{1},\ldots,\mbox{\boldmath{$v$}}_{k}$ from Step 1 to form orthogonal columns

\mbox{\boldmath{$u$}}_{j}=\left\{\begin{array}[]{ll}\mbox{\boldmath{$v$}}_{1},&j=1;\\ \mbox{\boldmath{$v$}}_{j}-\sum_{i=1}^{j-1}\frac{\mbox{\boldmath{$u$}}_{i}\mbox{\boldmath{$v$}}_{j}}{||\mbox{\boldmath{$u$}}_{i}||^{2}}\mbox{\boldmath{$u$}}_{i},&j=2,\ldots,k,\end{array}\right.

where $||\mbox{\boldmath{$u$}}_{i}||$ represents $l_{2}$ norm of $\mbox{\boldmath{$u$}}_{i}$ .

Step 3:

Scale $\mbox{\boldmath{$u$}}_{j}$ from Step 2 to the desired range and denote the resulting column by $\mbox{\boldmath{$x$}}_{j}$ . Set $\mbox{\boldmath{$X$}}=(\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{k})$ .

Any two columns of design $X$ constructed via the three steps above have zero correlation.

19.3.2 Low-discrepancy sequences and uniform designs

Many problems in various fields such as quantum physics and computational finance require calculating definite integrals of a function over a multi-dimensional unit cube. It is very common that the function may be so complicated that the integral cannot be obtained analytically and precisely, which calls for numerical methods of approximating the integral.

Recall the numerical integration problem discussed in Section 19.2.1. The quantity $\hat{\mu}$ in (19.3) is used to approximate $\mu$ in (19.2). Respecting the common notations, we use $s$ to denote the number of factors and $\mbox{ \boldmath${\chi}$}=[0,1]^{s}$ the design region in this section. Let $\mathcal{P}=(\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n})$ be a set of $n$ points in ${\chi}$ . The bound of the integration error is given by Koksma-Hlawka inequality,

\Big{|}\mu-\hat{\mu}\Big{|}\leq V(f)D^{*}(\mathcal{P}),

(19.20)

where $V(f)$ is the variation of $f$ in the sense of Hardy and Krause and $D^{*}(\mathcal{P})$ is the star discrepancy of the $n$ points $\mathcal{P}$ Weyl (1916) and described below. Motivated by the fact that the Koksma-Hlawka bound in (19.20) is proportional to the star discrepancy of the points, different methods for generating points in ${\chi}$ with as small a star discrepancy as possible have been proposed. Such methods are referred to as quasi-Monte Carlo methods Niederreiter (1992).

For each $\mbox{\boldmath{$x$}}=(x_{1},\ldots,x_{s})$ in ${\chi}$ , let $J_{\mbox{\boldmath{$x$}}}=[0,\mbox{\boldmath{$x$}})$ denote the interval $[0,x_{1})\times\cdots\times[0,x_{s})$ , $N(\mathcal{P},J_{\mbox{\boldmath{$x$}}})$ denote the number of points of $\mathcal{P}$ falling in $J_{\mbox{\boldmath{$x$}}}$ , and $\hbox{Vol}(J_{\mbox{\boldmath{$x$}}})$ be the volume of interval $J_{\mbox{\boldmath{$x$}}}$ . The star discrepancy $D^{*}(\mathcal{P})$ of $\mathcal{P}$ is defined by

D^{*}(\mathcal{P})=\max_{\mbox{\boldmath{$x$}}\in\mbox{ \boldmath${\chi}$}}\Big{|}\frac{N(\mathcal{P},J_{\mbox{\boldmath{$x$}}})}{n}-\hbox{Vol}(J_{\mbox{\boldmath{$x$}}})\Big{|}.

(19.21)

A sequence $S$ of points in ${\chi}$ is called a low-discrepancy sequence if its first $n$ points have

D^{*}(\mathcal{P})=O(n^{-1}(\hbox{log}\ n)^{s}),

where $O(\cdot)$ is big O notation. As a comparison, if the set $\mathcal{P}$ is chosen by the Monte Carlo method, that is, $\mbox{\boldmath{$x$}}_{1},\ldots,\mbox{\boldmath{$x$}}_{n}$ are independent random samples from the uniform distribution, then $D^{*}(\mathcal{P})=O(n^{-1/2})$ , which is considered too slow in many applications Niederreiter (2012).

Construction of low-discrepancy sequences is a very active research area in the study of quasi-Monte Carlo methods. There are many constructions available; examples are the good lattice point method, the good point method, Halton sequences, Faure sequences and $(t,s)$ -sequences. For a comprehensive treatment of low-discrepancy sequences, see Niederreiter (1992). Here we provide a brief account of two popular and most widely studied methods, $(t,s)$ -sequences and uniform designs.

$(t,m,s)$ -nets and $(t,s)$ -sequences

The definitions of $(t,m,s)$ -nets and $(t,s)$ -sequences require a concept of elementary intervals. An elementary interval in base $b$ is an interval $E$ in $[0,1)^{s}$ of the form

E=\prod_{i=1}^{s}\Big{[}\frac{a_{i}}{b^{d_{i}}},\frac{a_{i}+1}{b^{d_{i}}}\Big{)}

(19.22)

with integers $a_{i}$ and $d_{i}$ satisfying $d_{i}\geq 0$ and $0\leq a_{i}<b^{d_{i}}$ . For $i=1,\ldots,s$ , the $i$ th axis of an elementary interval has length $b^{-d_{i}}$ and thus an elementary interval has a volume $b^{-\sum_{i=1}^{s}d_{i}}$ .

For integers $b\geq 2$ and $0\leq t\leq m$ , a $(t,m,s)$ -net in base $b$ is a set of $b^{m}$ points in $[0,1)^{s}$ such that every elementary interval in base $b$ of volume $b^{t-m}$ contains exactly $b^{t}$ points. For given values of $b$ , $m$ and $s$ , a smaller value of $t$ leads to a smaller elementary interval, and thus a set of points with better uniformity. Consequently, a smaller value of $t$ in $(t,m,s)$ -nets in base $b$ is preferred.

An infinite sequence of points $\{\mbox{\boldmath{$x$}}_{n}\}$ in $[0,1)^{s}$ is a $(t,s)$ -sequence in base $b$ if for all $k\geq 0$ and $m>t$ , the finite sequence $\mbox{\boldmath{$x$}}_{kb^{m}+1},\ldots,\mbox{\boldmath{$x$}}_{(k+1)b^{m}}$ forms a $(t,m,s)$ -net in base $b$ . Example 19.3.2 illustrates both concepts.

Example 19.3.2.

Consider a $(0,2)$ -sequence in base $2$ . Its first 8 points form a $(0,3,2)$ -net in base $2$ and are displayed in Figure 19.5 with $t=0,m=3,s=2$ . There are four types of elementary intervals in base $2$ of volume $2^{-3}$ with $(d_{1},d_{2})$ ’s in (19.22) being $(0,3)$ , $(3,0)$ , $(1,2)$ , and $(2,1)$ . Figures 19.5(a) - 19.5(d) show a $(0,3,2)$ -net in base $2$ when elementary intervals are given by $(d_{1},d_{2})=(0,3)$ , $(d_{1},d_{2})=(3,0)$ , $(d_{1},d_{2})=(1,2)$ , and $(d_{1},d_{2})=(2,1)$ , respectively. Note that in every elementary interval of the form

\Big{[}\frac{a_{1}}{2^{d_{1}}},\frac{(a_{1}+1)}{2^{d_{1}}}\Big{)}\times\Big{[}\frac{a_{2}}{2^{d_{2}}},\frac{(a_{2}+1)}{2^{d_{2}}}\Big{)},\ 0\leq a_{i}<2^{d_{i}},\ i=1,2,

there is exactly one point. The next 8 points in this (0,2)-sequence in base 2 also form a $(0,3,2)$ -net in base 2. The totality of all 16 points is a $(0,4,2)$ -net in base 2. Analogous to Figure 19.5, Figure 19.6 exhibits the $(0,4,2)$ -net in base 2 when elementary intervals are given by all $(d_{1},d_{2})$ ’s that satisfy $d_{1}+d_{2}=m=4$ .

A general theory of $(t,m,s)$ -nets and $(t,s)$ -sequences was developed by Niederreiter (1987). Some special cases of $(t,s)$ -sequences are as follows. Sobol’ sequences Sobol’ (1967) are $(t,s)$ -sequences in base 2. Faure sequences Faure (1982) are $(0,s)$ -sequences in base $q$ where $q$ is a prime with $s\leq q$ . Niederreiter sequences Niederreiter (1987) are $(0,s)$ sequences in base $q$ where $q$ is a prime or a prime power with $s\leq q$ . Niederreiter-Xing sequences Niederreiter and Xing (1996) are $(t,s)$ -sequences in base $q$ for some certain $t$ where $q$ is a prime or a prime power with $s>q$ . For constructions of all these sequences, we refer the readers to Niederreiter (2008). Results on existing $(t,s)$ -sequences are available in Schürer and Schmid (2010).

Uniform designs

Motivated by the Koksma-Hlawka inequality in (19.20), Fang and Wang Fang (1980); Wang and Fang (1981) introduced uniform designs, and by their definition, a uniform design is a set of design points with the smallest discrepancy among all possible designs of the same run size. One choice of discrepancy is the star discrepancy in (19.21). More generally, one can use the $L_{p}$ discrepancy,

D_{p}(\mathcal{P})=\Big{[}\int_{\mbox{ \boldmath${\chi}$}}\big{|}\frac{N(\mathcal{P},J_{\mbox{\boldmath{$x$}}})}{n}-\hbox{Vol}(J_{\mbox{\boldmath{$x$}}})\big{|}^{p}d\mbox{\boldmath{$x$}}\Big{]}^{1/p},

where $N(\mathcal{P},J_{\mbox{\boldmath{$x$}}})$ and $\hbox{Vol}(J_{\mbox{\boldmath{$x$}}})$ are defined as in (19.21). Two special cases of the $L_{p}$ discrepancy are the $L_{\infty}$ discrepancy, which is the star discrepancy, and the $L_{2}$ discrepancy. While the $L_{\infty}$ discrepancy is difficult to compute, the $L_{2}$ discrepancy is much easier to evaluate because of a simple formula given by Warnock (1972),

D_{2}(\mathcal{P})=2^{-s}-\frac{2^{1-s}}{n}\sum_{i=1}^{n}\prod_{l=1}^{s}(1-x_{il}^{2})+\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\prod_{l=1}^{s}[1-\hbox{max}(x_{il},x_{jl})],

where $x_{il}$ is the setting of the $l$ th factor in the $i$ th run, $i=1,\ldots,n$ and $l=1,\ldots,s$ .

The $L_{p}$ discrepancy aims to achieve uniformity in the $s$ -dimensional design space. Designs with the smallest $L_{p}$ discrepancy do not necessarily perform well in terms of projection uniformity in low dimensions. Hickernell (1998) proposed three new measures of uniformity, the symmetric $L_{2}$ discrepancy ( $SL_{2}$ ), the centered $L_{2}$ discrepancy ( $CL_{2}$ ), and the modified $L_{2}$ discrepancy ( $ML_{2}$ ). They are all defined through

D_{mod}(\mathcal{P})=\sum_{u\neq\emptyset}\int_{\mbox{ \boldmath${\chi}$}_{u}}\Big{|}\frac{N(\mathcal{P}_{u},J_{\mbox{\boldmath{$x$}}_{u}})}{n}-\hbox{Vol}(J_{\mbox{\boldmath{$x$}}_{u}})\Big{|}^{2}du,

(19.23)

where $\emptyset$ represents the empty set, $u$ is a non-empty subset of the set $\{1,\ldots,s\}$ , $|u|$ denotes the cardinality of $u$ , $\mbox{ \boldmath${\chi}$}_{u}$ is the $|u|$ -dimensional unit cube involving the coordinates in $u$ , $\mathcal{P}_{u}$ is the projection of the set of points $\mathcal{P}$ on $\mbox{ \boldmath${\chi}$}_{u}$ , $J_{\mbox{\boldmath{$x$}}_{u}}$ is the projection of $J_{\mbox{\boldmath{$x$}}}$ on $\mbox{ \boldmath${\chi}$}_{u}$ , $N(\mathcal{P}_{u},J_{\mbox{\boldmath{$x$}}_{u}})$ denotes the number of points of $\mathcal{P}_{u}$ falling in $J_{\mbox{\boldmath{$x$}}_{u}}$ , and $\hbox{Vol}(J_{\mbox{\boldmath{$x$}}_{u}})$ represents the volume of $J_{\mbox{\boldmath{$x$}}_{u}}$ . The symmetric $L_{2}$ discrepancy chooses $J_{\mbox{\boldmath{$x$}}}$ such that it is invariant if $x_{il}$ is replaced by $1-x_{il}$ , $i=1,\ldots,n$ and $l=1,\ldots,s$ , and it has the formula

(SL_{2}(\mathcal{P}))^{2}=(\frac{4}{3})^{s}-\frac{2}{n}\sum_{i=1}^{n}\prod_{l=1}^{s}(1+2x_{il}-2x_{il}^{2})+\frac{2^{s}}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\prod_{l=1}^{s}\big{(}1-|x_{il}-x_{jl}|\big{)}.

The centered $L_{2}$ discrepancy chooses $J_{\mbox{\boldmath{$x$}}}$ such that it is invariant under the reflections of $\mathcal{P}$ around any hyperplane with the $l$ th coordinate being 0.5. Let $\mathcal{A}^{s}$ denote the set of $2^{s}$ vertices of the unit cube ${\chi}$ and $\mbox{\boldmath{$a$}}=(a_{1},\ldots,a_{s})\in\mathcal{A}^{s}$ be the closest one to $x$ . The centered $L_{2}$ discrepancy takes $J_{\mbox{\boldmath{$x$}}}$ in (19.23) to be

\{\mbox{\boldmath{$d$}}=(d_{1},\ldots,d_{s})\in\mbox{ \boldmath${\chi}$}\ |\ \hbox{min}(a_{j},x_{j})\leq d_{j}<\hbox{max}(a_{j},x_{j}),j=1,\ldots,s\}.

The formula for the centered $L_{2}$ discrepancy is given by

	$\displaystyle(CL_{2}(\mathcal{P}))^{2}$	$\displaystyle=$	$\displaystyle(\frac{13}{12})^{2}-\frac{2}{n}\sum_{i=1}^{n}\prod_{l=1}^{s}\big{(}1+\frac{1}{2}\|x_{il}-0.5\|-\frac{1}{2}\|x_{il}-0.5\|^{2}\big{)}$
			$\displaystyle+\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\prod_{l=1}^{s}\big{(}1+\frac{1}{2}\|x_{il}-0.5\|+\frac{1}{2}\|x_{jl}-0.5\|-\frac{1}{2}\|x_{il}-x_{jl}\|\big{)}.$

The modified $L_{2}$ discrepancy takes $J_{\mbox{\boldmath{$x$}}}=[\textbf{0},\mbox{\boldmath{$x$}})$ and has the formula

(ML_{2}(\mathcal{P}))^{2}=(\frac{4}{3})^{s}-\frac{2^{1-s}}{n}\sum_{i=1}^{n}\prod_{l=1}^{s}(3-x_{il}^{2})+\frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{s}\prod_{l=1}^{s}\big{[}2-\hbox{max}(x_{il},x_{il})\big{]}.

For other discrepancy measures such as the wrap-around discrepancy, see Fang et al. (2006).

Finding uniform designs based on a discrepancy criterion is an optimization problem. However, searching uniform designs in the entire unit cube is computationally prohibitive for large designs. Instead, it is convenient to find uniform designs within a class of U-type designs. Suppose that each of the $s$ factors in an experiment has $q$ levels, $\{1,\ldots,q\}$ . A U-type design, denoted by $U(n;q^{s})$ , is an $n\times s$ matrix in which the $q$ levels in each column appear equally often. Table 19.11 displays a $U(6;3^{2})$ and a $U(6;6^{2})$ . For $q=n$ , uniform $U$ -type designs can be constructed by several methods such as the good lattice method, the Latin square method, the expanding orthogonal array method and the cutting method Fang et al. (2006). For general values of $q$ , optimization algorithms have been considered, such as simulated annealing, genetic algorithm, and threshold accepting Bohachevsky et al. (1986); Holland (1975); Winker and Fang (1997). For more detailed discussions on the theory and applications of uniform designs, see Fang et al. (2000) and Fang and Lin (2003).

Table 19.11: Uniform designs

U(6;3^{2})

and

U(6;6^{2})

$U(6;3^{2})$		$U(6;6^{2})$
$\begin{array}[]{rr}1&1\\ 2&2\\ 3&3\\ 1&3\\ 2&1\\ 3&2\\ \end{array}$		$\begin{array}[]{rr}1&3\\ 2&5\\ 3&1\\ 4&6\\ 5&2\\ 6&4\\ \end{array}$

19.4 Concluding Remarks

We have provided an expository account of the constructions and properties of space-filling designs for computer experiments. Research in this area remains active and will continue to thrive. Recently, a number of new directions have been pursued. He and Tang (2013) introduced strong orthogonal arrays and associated Latin hypercubes, while Tang et al. (2012) studied uniform fractional factorial designs. Research has also been conducted to take advantage of many available results from other design areas such as factorial design theory, one such work being multi-layer designs proposed by Ba and Joseph (2011). Another important direction is to develop methodology for the design regions in which input variables have dependency or constraints; see Draguljic et al. (2012) and Bowman and Woods (2013) for more details.

References

Audze and Eglais (1977) Audze, P. and Eglais, V. (1977), “New approach to planning out of experiments,” Problems of Dynamics and Strength, 35, 104–107.
Ba and Joseph (2011) Ba, S. and Joseph, V. R. (2011), “Multi-layer designs for computer experiments,” Journal of the American Statistical Association, 106, 1139–1149.
Bingham et al. (2009) Bingham, D., Sitter, R. R., and Tang, B. (2009), “Orthogonal and nearly orthogonal designs for computer experiments,” Biometrika, 96, 51–65.
Bohachevsky et al. (1986) Bohachevsky, I. O., Johnson, M. E., and Stein, M. L. (1986), “Generalized simulated annealing for function optimization,” Technometrics, 28, 209–217.
Bowman and Woods (2013) Bowman, V. E. and Woods, D. C. (2013), “Weighted space-filling designs,” Journal of Simulation, 7, 249–263.
Carnell (2009) Carnell, R. (2009), “lhs: Latin hypercube samples,” R package version 0.5.
Chen et al. (2013) Chen, R., Hsieh, D., Hung, Y., and Wang, W. (2013), “Optimizing Latin hypercube designs by particle swarm,” Statistics and Computing, 23, 664–676.
Cioppa and Lucas (2007) Cioppa, T. M. and Lucas, T. W. (2007), “Efficient nearly orthogonal and space-filling Latin hypercubes,” Technometrics, 49, 45–55.
Conway et al. (1999) Conway, J. H., Sloane, N. J. A., and Bannai, E. (1999), Sphere Packings, Lattices, and Groups, vol. 290, Springer Verlag.
Dey and Mukerjee (1999) Dey, A. and Mukerjee, R. (1999), Fractional Factorial Plans, New York: Wiley.
Draguljic et al. (2012) Draguljic, D., Dean, A. M., and Santner, T. J. (2012), “Noncollapsing space-filling designs for bounded nonrectangular regions,” Technometrics, 54, 169–178.
Erkut (1990) Erkut, E. (1990), “The discrete p-dispersion problem,” European Journal of Operational Research, 46, 48–60.
Fang (1980) Fang, K. T. (1980), “The uniform design: application of number-theoretic methods in experimental design,” Acta Mathematicae Applicatae Sinica, 3, 363–372.
Fang et al. (2006) Fang, K. T., Li, R., and Sudjianto, A. (2006), Design and Modeling for Computer Experiments, vol. 6, CRC Press.
Fang and Lin (2003) Fang, K. T. and Lin, D. K. J. (2003), “Uniform experimental designs and their applications in industry,” Handbook of Statistics, 22, 131–170.
Fang et al. (2000) Fang, K. T., Lin, D. K. J., Winker, P., and Zhang, Y. (2000), “Uniform design: theory and application,” Technometrics, 42, 237–248.
Faure (1982) Faure, H. (1982), “Discrépance de suites associéesa un systeme de numération (en dimension s),” Acta Arith, 41, 337–351.
Forrester et al. (2008) Forrester, A., Sobester, A., and Keane, A. (2008), Engineering Design via Surrogate Modelling, Chichester: Wiley.
Georgiou (2009) Georgiou, S. D. (2009), “Orthogonal Latin hypercube designs from generalized orthogonal designs,” Journal of Statistical Planning and Inference, 139, 1530–1540.
Georgiou (2011) — (2011), “Orthogonal designs for computer experiments,” Journal of Statistical Planning and Inference, 141, 1519–1525.
Grosso et al. (2009) Grosso, A., Jamali, A., and Locatelli, M. (2009), “Finding maximin Latin hypercube designs by iterated local search heuristics,” European Journal of Operational Research, 197, 541–547.
He and Tang (2013) He, Y. and Tang, B. (2013), “Strong orthogonal arrays and associated Latin hypercubes for computer experiments,” Biometrika, 100, 254–260.
Hedayat et al. (1999) Hedayat, A., Sloane, N. J. A., and Stufken, J. (1999), Orthogonal Arrays: Theory and Applications, Springer Verlag.
Hickernell (1998) Hickernell, F. J. (1998), “A generalized discrepancy and quadrature error bound,” Mathematics of Computation, 67, 299–322.
Holland (1975) Holland, J. H. (1975), Adaptation in Natural and Artificial Systems, Ann Arbor, MI: University of Michigan Press.
Iman and Conover (1982) Iman, R. L. and Conover, W. J. (1982), “Distribution-free approach to inducing rank correlation among input variables,” Communications in Statistics - Simulation and Computation, 11, 311–334.
Jin et al. (2005) 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.
Johnson et al. (1990) Johnson, M. E., Moore, L. M., and Ylvisaker, D. (1990), “Minimax and maximin distance designs,” Journal of Statistical Planning and Inference, 26, 131–148.
Joseph and Hung (2008) Joseph, V. R. and Hung, Y. (2008), “Orthogonal-maximin Latin hypercube designs,” Statistica Sinica, 18, 171–186.
Joseph et al. (2008) Joseph, V. R., Hung, Y., and Sudjianto, A. (2008), “Blind kriging: a new method for developing metamodels,” Journal of Mechanical Design, 130, 31102.
Kuhfeld (2009) Kuhfeld, W. F. (2009), “Orthogonal arrays,” Website courtesy of SAS Institute.
Leary et al. (2003) Leary, S., Bhaskar, A., and Keane, A. (2003), “Optimal orthogonal-array-based Latin hypercubes,” Journal of Applied Statistics, 30, 585–598.
Leatherman et al. (2014) Leatherman, E., Dean, A., and Santner, T. (2014), “Designs for computer experiments that minimize the weighted integrated mean square prediction error,” Tech. Rep. 875, The Ohio State University, Columbus, Ohio.
Li and Wu (1997) Li, W. W. and Wu, C. F. J. (1997), “Columnwise-pairwise algorithms with applications to the construction of supersaturated designs,” Technometrics, 39, 171–179.
Liefvendahl and Stocki (2006) Liefvendahl, M. and Stocki, R. (2006), “A study on algorithms for optimization of Latin hypercubes,” Journal of Statistical Planning and Inference, 136, 3231–3247.
Lin (2008) Lin, C. D. (2008), New Development in Designs for Computer Experiments and Physical Experiments, Ph.D. thesis, Simon Fraser University.
Lin et al. (2010) Lin, C. D., Bingham, D., Sitter, R. R., and Tang, B. (2010), “A new and flexible method for constructing designs for computer experiments,” Annals of Statistics, 38, 1460–1477.
Lin et al. (2009) Lin, C. D., Mukerjee, R., and Tang, B. (2009), “Construction of orthogonal and nearly orthogonal Latin hypercubes,” Biometrika, 96, 243–247.
McKay et al. (1979) 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.
Melissen (1997) Melissen, J. B. M. (1997), Packing and Covering with Circles, Ph.D. thesis, Utrecht University, The Netherlands.
Moon et al. (2011) Moon, H., Dean, A. M., and Santner, T. J. (2011), “Algorithms for generating maximin orthogonal and Latin hypercube designs,” Journal of Statistical Theory and Practice, 5, 81–98.
Morris and Mitchell (1995) Morris, M. D. and Mitchell, T. J. (1995), “Exploratory designs for computer experiments,” Journal of Statistical Planning and Inference, 43, 381–402.
Mukerjee and Wu (2006) Mukerjee, R. and Wu, C. F. (2006), A Modern Theory of Factorial Designs, Springer Verlag.
Niederreiter (1987) Niederreiter, H. (1987), “Point sets and sequences with small discrepancy,” Monatshefte für Mathematik, 104, 273–337.
Niederreiter (1992) — (1992), Random Number Generation and Quasi-Monte Carlo Methods, SIAM CBMS-NSF Regional Conference Series in Applied Mathematics. Philadelphia.
Niederreiter (2008) — (2008), “Nets,(t, s)-sequences, and codes,” Monte Carlo and Quasi-Monte Carlo Methods 2006, 1, 83–100.
Niederreiter (2012) — (2012), “Low-discrepancy simulation,” Handbook of Computational Finance, 703–729.
Niederreiter and Xing (1996) Niederreiter, H. and Xing, C. (1996), “Low-discrepancy sequences and global function fields with many rational places,” Finite Fields and Their Applications, 2, 241–273.
Owen (1992) Owen, A. B. (1992), “A central limit theorem for Latin hypercube sampling,” Journal of the Royal Statistical Society. Series B (Methodological), 54, 541–551.
Owen (1994) — (1994), “Controlling correlations in Latin hypercube samples,” Journal of the American Statistical Association, 89, 1517–1522.
Pang et al. (2009) Pang, F., Liu, M. Q., and Lin, D. K. J. (2009), “A construction method for orthogonal Latin hypercube designs with prime power levels,” Statistica Sinica, 19, 1721–1728.
Patterson (1954) Patterson, H. D. (1954), “The errors of lattice sampling,” Journal of the Royal Statistical Society. Series B (Methodological), 16, 140–149.
Sacks et al. (1989) Sacks, J., Welch, W. J., Mitchell, T. J., and Wynn, H. P. (1989), “Design and analysis of computer experiments,” Statistical Science, 4, 409–423.
Santner et al. (2003) Santner, T. J., Williams, B. J., and Notz, W. (2003), The Design and Analysis of Computer Experiments, Springer Verlag.
Schürer and Schmid (2010) Schürer, R. and Schmid, W. C. (2010), “MinT – Architecture and applications of the ( $t,m,s$ )-net and OOA database,” Mathematics and Computers in Simulation, 80, 1124–1132.
Shewry and Wynn (1987) Shewry, M. C. and Wynn, H. P. (1987), “Maximum entropy sampling,” Journal of Applied Statistics, 14, 165–170.
Sobol’ (1967) Sobol’, I. M. (1967), “On the distribution of points in a cube and the approximate evaluation of integrals,” USSR Computational Mathematics and Mathematical Physics, 7, 86–112.
Stein (1987) Stein, M. (1987), “Large sample properties of simulations using Latin hypercube sampling,” Technometrics, 29, 143–151.
Steinberg and Lin (2006) Steinberg, D. M. and Lin, D. K. J. (2006), “A construction method for orthogonal Latin hypercube designs,” Biometrika, 93, 279–288.
Sun et al. (2009) Sun, F., Liu, M. Q., and Lin, D. K. J. (2009), “Construction of orthogonal Latin hypercube designs,” Biometrika, 96, 971–974.
Sun et al. (2010) — (2010), “Construction of orthogonal Latin hypercube designs with flexible run sizes,” Journal of Statistical Planning and Inference, 140, 3236–3242.
Tang (1993) Tang, B. (1993), “Orthogonal array-based Latin hypercubes,” Journal of the American Statistical Association, 88, 1392–1397.
Tang (1994) — (1994), “A theorem for selecting OA-based Latin hypercubes using a distance criterion,” Communications in Statistics - Theory and Methods, 23, 2047–2058.
Tang (1998) — (1998), “Selecting Latin hypercubes using correlation criteria,” Statistica Sinica, 8, 965–978.
Tang et al. (2012) Tang, Y., Xu, H., and Lin, D. K. J. (2012), “Uniform fractional factorial designs,” Annals of Statistics, 40, 891–907.
van Dam et al. (2007) van Dam, E. R., Husslage, B., den Hertog, D., and Melissen, H. (2007), “Maximin Latin hypercube designs in two dimensions,” Operations Research, 55, 158–169.
Viana et al. (2010) Viana, F. A. C., Venter, G., and Balabanov, V. (2010), “An algorithm for fast optimal Latin hypercube design of experiments,” International Journal for Numerical Methods in Engineering, 82, 135–156.
Wang and Fang (1981) Wang, Y. and Fang, K. T. (1981), “A note on uniform distribution and experimental design,” KeXue TongBao, 26, 485–489.
Warnock (1972) Warnock, T. T. (1972), “Computational investigations of low-discrepancy point sets,” Applications of Number Theory to Numerical Analysis, 319–343.
Weyl (1916) Weyl, H. (1916), “Über die Gleichverteilung von Zahlen mod. eins,” Mathematische Annalen, 77, 313–352.
Winker and Fang (1997) Winker, P. and Fang, K. T. (1997), “Application of threshold accepting to the evaluation of the discrepancy of a set of points,” SIAM Journal on Numerical Analysis, 34, 2028–2042.
Wu and Hamada (2011) Wu, C. F. J. and Hamada, M. S. (2011), Experiments: Planning, Analysis, and Optimization, Wiley Series in Probability and Statistics, John Wiley & Sons.
Yang and Liu (2012) Yang, J. and Liu, M. Q. (2012), “Construction of orthogonal and nearly orthogonal Latin hypercube designs from orthogonal designs,” Statistica Sinica, 22, 433–442.
Ye (1998) Ye, K. Q. (1998), “Orthogonal column Latin hypercubes and their application in computer experiments,” Journal of the American Statistical Association, 93, 1430–1439.
Ye et al. (2000) Ye, K. Q., Li, W., and Sudjianto, A. (2000), “Algorithmic construction of optimal symmetric Latin hypercube designs,” Journal of Statistical Planning and Inference, 90, 145–159.
Zhu et al. (2011) Zhu, H., Liu, L., Long, T., and Peng, L. (2011), “A novel algorithm of maximin Latin hypercube design using successive local enumeration,” Engineering Optimization, 1, 1–14.

Latin Hypercubes and Space-filling Designs

19.1 Introduction

19.2 Latin Hypercube Designs

19.2.1 Introduction and examples

Theorem 19.2.1.

Theorem 19.2.2.

19.2.2 Latin hypercube designs based on measures of distance

Maximin Latin hypercube designs

19.2.3 Orthogonal array-based Latin hypercube designs

Example 19.2.1.

Theorem 19.2.3.

19.2.4 Orthogonal and nearly orthogonal Latin hypercube designs

Theorem 19.2.4.

A construction method based on an orthogonal array and a small orthogonal Latin hypercube

Theorem 19.2.5.

Corollary 19.2.1.

Example 19.2.2.

Example 19.2.3.

A recursive construction method

Example 19.2.4.

A construction method based on small orthogonal designs and small orthogonal Latin hypercubes

Theorem 19.2.6.

Example 19.2.5.

Theorem 19.2.7.

Example 19.2.6.

Theorem 19.2.8.

Theorem 19.2.9.

Example 19.2.7.

Existence of orthogonal Latin hypercubes

Theorem 19.2.10.

19.3 Other Space-filling Designs

19.3.1 Orthogonal designs with many levels

Example 19.3.1.

Theorem 19.3.1.

19.3.2 Low-discrepancy sequences and uniform designs

(t,m,s)(t,m,s)-nets and (t,s)(t,s)-sequences

Example 19.3.2.

Uniform designs

19.4 Concluding Remarks

References

$(t,m,s)$ -nets and $(t,s)$ -sequences