A new and flexible method for constructing designs for computer experiments

Let $A=(a_{ij})_{n_{1}\times m_{1}}$ be a matrix with entries $a_{ij}=\pm 1$, $B=(b_{ij})_{n_{2}\times m_{2}}$ be a $D(n_{2},s_{2}^{m_{2}})$, $C=(c_{ij})_{n_{1}\times m_{1}}$ be a $D(n_{1},s_{1}^{m_{1}})$ and $D=(d_{ij})_{n_{2}\times m_{2}}$ be a matrix with entries $d_{ij}=\pm 1$. Let $\gamma$ be a real number. New designs are found using the following construction:

where the Kronecker product $A\otimes B$ is the $n_{1}n_{2}\times m_{1}m_{2}$ matrix,

with $a_{ij}B$ itself being an $n_{2}\times m_{2}$ matrix. The resulting design $L$ in (1) has $n=n_{1}n_{2}$ runs and $m=m_{1}m_{2}$ factors.

The above construction has an interesting interpretation. As an illustration, consider a simple case in which $A=(1,1)^{T}$ and $C=(1/2,-1/2)^{T}$. Design $L$ in (1) has a column,

where $b$ is a column of $B$ and $d$ is a column of $D$. Further let $d=(d_{1},\ldots,d_{n_{2}})^{T}$. Since $d_{i}=\pm 1$, the column (2) can be viewed as simultaneously shifting each level in $b$ to the left and the right by $\gamma/2$. If we view $b$ as a block of level settings, then we are shifting two identical blocks $b$, one to the left and the other to the right. We will show in Section 2.2 that with the appropriate choices of $A$, $B$, $C$, $D$ and $\gamma$, the levels in each column of $L$ in (1) are equally spaced and unreplicated, thus resulting in a Latin hypercube.

Now consider all $m$ columns of $L$ under this simple case. Each one-dimensional block $b$ becomes an $m$-dimensional stratum, $B$. Suppose $D$ is a matrix of plus ones. Then the design points in $B+\gamma D/2$ can be obtained by shifting the entire stratum $B$ to the right by $\gamma/2$. Similarly, the design points in $B-\gamma D/2$ can be obtained by shifting the entire stratum $B$ to the left by $\gamma/2$. In this case, closely clustered points in each stratum are expected. This feature can be utilized to construct cascading Latin hypercubes [Lin (2008)].

We shall see that the orthogonality or near orthogonality of $L$ in (1) is determined by the orthogonality or near orthogonality of $A$, $B$, $C$ and $D$, the correlations between the columns in $A$ and those in $C$, and the correlations between the columns in $B$ and those in $D$. As a result, the method allows orthogonal and nearly orthogonal Latin hypercubes to be easily constructed.

Vartak (1955) appears to be the first to use the Kronecker product systematically to construct statistical experimental designs. In a recent work, Bingham, Sitter and Tang (2009) introduced a method for constructing a rich class of designs that are suitable for use in computer experiments. Their approach corresponds to $\gamma=0$ in the general construction given in (1). The designs in that paper have many levels and are not Latin hypercubes in general.

The following result shows how to obtain Latin hypercubes from the construction in (1).

The proof is given in the Appendix. Just in terms of constructing Latin hypercubes, Lemma 1 is not of much significance in itself as one can easily obtain a Latin hypercube simply by combining several permutations of the set of levels. The significance of Lemma 1 lies in the fact that it produces Latin hypercubes with the structure in (1) and thus provides a path to the construction of orthogonal and cascading Latin hypercubes.

Condition (i) in Lemma 1 is not really a condition, and it simply tells us to choose $B$ and $C$ to be Latin hypercubes. In order for $L$ to be a Latin hypercube, the only mild condition is that in (ii) of Lemma 1. Two situations where condition (ii) is obviously met are as follows: ($\alpha$) $C$ has a foldover structure in the sense that $C=(C_{0}^{T},-C_{0}^{T})^{T}$, and $A$ has the form $A=(A_{0}^{T},A_{0}^{T})^{T}$; ($\beta$) $A$ or $D$ is a matrix of all plus ones. Both situations are useful. Theorem 3 of Section 3.3 is derived under situation ($\alpha$). Situation ($\beta$) can be used for constructing cascading Latin hypercubes. We now give an example to illustrate Lemma 1.

A design or matrix $X=(x_{1},\ldots,x_{m})$ is said to be orthogonal if the inner product of any two columns is zero, that is, $x_{i}^{T}x_{j}=0$ for all $i\neq j$. The next result provides a set of sufficient conditions for design $L$ in (1) to be orthogonal.

The proof is simple, making use of the following properties of the Kronecker product:

Lemma 1 tells how to make $L$ in (1) a Latin hypercube whereas Lemma 2 tells how to make it orthogonal. When the two lemmas are combined, we have a way of obtaining orthogonal Latin hypercubes.

The role played by $A$ and $D$ is very different from that of $B$ and $C$ in Theorem 1. To help understand Theorem 1, one may think that $B$ and $C$ are the building material while $A$ and $D$ provide a blueprint for the construction. Small orthogonal Latin hypercubes $B$ and $C$ are used to construct a large orthogonal Latin hypercube $L$ in Theorem 1. Exactly how the construction is accomplished is guided by $A$ and $D$ which are orthogonal matrices of $\pm 1$. In addition to the right blueprint and building material, a considerable amount of care is necessary for the final structure to be right. This is achieved via $\gamma=n_{2}$ and conditions (iii) and (iv) in Theorem 1.

Note that $A$ and $D$ may or may not be square matrices, and the orthogonality of $A$ and $D$ is imposed on their columns. In some mathematics literature, such matrices are called Hadamard submatrices. For convenience, we simply call $A$ or $D$ an orthogonal matrix when its columns are orthogonal. Hadamard matrices and orthogonal arrays with levels $\pm 1$ are all such orthogonal matrices in our terminology. A Hadamard matrix is a square orthogonal matrix of $\pm 1$. An orthogonal array with two levels $\pm 1$ requires that each of the four combinations $(-1,-1)$, $(-1,+1)$, $(+1,-1)$ and $(+1,+1)$ occurs the same number of times in every two columns. For some comprehensive discussion on these and other topics in the theory of factorial designs, we refer to Dey and Mukerjee (1999), Hedayat, Sloane and Stufken (1999) and Mukerjee and Wu (2006).

Because of the orthogonality of $A$ and $D$, we must have that $n_{1}$ and $n_{2}$ are equal to two or multiples of four. The case where $n_{1}=n_{2}=2$ is trivial. Consequently, Theorem 1 can be used to construct orthogonal Latin hypercubes of $n=8k$ runs, thereby providing designs that are unavailable in Ye (1998) and Steinberg and Lin (2006). When $n=n_{1}n_{2}$ is a multiple of 16, Theorem 1 becomes more powerful. This point will be highlighted in Section 3.3. We now revisit Example 1 for an illustration of Theorem 1.

We now discuss how to choose $A,B,C,D$ and $\gamma$ to construct orthogonal Latin hypercubes. According Theorem 1, we have that $\gamma=n_{2}$. Matrices $A$ and $D$ need to be orthogonal with entries of $\pm 1$. As discussed earlier, two level orthogonal arrays and Hadamard matrices are all such orthogonal matrices. Theorem 1 requires that designs $B$ and $C$ be orthogonal Latin hypercubes. All known orthogonal Latin hypercubes from the existing literature can be used here. Later in this paper (see Table 3), we obtain a collection of small orthogonal Latin hypercubes through a computer search for this purpose. So far, all are straightforward. The nontrivial aspect from applying Theorem 1 is to satisfy conditions (iii) and (iv) which require that $A$ and $C$ (or $B$ and $D$) jointly have certain properties. In this paper, we satisfy these two conditions by choosing $A$ of form $A=(A_{0}^{T},A_{0}^{T})^{T}$ and $C$ of form $C=(C_{0}^{T},-C_{0}^{T})^{T}$ where $A_{0}$ and $C_{0}$ are such that all the columns in the matrix,

are mutually orthogonal. In Section 3.3 we provide a method of finding such orthogonal matrices with the structure in (4) when proving Theorem 3. Comments similar to those in this paragraph can also be made regarding the application of Proposition 2 in Section 3.2.

Consider an orthogonal Latin hypercube of $n$ runs with $m\geq 2$ factors. Trivially, run size $n$ cannot be two or three. So we must have $n\geq 4$. The next result provides a complete characterization of the existence of an orthogonal Latin hypercube in terms of run size $n$.

The Appendix contains a proof for Theorem 2. Equivalently, Theorem 2 says that the run size of an orthogonal Latin hypercube has to be odd or a multiple of 4. Theorem 1 provides a method for constructing orthogonal Latin hypercubes of $n=8k$ runs. The present section examines how to construct orthogonal Latin hypercubes of other run sizes.

The basic idea of our method is quite simple. To obtain an orthogonal Latin hypercube, we stack up two orthogonal designs with mutually exclusive and exhaustive sets of levels. To make it precise, we use $\mathcal{S}$ to denote the set of $n$ levels of a Latin hypercube of $n$ runs. Let $\mathcal{S}=\mathcal{S}_{a}\cup\mathcal{S}_{b}$ where $\mathcal{S}_{a}\cap\mathcal{S}_{b}=\phi$, and let $n_{a}$ and $n_{b}$ be the numbers of levels in $\mathcal{S}_{a}$ and $\mathcal{S}_{b}$, respectively. Suppose that there exist an $n_{a}\times m$ orthogonal design $D_{a}$ with levels in $\mathcal{S}_{a}$ and an $n_{b}\times m$ orthogonal design $D_{b}$ with levels in $\mathcal{S}_{b}$, where for both $D_{a}$ and $D_{b}$, each level appears precisely once within each column. Then

is an $n\times m$ orthogonal Latin hypercube with $n=n_{a}+n_{b}$. Note that $D_{a}$ and $D_{b}$ themselves are not necessarily Latin hypercubes.

We consider two special choices for $\mathcal{S}_{a}$ and $\mathcal{S}_{b}$. For easy reference later in the paper, we call them two stacking methods. Our first stacking method chooses $n_{a}$ and $n_{b}$ such that $|n_{a}-n_{b}|=1$ with the corresponding $\mathcal{S}_{a}=\{-(n_{a}-1),-(n_{a}-3),\ldots,n_{a}-3,n_{a}-1\}$ and $\mathcal{S}_{b}=\{-(n_{b}-1),-(n_{b}-3),\ldots,\break n_{b}-3,n_{b}-1\}$. This implies that both $D_{a}/2$ and $D_{b}/2$ in (5) are orthogonal Latin hypercubes. We may assume that $n_{a}$ is odd and $n_{b}$ is even in the above. By Theorem 2, we know that $n_{b}$ has form $4k$. It follows that $n_{a}$ has form $4k-1$ or $4k+1$. Thus the first stacking method allows orthogonal Latin hypercubes of run sizes $8k-1$ and $8k+1$ to be constructed.

The second stacking method is more generally applicable and it chooses $\mathcal{S}_{a}=\{-(n_{a}-1)/2,-(n_{a}-3)/2,\ldots,(n_{a}-3)/2,(n_{a}-1)/2\}$ and

where $n=n_{a}+n_{b}$. For this choice, $D_{a}$ is an orthogonal Latin hypercube while $D_{b}$ is not. We examine how to construct an orthogonal design $D_{b}$ with level set $\mathcal{S}_{b}$ given in (6). Now consider the matrices in Table 2. Each of the four matrices in Table 2 has the following properties: (i) it has real entries $\pm x_{1},\ldots,\pm x_{n/2}$; (ii) both $x_{i}$ and $-x_{i}$ occur exactly once in each column; (iii) every two columns are orthogonal. We note that the matrices in Table 2 are related to but different from orthogonal designs in the combinatorics literature [Geramita and Seberry (1979)].

The matrices in Table 2 can be used to construct orthogonal Latin hypercubes of $n$ runs by setting $x_{i}=(2i-1)/2$ for $i=1,\ldots,n/2$. They also provide a direct construction of orthogonal designs $D_{b}$ with level set $\mathcal{S}_{b}$ in (6) by choosing $x_{i}=(n_{a}+2i-1)/2$ for $i=1,\ldots,n_{b}/2$. Most importantly, they are useful in the following result that allows us to construct $D_{b}$ with level set $\mathcal{S}_{b}$ in (6) for more general $n_{b}$.

Orthogonality of design $L$ follows from Lemma 2. That $L$ has a desired set of levels can easily be established which follows a similar path to that for Lemma 1. Comparing Proposition 2 with Theorem 1, we see that the only changes are those made to $\gamma$ and $C$. Mathematically, Theorem 1 is a special case of Proposition 2 as one can obtain the former from the latter by setting $n_{a}=0$. We present them separately because they carry different messages and serve different purposes in this paper.

Design $C$ required in Proposition 2 can easily be obtained from the matrices in Table 2. By letting $n=n_{b}$ in Proposition 2, design $L$ in Proposition 2 can then used as our $D_{b}$ as it has desired level set $\mathcal{S}_{b}$ in (6). The run size $n_{b}$ of such $D_{b}$ has form $n_{b}=8k$. Since there is no restriction in the run size $n_{a}$ of $D_{a}$, other than that $D_{a}$ is an orthogonal Latin hypercube, this second stacking method allows orthogonal Latin hypercubes of any run size $n\neq 4k+2$ to be constructed.

The methods in Sections 3.1 and 3.2 both build large orthogonal Latin hypercubes from small ones. To apply the methods, we need to find orthogonal Latin hypercubes with small runs. Various efficient algorithms can be helpful in this regard. Lin (2008) reported an algorithm adapted from that of Xu (2002). The key idea of the algorithm is to add columns sequentially to an existing design. To add a column, two operations, pairwise switch and exchange, are used. A pairwise switch switches a pair of distinct levels in a column. For a candidate column, the algorithm searches for all possible pairwise switches and makes the pairwise switch that achieves the best improvement. This search and pairwise switch procedure is repeated until an orthogonal Latin hypercube is found. An exchange replaces the candidate column by a randomly generated column. The exchange step is repeated at most $T_{1}$ (user-specified) times if no orthogonal Latin hypercube is obtained. Since the procedure relies on the initial random columns, the entire procedure is repeated $T_{2}$ times. Apart from the sequential idea, the efficiency of the algorithm benefits from its fast updates of orthogonality. An update is needed when a pairwise switch is applied. The maximum number $m$ of the columns in orthogonal Latin hypercubes of $n$ runs found by the algorithm is given in Table 3 for $4\leq n\leq 21$ except for $n=16$, in which case, our algorithm finds $m=6$. The entry $m=12$ for $n=16$ in Table 3 is due to Steinberg and Lin (2006). The detailed design matrices for the orthogonal Latin hypercubes in Table 3 are presented in Lin (2008) and also available from the authors.

For a concise presentation of the results in this section, we use $\operatorname{OLH}(n,m)$ to denote an orthogonal Latin hypercube of $n$ runs for $m$ factors. We now present a general result from the application of Theorem 1 in Section 3.1 and the second stacking method in Section 3.2.

We give a proof for Theorem 3. The proof in fact provides a detailed procedure for the actual construction of these orthogonal Latin hypercubes. Part (i) of Theorem 3 results from an application of Theorem 1 in Section 3.1. In the general construction (1), we choose $B$ to be the given $\operatorname{OLH}(n,m)$. Matrix $D$ is obtained by taking $m$ columns from a Hadamard matrix of order $n$. Design $C$ is chosen to be an orthogonal Latin hypercube derived from a matrix in Table 2. Note that each of the four matrices in Table 2 has a fold-over structure in that it can be written as $(X^{T},-X^{T})^{T}$. Now let $A=(S^{T},S^{T})^{T}$ where $S$ is obtained from $X$ by setting $x_{i}=1$ for all $i$. With the above choices for $A$, $B$, $C$ and $D$, conditions (i), (ii), (iii) and (iv) in Theorem 1 are all satisfied. This proves part (i) of Theorem 3. The proof for part (ii) of Theorem 3 is similar, involving the second stacking method with $n_{a}=1$ and an application of Proposition 2.

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

One important problem in the study of orthogonal Latin hypercubes is to determine the maximum number $m^{*}$ of factors for an $\operatorname{OLH}(n,m^{*})$ to exist. Theorem 2 says that $m^{*}=1$ if $n$ is 3 or has form $n=4k+2$ and that $m^{*}\geq 2$ otherwise. This result is now strengthened below.

Part (i) of Proposition 3 is obtained as follows. By our second stacking method with the use of the $16\times 8$ matrix in Table 2, we can construct an $\operatorname{OLH}(n+16,m)$ where $m\leq 8$ if an $\operatorname{OLH}(n,m)$ is available. Part (i) of Proposition 3 will be true if we can claim that an $\operatorname{OLH}(n,6)$ exists for all $17\leq n\leq 32$ except for $n=18,22,26,30$. We already know that the claim is true for $n=17,19,20,21$ from Table 3 and for $n=32$ from Example 2. Note that an $\operatorname{OLH}(11,6)$ can be obtained by choosing any six columns from the $\operatorname{OLH}(11,7)$ in Table 3. For $n=23$, we use the first stacking method by choosing $n_{a}=11$ and $n_{b}=12$ and using an $\operatorname{OLH}(11,6)$ and the $\operatorname{OLH}(12,6)$ in Table 3. The case $n=24$ follows from applying part (i) of Theorem 3 to the $\operatorname{OLH}(12,6)$ in Table 3. For $n=25$, an $\operatorname{OLH}(25,6)$ can be constructed using the first stacking method with $n_{a}=13$ and $n_{b}=12$. For $n=27$, we apply the second stacking method by choosing $n_{a}=11$ and $n_{b}=16$. The second stacking method also allows the construction of an $\operatorname{OLH}(28,6)$, an $\operatorname{OLH}(29,6)$ and an $\operatorname{OLH}(31,6)$. We choose $n_{a}=12$ and $n_{b}=16$ for $n=28$, $n_{a}=13$ and $n_{b}=16$ for $n=29$, and $n_{a}=15$ and $n_{b}=16$ for $n=31$. Part (ii) follows from the existence of an $\operatorname{OLH}(11,7)$ in Table 3. Parts (iii), (iv) and (v) follows from an application of Theorem 3.

The following remarks are in order regarding Proposition 3. If we wish, we can obtain sharper lower bounds on $m^{*}$ for certain values of $n$ by applying Theorem 3. For example, using the $\operatorname{OLH}(12,6)$ in Table 3, we can establish that $m^{*}\geq 6\times 8^{k}$ for $n=12\times 16^{k}$. We will not dwell further on this issue but are satisfied with the general lower bound in Proposition 3. The lower bound in Proposition 3 is derived from the small orthogonal Latin hypercubes found by our algorithm. Therefore, improved bounds will be naturally available in the future if better results are obtained from computer search.