Optimal Sliced Latin Hypercube Designs with Slices of Arbitrary Run Sizes

TL;DR

This paper introduces a new method for constructing optimal sliced Latin hypercube designs with arbitrary run sizes, improving space-filling properties and flexibility in experimental design.

Contribution

It proposes a novel construction method and a combined space-filling measurement for SLHDs with arbitrary run sizes, along with algorithms to find optimal designs.

Findings

Effective construction of SLHDs with arbitrary run sizes demonstrated

New space-filling measurement improves design quality

Algorithms successfully identify optimal designs in examples

Abstract

Sliced Latin hypercube designs (SLHDs) are widely used in computer experiments with both quantitative and qualitative factors and in batches. Optimal SLHDs achieve better space-filling property on the whole experimental region. However, most existing methods for constructing optimal SLHDs have restriction on the run sizes. In this paper, we propose a new method for constructing SLHDs with arbitrary run sizes, and a new combined space-filling measurement describing the space-filling property for both the whole design and its slices. Furthermore, we develop general algorithms to search the optimal SLHD with arbitrary run sizes under the proposed measurement. Examples are presented to illustrate that effectiveness of the proposed methods.

Tables

Table 1. Table 1: Performance of the efficient two-part algorithm for repeating 100 times

Algorithm	Design	Min	Mean	Max	SD	AT
SESE	FSLHD$(15,30;2,2)$	7.8674	8.3100	8.7988	0.0172	406 seconds
Part-I	FSLHD$(15,30;2,2)$	9.0118	10.3349	14.5465	0.5958	12 seconds
Part-I +Part-II	FSLHD$(15,30;2,2)$	8.2610	9.1712	11.2161	0.1276	18 seconds
SESE	FSLHD$(5,10,15,30;4,6)$	1.8614	2.0823	2.4874	0.0103	1110 seconds
Part-I	FSLHD$(5,10,15,30;4,6)$	2.0474	2.2545	2.8380	0.0122	51 seconds
Part-I +Part-II	FSLHD$(5,10,15,30;4,6)$	1.9041	2.2424	2.0394	0.0031	71 seconds

Equations

$$d_s^i=(m_s^i-{\epsilon }_s^i)/L,s=1,\cdots ,n_i,$$

$$\underset{\forall 1\le i,j\le n,i\ne j}{\min }d({\mathbf{x}}_i,{\mathbf{x}}_j),$$

$$d({\mathbf{x}}_i,{\mathbf{x}}_j)=d_{ij}={\left(\sum _{k=1}^q|x_{ik}-x_{jk}{|}^m\right)}^{1/m},\text{ m=1 or 2}.$$

$${\varphi }_t={\left(\sum _{1\le i<j\le n}(d_{ij}{)}^{-t}\right)}^{1/t},$$

$${\mathbf{D}}^{\ast }=\arg \underset{\mathbf{D}}{\min }{\varphi }_t(\mathbf{D}).$$

$${\varphi }_{\text{CSM}}(\mathbf{D})=w{\varphi }_t(\mathbf{D})+(1-w)\left(\sum _{i=1}^u{\lambda }_i{\varphi }_t({\mathbf{D}}^{(i)})\right),$$

$${\varphi }_{\text{CSM}}(\mathbf{D})=w{\varphi }_{{\text{CD}}_2}(\mathbf{D})+(1-w)\left(\sum _{i=1}^u{\lambda }_i{\varphi }_{{\text{CD}}_2}({\mathbf{D}}^{(i)})\right),$$

$${\varphi }_{{\text{CD}}_2}=$$

$${+\frac{1}{n^2}\sum _{i=1}^n\sum _{j=1}^n\sum _{k=1}^m\left(1+\frac{1}{2}|x_{ik}-0.5|+\frac{1}{2}|x_{jk}-0.5|-|x_{ik}-x_{jk}|\right))}^{1/2}$$

$$P(\mathbf{A})=\sum _{1\le i<j\le n}1\{d({\mathbf{a}}_i,{\mathbf{a}}_j)=0\},$$

$$\bm{M}=\left(\begin{array}[]{cccccccccc}54&12&24&42&60&30&6&18&48&36\\ 54&42&12&24&18&6&36&48&60&30\end{array}\right)^{T}.$$

$$\bm{M}^{1}=\left(\begin{array}[]{cccccccccc}4&3&1&2&2&1&3&4&4&2\\ 4&1&2&3&4&2&1&2&4&3\end{array}\right)^{T},$$

$$\bm{M}^{2}=\left(\begin{array}[]{cccccccccc}6&5&2&3&2&1&4&5&6&3\\ 6&2&3&5&6&3&1&2&5&4\end{array}\right)^{T}.$$

$$d^{\prime }({\mathbf{x}}_r,{\mathbf{x}}_v)=((d({\mathbf{x}}_r,{\mathbf{x}}_v){)}^m+h(r,s,k,v){)}^{1/m},$$

$$d^{\prime }({\mathbf{x}}_s,{\mathbf{x}}_v)=((d({\mathbf{x}}_r,{\mathbf{x}}_v){)}^m-h(r,s,k,v){)}^{1/m},$$

$${\varphi }_{\text{CSM}}^{\prime }(\mathbf{D})=w{\varphi }_t^{\prime }(\mathbf{D})+(1-w)\left(\sum _{i=1}^u{\lambda }_i{\varphi }_t^{\prime }({\mathbf{D}}^{(i)})\right).$$

$$\begin{array}{rl}{\varphi }_t^{\prime }(\mathbf{D})=& (({\varphi }_t(\mathbf{D}){)}^t+\sum _{1\le v\le n,v\ne r,s}\left(d^{{}^{\prime }}({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}\right)\\ & {+\sum _{1\le v\le n,v\ne r,s}\left(d^{{}^{\prime }}({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t}\right))}^{1/t},\end{array}$$

$$\begin{gathered} {\varphi }_t^{\prime }({\mathbf{D}}^{(i)})= \\ \{\begin{array}{rl}& {\varphi }_t({\mathbf{D}}^{(i)}),\text{if i=e};\\ & ({\left({\varphi }_t({\mathbf{D}}^{(i)})\right)}^t+\sum _{v\in J_e,v\ne r,s}(d^{\prime }({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t})\\ & {+\sum _{v\in J_e,v\ne r,s}(d^{\prime }({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t})}^{1/t},\text{if i=e}.\end{array} \end{gathered}$$

$$\begin{array}{rl}{\varphi }_t^{\prime }(\mathbf{D})=& (({\varphi }_t(\mathbf{D}){)}^t+\sum _{1\le v\le n,v\ne r,s}(d^{{}^{\prime }}({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t})\\ & {+\sum _{1\le v\le n,v\ne r,s}(d^{{}^{\prime }}({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_s,{\mathbf{x}}_v{)}^{-t}))}^{1/t},\end{array}$$

$$\phi_{t}^{\prime}(\bm{D}^{(i)})=\left\{\begin{split}&\phi_{t}(\bm{D}^{(i)}),&\text{$if$ $i\neq e,e^{\prime}$};\\ &\left(\left(\phi_{t}(\bm{D}^{(i)})\right)^{t}+\sum_{v\in J_{e},v\neq r}\left(d^{\prime}(\bm{x}_{r},\bm{x}_{v})^{-t}-d(\bm{x}_{r},\bm{x}_{v})^{-t}\right)\right)^{1/t},&\text{$if$ $i=e$};\\ &\left(\left(\phi_{t}(\bm{D}^{(i)})\right)^{t}+\sum_{v\in J_{e^{\prime}},v\neq s}\left(d^{\prime}(\bm{x}_{s},\bm{x}_{v})^{-t}-d(\bm{x}_{s},\bm{x}_{v})^{-t}\right)\right)^{1/t},&\text{$if$ $i=e^{\prime}$}.\end{split}\right.$$

$$\begin{array}{rl}{\varphi }_t^{\prime }(\mathbf{D})=& {\left({\left({\varphi }_t(\mathbf{D})\right)}^t+\sum _{1\le v\le n,v\ne r}\left(d^{{}^{\prime }}({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}-d({\mathbf{x}}_r,{\mathbf{x}}_v{)}^{-t}\right)\right)}^{1/t},\end{array}$$

$$\phi_{t}^{\prime}(\bm{D}^{(i)})=\left\{\begin{split}&\phi_{t}(\bm{D}^{(i)}),&\text{$if$ $i\neq e$};\\ &\left(\left(\phi_{t}(\bm{D}^{(i)})\right)^{t}+\sum_{v\in J_{e},v\neq r}\left(d^{\prime}(\bm{x}_{r},\bm{x}_{v})^{-t}-d(\bm{x}_{r},\bm{x}_{v})^{-t}\right)\right)^{1/t},&\text{$if$ $i=e$}.\end{split}\right.$$

Full text

Keywords: computer experiment; optimal design; space-filling design; maximin distance criterion

1 Introduction

Computer experiments are becoming increasingly significant in many fields, such as finite element analysis and computational fluid dynamics. Latin hypercube designs (LHDs) McKay et al. (1979) are widely used in computer experiments because of their optimal univariate uniformity. A design with $n$ runs and $q$ factors is called an LHD, if the design is projected onto any one dimension, there is precisely one point lying within one of the $n$ intervals $(0,1/n],(1/n,2/n],\cdots,((n-1)/n,1]$. Such an LHD is said to have optimal univariate uniformity. Sliced Latin hypercube designs (SLHDs) are LHDs that can be partitioned into some LHD slices Qian (2012), which means that the SLHDs have the optimal univariate uniformity for both the whole design and their slices. In He and Qian (2016), a central limit theorem for SLHDs is proposed. SLHDs are popular for computer experiments with both qualitative and quantitative variables; see Qian et al. (2008); Gang et al. (2009); Deng et al. (2017) and the references therein. Each slice of an SLHD can be used under one level-combination of the qualitative factors. However, the original SLHDs and almost all existing methods for constructing variants of SLHDs requires that the run sizes of each slice are equal; see Hwang et al. (2016); Yin et al. (2014); Xie et al. (2014); Yang et al. (2016).

An SLHD is called desirable if its design points are well spread out for both the whole design and its slices. Randomly generated SLHDs usually have a poor space-filling property in the entire experimental region, i.e., randomly generated SLHDs may not be desirable. There are a lot of methods that aim to improve the space-filling property of an SLHD. For instance, the method proposed by Huang et al. (2015) can be used to generate an optimal clustered-sliced Latin hypercube design (OCSLHD) which has good space-filling property in the whole experimental region. In a multi-fidelity computer experiment, each slice of an OCSLHD can be used for each accuracy of the computer code Huang et al. (2015). Generally, we want to use more design points for the lower-accuracy experiments than those of the higher-accuracy experiments, since the lower the accuracy is, the faster it runs Kennedy and O’Hagan (2000); Qian et al. (2006). However, a lot of existing method for constructing optimal SLHDs can only generate SLHDs with equal run sizes of each slices, e.g., Huang et al. (2015); Ba et al. (2015); Chen et al. (2016). To overcome this restriction, we need a method that can construct SLHDs with slices of arbitrary run sizes, and with a good space-filling property over the whole experimental region. For example, reference Kong et al. (2018) gave flexible sliced designs, but such designs are not LHDs. The method given in Xu et al. (2018) provided SLHDs with unequal batch sizes, but this type of design only accommodates two different run sizes. An algorithm is proposed in Xu et al. (2019) to construct an SLHD with unequal run sizes, but this construction method is difficult to search the optimal design.

In this paper, we propose a method to construct SLHDs with slices of arbitrary run sizes, which are called flexible sliced Latin hypercube designs (FSLHDs). The new construction method can be easily adapted to generate the optimal design. Furthermore, we provide a combined space-filling measurement (CSM) to descibe the space-filling properties of both the whole design and each of slices. Based on an optimization algorithm called the enhanced stochastic evolutionary algorithm (ESE), we propose a sliced ESE (SESE) algorithm to find the optimal FSLHDs. We further develop an efficient two-part algorithm to improve the efficiency in generating space-filling FSLHDs with large runs and factors. The generated optimal FSLHDs have three attractive features: (i) arbitrary run sizes of all slices, (ii) optimal univariate uniformity in the whole design and each slice, (iii) good space-filling property in the experimental region. We believe that they are suitable for many multi-fidelity computer experiments in practice.

The remainder of this paper is organized as follows. The construction of FSLHDs is provided in Section 2. In Section 3, an CSM are given to descibe the space-filling properties of both the whole design and each of slices, and then we develop an SESE algorithm to obtain optimal FSLHDs based on the CSM and a two-part algorithm to improve efficiency. Some simulation results are illustrated in Section 4. Section 5 provides some discussions. Section 6 concludes this paper.

2 Construction of SLHDs with slices of arbitrary run sizes

For a real number $a$, let $\lceil a\rceil$ denote the smallest integer not smaller than $a$. Given $u$ positive integers $n_{1},\cdots,n_{u}$, let $n=\sum_{i=1}^{u}n_{i}$ and let $L={\rm lcm}(n_{1},\cdots,n_{u},n)$ be the least common multiple of $n_{1},\cdots,n_{u}$, and $n$. Suppose that ${\rm FSLHD}(n_{1},\cdots,n_{u};u,q)$ is an FSHLD with $u$ slices of run sizes $n_{1},\cdots,n_{u}$ and $q$ factors. Each column of the FSLHD is generated independently by the following algorithm.

1. Step 1.

Let $H^{i}=\varnothing$ for $i=1,\cdots,u$, and $R_{0}=\varnothing$. 2. Step 2.

For $j=1,\cdots,n$, let $R_{j,0}=R_{j-1}\bigcup\{j\}$ and calculate

$\theta_{j}=\sum_{i=1}^{u}\left(\lceil n_{i}(j+1)/n\rceil-\lceil n_{i}j/n\rceil\right)$.

If $\theta_{j}>0$, for $k=1,\cdots,\theta_{j}$, let $l$ denote the $k$th smallest integer of the set $\{p|\lceil n_{p}(j+1)/n\rceil-\lceil n_{p}j/n\rceil=1\}$ and $r=\min\{r|\lceil n_{l}r/n\rceil=\lceil n_{l}j/n\rceil,r\in R_{j,k-1}\}$ add $r$ to $H^{l}$ and let $R_{j,k}=R_{j,k-1}\setminus\{r\}$. Let $R_{j}=R_{j,\theta_{j}}$ and go to the next $j$. 3. Step 3.

For $i=1,\cdots,u$, generate a vector $\bm{h}^{i}$ by randomly permuting $H^{i}$. 4. Step 4.

For $i=1,\cdots,u$, calculate $\bm{m}^{i}=L\bm{h}^{i}/n$, where $\bm{m}^{i}=(m^{i}_{1},\cdots,m^{i}_{n_{i}})$. Combine $\bm{m}^{1},\cdots,\bm{m}^{u}$ to obtain an $n$-dimensional column vector $\bm{m}=(m_{1},\cdots,m_{n})^{T}$, then let $\bm{d}^{i}=(d^{i}_{1},\cdots,d^{i}_{n_{i}})^{T}$ be constructed by

[TABLE]

where $\varepsilon^{i}_{s}\sim U(0,1)$. Combine $\bm{d}^{1},\cdots,\bm{d}^{u}$ to obtain an $n$-dimensional column vector $\bm{d}=(d_{1},\cdots,d_{n})^{T}$, and $\bm{d}$ is one column of the design.

In the above algorithm $\bm{m}$ is called a column of the flexible sliced Latin hypercube (FSLH). The following theorem shows that both the whole FSLHD and its slices are LHDs.

Theorem 1

*Let $\bm{d}=(d_{1},\cdots,d_{n})^{T}$ denote an arbitrary column of FSLHD($n_{1},\cdots,n_{u};u,q$) generated by the above method. Let $\bm{d}^{1},\cdots,\bm{d}^{u}$ denote each slice. For $i=1,\cdots,u$, let $t^{i}=L/n_{i}$ and $t^{\prime}=L/n$.

(i)Precisely one point of $\bm{d}=(d_{1},\cdots,d_{n})^{T}$ lies within one of the $n$ intervals $(0,1/n],(1/n,2/n],\cdots,((n-1)/n,1]$.

(ii)Precisely one point of $\bm{d}^{i}=(d^{i}_{1},\cdots,d^{i}_{n_{i}})$ lies within one of the $n_{i}$ intervals $(0,1/n_{i}],(1/n_{i},2/n_{i}]$,$\cdots,((n_{i}-1)/n_{i},1]$.*

Proof 1

(i) Combine $\bm{h}^{1},\cdots,\bm{h}^{u}$ to obtain $\bm{h}=(h_{1},\cdots,h_{n})^{T}$ that is a permutation of $\{1,\cdots,n\}$. Combine $\bm{m}^{1},\cdots,\bm{m}^{u}$ to obtain $\bm{m}=(m_{1},\cdots,m_{n})^{T}$. Therefore, $\bm{m}=L\bm{h}/n$. For $t^{\prime}=L/n$, because $\lceil\bm{m}/t^{\prime}\rceil=\lceil(L\bm{h}/n)/t^{\prime}\rceil=\lceil\bm{h}\rceil$, $\lceil\bm{m}/t^{\prime}\rceil$ is a permutation of $\{1,\cdots,n\}$. Therefore, precisely one point of $\bm{d}=(d_{1},\cdots,d_{n})^{T}$ lies within one of the $n$ intervals $(0,1/n],(1/n,2/n],\cdots,((n-1)/n,1]$.

(ii) According to Step 2, for $i=1,\cdots,u$, it is clear that $card(H^{i})=\sum_{j=1}^{n}(\lceil n_{i}(j+1)/n\rceil-\lceil n_{i}j/n\rceil)=\lceil n_{i}(n+1)/n\rceil-\lceil n_{i}/n\rceil=n_{i}$, and for $j=1,\cdots,n$, $\lceil n_{i}j/n\rceil<\lceil n_{i}(j+1)/n\rceil$. For any $i$,$j$, there is an integer $h\in H^{i}$ that satisfies $\lceil n_{i}h/n\rceil=\lceil n_{i}j/n\rceil$. Therefore, we have $\{m|m=\lceil n_{i}h/n\rceil,h\in H^{i}\}=\{1,\cdots,n_{i}\}$, which means that $\lceil n_{i}\bm{h}^{i}/n\rceil$ is a permutation of $\{1,\cdots,n_{i}\}$. Since $\bm{m}^{i}=L\bm{h}^{i}/n$, we have $\lceil\bm{m}^{i}/t^{i}\rceil=\lceil(L\bm{h}^{i}/n)/(L/n_{i})\rceil=\lceil n_{i}\bm{h}^{i}/n\rceil$. Thus, $\lceil\bm{m}^{i}/t^{i}\rceil$ is a permutation of $\{1,\cdots,n_{i}\}$. Therefore, precisely one point of $\bm{d}^{i}$ lies within one of the $n_{i}$ intervals $(0,1/n_{i}],(1/n_{i},2/n_{i}],\cdots,((n_{i}-1)/n_{i},1]$.

We give an example to illustrate the process of the above method.

Example 2.1

Consider $n_{1}=3$, $n_{2}=4$, $n_{3}=5$, $u=3$, $n=12$, and $L=60$.

1. Step 1.

$H^{1}=H^{2}=H^{3}=R_{0}=\emptyset$. 2. Step 2.

Calculate $(\theta_{1},\cdots,\theta_{n})=(0,1,1,2,0,1,1,1,2,0,0,3)$. For $j=1$, then $R_{1,0}=\{1\}$, since $\theta_{1}=0$, we obtain $R_{1}=R_{1,0}=\{1\}$. For $j=2$, $R_{2,0}=R_{1}\cup\{2\}=\{1,2\}$, $\theta_{2}=1$, only an integer $l=3$ satisfies $\lceil n_{l}(j+1)/n\rceil-\lceil n_{L}j/n\rceil=1$, and $r=$min$\{r|\lceil n_{3}r/n\rceil=\lceil n_{3}j/n\rceil$, $r\in R_{2,0}\}=$ min$\{1,2\}=1$. Hence, we add $r=1$ to $H^{3}$, $R_{2,1}=R_{2,0}\backslash\{1\}=\{2\}$, and $R_{2}=R_{2,1}=\{2\}$. For $j=3$, $R_{3,0}=R_{2}\cup\{3\}=\{2,3\}$, $\theta_{3}=1$, only an integer $l=2$ satisfies $\lceil n_{2}(j+1)/n\rceil-\lceil n_{2}j/n\rceil=1$, and $r=$min$\{r|\lceil n_{2}r/n\rceil=\lceil n_{2}j/n\rceil$, $r\in R_{3,0}\}=$ min$\{2,3\}=2$. Therefore, we add $r=2$ to $H^{2}$ , $R_{3,1}=R_{3,0}\backslash\{2\}=\{3\}$, and $R_{3}=R_{3,1}=\{3\}$. After passing all $j$, we can get $R_{12}=\emptyset$, $H^{1}=\{3,7,10\}$, $H^{2}=\{2,5,8,11\}$, and $H^{3}=\{1,4,6,9,12\}$. 3. Step 3.

We get $\bm{h}^{1}=(10,7,3)$ , $\bm{h}^{2}=(5,8,2,11)$, and $\bm{h}^{3}=(6,9,12,1,4)$ by randomly permuting $H^{1}$ , $H^{2}$, and $H^{3}$. 4. Step 4.

We obtain $\bm{m}^{1}=(50,35,15)$, $\bm{m}^{2}=(25,40,10,55)$, and $\bm{m}^{3}=(30,45,60,5,20)$. Then $\bm{d}^{i}=(d^{i}_{1},\cdots,d^{i}_{n_{i}})^{T}$ is constructed through $d^{i}_{s}=(m^{i}_{s}-\varepsilon^{i}_{s})/60$, where $i=1,\cdots,3$, $s=1,\cdots,n_{i}$, and $\varepsilon^{i}_{s}\sim U(0,1)$. Thus, we obtain an arbitrary column $\bm{d}=(d_{1},\cdots,d_{n})^{T}$ of the design.

3 Optimal SLHDs with slices of arbitrary run sizes

Given $n_{1},\cdots,n_{u},u,q$, a number of possible FSLHDs can be generated through the proposed method in Section 2. Among such FSLHDs, we can find the optimal FSLHD through a given space-filling criterion. We first propose a combined space-filling measurement (CSM) to evaluate space-filling property of FSLHD in Subsection 3.1. Then, to keep the structure of the design during the optimization process, three methods are proposed to change position of the elements in one column in Subsection 3.2. Finally, we present a sliced ESE algorithm to optimize FSLHD in Subsection 3.3. An efficient two-part algorithm for generating the space-filling FSLHD is given in Subsection 3.4.

3.1 A combined space-filling measurement for FSLHDs

Various space-filling criteria are used to evaluate the LHDs, such as the maximin distance criterion Johnson et al. (1990); Grosso et al. (2009); Dam et al. (2007, 2009), the $\phi_{t}$ criterion Jin et al. (2016); Morris and Mitchell (1995); Ye et al. (2000); Viana et al. (2010), and the centered $L_{2}$-discrepancy (${\textrm{C}D}_{2}$) criterion Hickernell (1998); Fang et al. (2002). All the space-filling criteria can be extend to describe space-filling propert of the FSLHDs. We mainly focus on the $\phi_{t}$ criterion which is an attractive extension of maximin distance criterion.

The maximin distance criterion is a popular space-filling criterion introduced in Johnson et al. (1990). Let $\bm{D}=[\bm{x}_{1},\cdots,\bm{x}_{n}]^{T}$ denote a design matrix with $n$ runs and $q$ factors, where each row $\bm{x}_{i}^{T}=(x_{i1},\cdots,x_{iq})$ is a design point and each column is a factor with $i=1,\cdots,n$. A maximin distance design is generated by maximizing the minimum inter-site distance, which is expressed as

[TABLE]

where $d(\bm{x}_{i},\bm{x}_{j})$ is the distance between the design points $\bm{x}_{i}$ and $\bm{x}_{j}$ given by:

[TABLE]

Here $m=1$ and $m=2$ are the rectangular and Euclidean distances, respectively. In this article, we use the Euclidean distance. An extension of the maximin distance criterion Jin et al. (2016) is given by

[TABLE]

where $t$ is a positive integer. It is obviously that as $t\longrightarrow\infty$, minimizing (4) is equivalent to maximizing (2). The calculation of $\phi_{t}$ is simpler compared with the maximin distance criterion.

We search an optimal design by minimizing $\phi_{t}$, i.e.

[TABLE]

Suppose that $\bm{D}$ is the design matrix of an ${\rm FSLHD}(n_{1},\ldots,n_{u};u,q)$. For $i=1,\ldots,u$, let $\bm{D}^{(i)}$ denote each slice of $\bm{D}$. We need to consider both the space-filling properties of the whole FSLHD and that of its slices. Consequently, our goal is to find a maximin FSLHD that minimizes $\phi_{t}(\bm{D})$ for the entire design as well as $\phi_{t}(\bm{D}^{(i)})$ for each slice of $\bm{D}$ ($i=1,\ldots,u$). This is a multi-objective optimization problem. It is a common method in multi-objective problem to use a weighted average of all individual objectives. It motivates us to develop a combined space-filling measurement (CSM) as follows:

[TABLE]

where $\lambda_{i}=n_{i}/n$ , $\sum_{i=1}^{u}\lambda_{i}=1$, and $w\in(0,1)$. Since run sizes of slices are $n_{1},\ldots,n_{u}$, respectively, it makes sense that we take the weight of each slice to be $\lambda_{i}=n_{i}/n$, for $i=1,\ldots,u$. The weight $w$ is selected flexibly. The space-filling property of the whole FSLHD is more important, hence we set $w=1/2$ in general. We can define a maximin distance FSLHD with respect to the CSM as the one which minimizes (6).

Note that other space-filing criteria can also evaluate the FSLHD. For instance, we can obtain an uniform FSLHD by minimizing a similar CSM given by

[TABLE]

where $\phi_{\textrm{CD}_{2}}$ is the centered $L_{2}$-discrepancy defined as

[TABLE]

proposed in Hickernell (1998).

3.2 Exchange procedures for FSLHDs

In the literature, some optimization algorithms have widely used to construct an optimal LHD. They utilize an exchange procedure to iteratively search the optimal LHD in the design space. In this way, two randomly selected elements in an arbitrary column of an LHD are exchanged to generate a new design. The exchange procedure for an FSLHD is more complex since the design should keep the sliced structure. In this subsection, in the optimization process of an FSLHD, we present three exchange procedures to generate a neighbor of the design which do not change the sliced structure of the design. A neighbor of an FSLH corresponds to a neighbor of an FSLHD. Let $\bm{M}$ be the FSLH($n_{1},\cdots,n_{u};u,q$) constructed in Section 2. Let $\bm{M}_{\textrm{N}}$ denote a neighbor of an FSLH and let $\bm{D}_{\textrm{N}}$ denote a neighbor of an FSLHD.

3.2.1 The within-slice exchange procedure

Given an FSLH($n_{1},\cdots,n_{u};u,q$)$(\bm{M})$, let $n_{0}=0$, $r_{0}=0$, and $r_{i}=\sum_{k=0}^{i}n_{k}$, for $i=1,\cdots,u$. The within-slice exchange procedure in the $i$th slice of $\bm{M}$ is to draw an $\bm{M}_{\textrm{N}}$ by the following four steps:

1. Step 1.

Randomly select a column of $\bm{M}$. 2. Step 2.

Select any two different elements $d_{j},d_{k}$ in $i$th slice of the column, where $r_{i-1}+1\leq j,k\leq r_{i}$. 3. Step 3.

Exchange $d_{j}$ and $d_{k}$ in the same slice. 4. Step 4.

Generate $\bm{M}_{\textrm{N}}$.

After this procedure, the neighbor design $\bm{M}_{\textrm{N}}$ still keeps the sliced structure. The within-slice exchange procedure is explained by an example about FSLH(4,6;2,2) illustrated in Figure 1.

3.2.2 The different-slice exchange and the out-slice exchange procedures

We first give some notations. Given an FSLH($n_{1},\cdots,n_{u};u,q$)$(\bm{M})$, let $\bm{M}(l:m,j)$ denote the $l$th to $m$th rows of the $j$th column, and $\bm{M}(l,j)$ denotes its $(l,j)$ element. For $i=1,\cdots,u$, let $r_{i}=\sum_{k=0}^{i}n_{k}$, and $B_{ij}=\bm{M}(r_{i-1}+1:r_{i},j)$ denotes $i$th slice in the $j$th column of $\bm{M}$ , where $n_{0}=0$, $r_{0}=0$. Define $E_{ij}=\{\bm{M}(r_{i}+1,j),\cdots,\bm{M}(n,j)\}$, where $i=1,\cdots,u-1$ and $n=\sum_{i=1}^{u}n_{i}$. Let $A_{L}=\{1,\cdots,L\}$ denote a set of integers from 1 to $L$, where $L=$ lcm$(n_{1},\cdots,n_{u},n)$. Set $B=\{\bm{M}(1,j),\cdots,\bm{M}(n,j)\}$. Let $C=A_{L}\setminus B$ denote $A$ minus $B$.

It is observed that elements of each slice on an FSLHD are fixed by the construction method in Section 2. There are two situations. On the one hand, some elements in an arbitrary column of an FSLH from different slices are exchanged, and the resulting FSLH does not change the sliced structure. On the other hand, some elements which are used to construct an column of an FSLH are not selected in $C$, besides, we exchange some elements between $B_{ij}$ and $C$, and the resulting FSLH still keeps the sliced structure. It motivates us to propose a different-slice exchange procedure and an out-slice exchange procedure to generate more diverse neighbors of the design. By the above ways, we can more easily find the optimal design. The detailed process of the two procedures are as follows.

The different-slice exchange procedure in the $i$th slice: we select any element $b$ of $B_{ij}$. Let $\rho(b)$ be a subset of $E_{ij}$ satisfying that the generated FSLH still keeps the sliced structure by exchanging $b$ with arbitrary $c$ in $\rho(b)$, where $i=1,\cdots,u-1$.

The out-slice exchange procedure in the $i$th slice: the elements in $C$ are called out-slice elements in a column of the design. For the same $b$, let $\sigma(b)$ be a subset of $C$ satisfying that the obtained FSLH still maintains the sliced structure through exchanging $b$ with arbitrary $c$ in $\sigma(b)$, where $i=1,\cdots,u$. Let $\tau(b)=\rho(b)\cup\sigma(b)$. In the last slice, we only consider the out-slice exchange procedure, thus $\tau(b)=\sigma(b)$. For a set $R$, $R_{k}$ denotes the $k$th smallest element of $R$. Suppose that $\bm{M}_{\textrm{N}}(1:n,j)$ is a new column generated from $\bm{M}(1:n,j)$. Here, for $i=1,\cdots,u$, recall that $t^{i}=L/n_{i}$. We provide a method to generate $\tau(b)$ in the $i$th slice of $\bm{M}$ by the following steps:

1. Step 1.

Randomly select an element $b$ in $\bm{M}(r_{i-1}+1:r_{i},j)$. 2. Step 2.

Generate a set $R=\{(\lceil b/t^{i}\rceil-1)\times t^{i}+1,(\lceil b/t^{i}\rceil-1)\times t^{i}+2,\cdots,\lceil b/t^{i}\rceil\times t^{i}\}\backslash\{b\}$. 3. Step 3.

If $i<u$, go to Step 4; else, go to Step 5. 4. Step 4.

For $k$ from $1$ to $t^{i}-1$, if $R_{k}$ belongs to $\bm{M}(r_{i}+1:n,j)$, go to Step 5; else, go to Step 6. 5. Step 5.

Generate $\bm{M}_{\textrm{N}}(1:n,j)$ by exchanging $b$ with $R_{k}$. If $\bm{M}_{\textrm{N}}(1:n,j)$ still satisfies Theorem 1(ii), go to Step 7. 6. Step 6.

Generate $\bm{M}_{\textrm{N}}(1:n,j)$ by exchanging $b$ with $R_{k}$. If $\bm{M}_{\textrm{N}}(1:n,j)$ still satisfies Theorem 1(i), go to Step 7. 7. Step 7.

Add $R_{k}$ to $\tau(b)$.

Step 5 and Step 6 are critical for generating $\tau(b)$. In Step 5, since both $b$ and $R_{k}$ are in $\bm{M}(1:n,j)$, $\bm{M}_{\textrm{N}}(1:n,j)$ still satisfies Theorem (i), when we exchange $b$ with $R_{k}$. Thus, we just guarantee that $\bm{M}_{\textrm{N}}(1:n,j)$ still satisfies Theorem 1(ii). In Step 6, it is clear that changing $b$ with any element of $R$ can guarantee that $\bm{M}_{\textrm{N}}(1:n,j)$ still satisfies Theorem 1(ii), therefore, we only ensure that $\bm{M}_{\textrm{N}}(1:n,j)$ satisfies Theorem 1(i).

We introduce the different-slice exchange and the out-slice exchange procedures in Figure 2. For an FSLH(4,6; 2,2)$(\bm{M})$, we randomly select $b=54$ in $\bm{M}(1:4,1)$ in Figure 2(a) , then $t^{1}=60/4=15$ and $R=\{45,46,\cdots,53,55,\cdots,60\}$. We obtain $\tau(b)=\{49,50,51,52,53,60\}$ after conducting the above steps. In the different-slice exchange procedure, we can exchange 54 with 60 of $\tau(b)$ in Figure 2(a). In the out-slice exchange procedure, we can replace 54 with 49 of $\tau(b)$ in Figure 2(b). It can be seen that the two resulting designs still keep the sliced structure.

3.3 A sliced ESE algorithm for generating optimal FSLHDs

Researchers utilize various optimization algorithms to construct optimal LHDs, such as the enhanced stochastic evolutionary (ESE) algorithm Jin et al. (2016), the simulated annealing search algorithm Morris and Mitchell (1995), the column wise-pairwise swap algorithm Ye et al. (2000), the threshold accepting algorithm Fang et al. (2002), the particle swarm algorithm Chen et al. (2013); Kennedy and Eberhart (1995), and the genetic algorithm Liefvendahl and Stocki (2006); Bates et al. (2004). All the above algorithms can be extended to optimize FSLHDs. In this paper, we choose the ESE algorithm as a basic algorithm to find optimal FSLHDs.

The ESE algorithm can quickly construct an optimal LHD in a limited calculative resource and it can also move from a locally optimal LHD. The ESE algorithm includes double loops, i.e., an inner loop and an outer loop. The inner loop randomly generates neighbors of the design by the exchange procedures and decides whether to accept them on the basis of an acceptance criterion. The outer loop aims to adjust the threshold $T_{\textrm{h}}$ in the acceptance criterion through the performance of the inner loop, so the outer loop can control the whole optimization process. When extending the ESE algorithm for searching an optimal FSLHD, we need to consider the sliced structure of an FSLHD. Thus, based on the three exchange procedures in Section 3.2, we develop a sliced enhanced stochastic evolutionary (SESE) algorithm which contains double loops in Jin et al. (2016) and the slice by slice loop proposed in this article. Such a combined algorithm can suit the sliced structure of the FSLHD. It is a dynamic optimization approach to optimize the FSLHD slice by slice. This algorithm can search the optimal FSLHD by minimizing the CSM. Algorithm 1 describes the SESE algorithm.

The slice by slice loop: we start with an initial FSLHD denoted by $\bm{D}_{0}$. When we optimize the first slice of the design, $\bm{D}_{0}$ is an initial design in the outer loop. When optimizing the $i$th ($i\geq 2$) slice of the design, we make $\bm{D}_{\textrm{best}}$, generated from outer loop in the $(i-1)$th slice optimization, as the initial FSLHD. It means that a new slice optimization is based on the previous slice optimization until the last slice. The parameter settings of the inner loop and the outer loop have been discussed in Jin et al. (2016). The parameter settings are similar in Jin et al. (2016) for the construction method of an FSLHD.

The inner loop: the iterations $P$ should be set larger for larger problems but no larger than 100. The acceptance criterion is $\phi_{\textrm{CSM}}(\bm{D}_{\textrm{N}})-\phi_{\textrm{CSM}}(\bm{D})\leq T_{h}\cdot random(0,1)$, where $random(0,1)$ generates uniform numbers between 0 and 1. According to the discussion in Jin et al. (2016), if the settings of $I_{1},I_{2}$, and $I_{3}$ are too large, it can appear the locally optimal design for designs with small run sizes and low efficiency for designs with large run sizes. Let $I_{1}=$ min $(n_{\textrm{in-slice}}/5,50)$, where $n_{\textrm{in-slice}}$ is the number of all possible neighbors of the design in within-slice exchange procedure. Let $n_{\textrm{diff-slice}}$ and $n_{\textrm{out-slice}}$ be the number of all possible neighbors of the design for the different-slice exchange procedure and the out-slice exchange procedure, respectively. According to the construction method of the FSLHD, we can clearly know that $n_{\textrm{diff-slice}}$ and $n_{\textrm{out-slice}}$ are usually small, therefore it is reasonable to set $I_{2}+I_{3}=$ min$(n_{\textrm{diff-slice}}+n_{\textrm{out-slice}},50)$.

The outer loop: The setting of $T_{\textrm{h}}$ is a small value, i.e., $T_{\textrm{h}_{0}}=0.005\times$ (criterion value of the initial design). The threshold $T_{\textrm{h}}$ is adjusted by an improving process and an exploration process. After the Inner Loop, if the search process has improvement, then go to the improving process, while if the search process has no improvement, then go to the exploration process. We adjust $T_{\textrm{h}}$ by the same way in Jin et al. (2016) as follows. In the improving process, when $T_{\textrm{h}}$ keeps on a small value, only slightly worse design or better design will be accepted. The parameter $P$ is the number of tries in the inner loop. The threshold $T_{\textrm{h}}$ is adjusted by the acceptance ratio $p_{\textrm{ac}}=n_{\textrm{ac}}/P$ ($n_{\textrm{ac}}$, the number of the accepted designs) and the improvement ratio $p_{\textrm{im}}=n_{\textrm{im}}/P$ ($n_{\textrm{im}}$, the number of the improved designs). For $flag_{\textrm{im}}=1$, if $p_{\textrm{ac}}>0.1$ and $p_{\textrm{im}}<p_{\textrm{ac}}$, let $T_{\textrm{h\_try}}=\beta_{1}T_{h}$, where $0<\beta_{1}<1$; if $p_{\textrm{ac}}>0.1$ and $p_{\textrm{im}}=p_{\textrm{ac}}$, let $T_{\textrm{h\_try}}=T_{\textrm{h}}$; otherwise, $T_{\textrm{h\_try}}=T_{\textrm{h}}/\beta_{1}$. We set $\beta_{1}=0.8$, since it appears to do well in all tests. In the exploration process, $T_{\textrm{h}}$ is adjusted by $p_{\textrm{ac}}$. For $flag_{\textrm{im}}=0$, let $T_{\textrm{h\_try}}=T_{\textrm{h}}/\beta_{2}$ and $T_{\textrm{h}}$ will be quickly increased until $p_{\textrm{ac}}>0.8$; if $p_{\textrm{ac}}>0.8$, let $T_{\textrm{h\_try}}=T_{\textrm{h}}\beta_{3}$ and $T_{\textrm{h}}$ will be quickly decreased until $p_{\textrm{ac}}<0.1$, where $0<\beta_{2},\beta_{3}<1$. On the basis of some tests, the settings of $\beta_{2}=0.7$ and $\beta_{3}=0.9$ perform well. Increasing rapidly $T_{\textrm{h}}$ (more worse designs can be accepted) is useful to go away from a locally optimal design. After going away from a locally optimal design, decreasing slowly $T_{\textrm{h}}$ helps to search better designs. An improved design is found by repeating the exploration process, then we go into the improving process. The $tol$ is a small fixed value, i.e., $tol=0.1$. The stopping criterion $N$ is set to be 10 in our procedure, which is selected flexibly.

3.4 Efficient two-part algorithm for generating space-filing FSLHDs

For an FSLHD with $n$ runs and $q$ factors, when $n$ and $q$ are small, the SESE algorithm is more efficient and provides much better resulting designs. However, if $n$ and $q$ are getting larger, the convergence of the SESE algorithm may be slow because of the large number of neighbors of the design. In this subsection, we consider a similar strategy which is broadly applied in Ba et al. (2015); Chen and Xiong (2017) to avoid the poor space-filling designs and improve the efficiency when $n$ and $q$ are large.

We first give the strategy for our proposed design as follows: for an FSLHD$(n_{1},\cdots,n_{u};u,q)$ and $n=\sum_{i=1}^{u}n_{i}$, the $q$-dimensional input region in the $i$th slice of FSLHD is partitioned into $n_{i}^{q}$ cells through the $\underbrace{n_{i}\times,\cdots,\times n_{i}}_{q}$ coarser grid ($i=1,\cdots,u$). Since run sizes $n_{i}$ of each slice are different, the number $n_{i}^{q}$ of divided cells is different. It is possible that some of $n$ design points sampled from the $n_{i}^{q}$ cells can fall into the same cell. If $n_{i}^{q}>n$, we need to avoid design points falling into in the same cells and ensure the design still an FSLHD.

We give a detailed process of the above strategy. Let $\bm{1}\{\cdot\}$ denote the indicator function. For an $n\times q$ matrix $\bm{A}=[\bm{a}_{1},\cdots,\bm{a}_{n}]^{\textrm{T}}$, denote

[TABLE]

where $\bm{1}\{d(\bm{a}_{i},\bm{a}_{j})=0\}=1$ if $d(\bm{a}_{i},\bm{a}_{j})=0$ is true and $\bm{1}\{d(\bm{a}_{i},\bm{a}_{j})=0\}=0$ otherwise. It is clear that some rows of matrix $A$ are the same if $P(\bm{A})>0$. We call the same rows as repeating rows which fall into the same cell. We can find repeating rows of a design by (9). For FSLH($n_{1},\cdots,n_{u};u,q$) ($\bm{{M}}$), recall that $t^{i}=$lcm $(n_{1},\cdots,n_{u},n)/n_{i}$, for $i=,\dots,u-1$. Let $\bm{M}^{i}=\lceil\bm{M}/t_{i}\rceil$. If $n_{i}^{q}>n$ and $P(\bm{M}^{i})=\sum_{1\leq i<j<n}\textbf{1}\{d(\bm{a}_{i},\bm{a}_{j})=0\}>0$, then the matrix has repeating rows.

Let us look at the following example of a design matrix FSLH (4,6;2,2) ($\bm{M}$)

[TABLE]

By (10) and (11), both $\bm{M}^{1}=\lceil\bm{M}/15\rceil$ and $\bm{M}^{2}=\lceil\bm{M}/10\rceil$ have repeating rows, which indicates that $P(\bm{M}^{1})=4>0$ and $P(\bm{M}^{2})=1>0$. The FSLH corresponding to the design under different divided cells is depicted in Figure 3 and Figure 3, respectively. The design points of repeating rows fall into the same cell (filled with blue).

[TABLE]

[TABLE]

To make the design with better space-filling properties, we consider to put all the points into the different cells. Therefore, we can select randomly a column of the repeating rows, and conduct a within-slice exchange procedure in the randomly chosen column of the same slice, until $P(\bm{M}^{1})=0$ and $P(\bm{M}^{2})=0$. The resulting design are shown in Figure 4 and Figure 4, respectively, in which all the points fall into the different cells. In summary, the above strategy can quickly eliminate the undesirable designs that contain repeating rows.

Given an FSLHD with large $n$ runs and $q$ factors, we develop an efficient two-part algorithm for finding the space-filling FSLHDs based on the above strategy. Without loss of generality, assume $n_{1},\cdots,n_{u}$ with $n_{1}\leq n_{2}\leq,\cdots,\leq n_{u}$. Recall that $\bm{D}_{\textrm{N}}$ denotes a neighbor of FSLHD($\bm{D}$) and $\bm{M}_{\textrm{N}}$ denotes a neighbor of FSLH($\bm{M}$). This algorithm is provided as follows:

Part-I algorithm

The Part-I algorithm is useful for speeding up by removing some undesirable designs from neighbors of the design. It starts with an initial FSLH ($n_{1},\cdots,n_{u};u,q$)$(\bm{M}_{0})$. According to the run sizes of the design, it can be stopped by some flexible stopping criterions. In our proposed algorithm, when 100 iterations have been operated, we stop the program. The algorithm is given below:

1. Step 1.

Let $\bm{M}=\bm{M}_{0}$, and set the index $i=1$. 2. Step 2.

If $P(\lceil\bm{M}/t_{i}\rceil)=0$, compute $\phi_{\textrm{CSM}}(\bm{D})$, go to Step 5. 3. Step 3.

If $n_{i}^{q}>n$, randomly choose a repeating row of $\lceil\bm{M}/t_{i}\rceil$, and randomly choose another row in the same slice. We exchange two elements which corresponds to a randomly selected column of the two rows. Generate an $\bm{M}_{\textrm{N}}$; else, go to Step 5. 4. Step 4.

If $P(\lceil\bm{M}_{\textrm{N}}/t_{i}\rceil)<P(\lceil\bm{M}/t_{i}\rceil)$, $\bm{M}=\bm{M}_{\textrm{N}}$, go back to Step 2; else, go back to Step 3. 5. Step 5.

Under the condition of $P(\lceil\bm{M}_{\textrm{N}}/t_{i}\rceil)=0$, generate an $\bm{M}_{\textrm{N}}$ by the within-slice procedure in the $i$th slice of $\bm{M}$, then calculate $\phi_{\textrm{CSM}}(\bm{D}_{\textrm{N}})$. 6. Step 6.

If $\phi_{\textrm{CSM}}(\bm{D}_{\textrm{N}})<\phi_{\textrm{CSM}}(\bm{D})$, then replace $\bm{M}$ by $\bm{M}_{\textrm{N}}$; else, go back to Step 3. 7. Step 7.

Repeat Step 4 and Step 5 until meeting the stopping criterion. 8. Step 8.

Update $i=i+1$, if $i<u$, go to Step 2; else, output $\bm{M}_{\textrm{best}}$ = $\bm{M}$.

Part-II algorithm

We take $\bm{M}_{\textrm{best}}$ from the Part-I algorithm as an initial design in the Part-II algorithm. We generate a neighbor of FSLHD based on the different-slice or the out-slice exchange procedures in the Part-II algorithm. For $i=1,\cdots,u$, if $q$ is large and $n_{i}^{q}>>n$, then the $n$ design points is very sparse by the Part-I algorithm, consequently, the Part-II algorithm brings smaller effect for the space-filing properties of the design $\bm{D}$. Therefore, in this case, the Part-I algorithm is more important, and we can skip the Part-II algorithm and focus on the Part-I algorithm. We also can stop the running of Part-II algorithm when the repeating times arrive 100.

1. Step 1.

Let $\bm{M}=\bm{M}_{\textrm{best}}$, and set the index $i=1$. 2. Step 2.

In the $i$th slice of $\bm{M}$, generate an $\bm{M}_{\textrm{N}}$ by the different-slice or the out-slice exchange procedures under the condition of $P(\lceil\bm{M}_{\textrm{N}}/t_{i}\rceil)=0$. 3. Step 3.

If $\phi_{\textrm{CSM}}(\bm{D}_{\textrm{N}})<\phi_{\textrm{CSM}}(\bm{D})$, replace $\bm{M}$ by $\bm{M}_{\textrm{N}}$. 4. Step 4.

Repeat Step 2 and Step 3 until meeting the stopping criterion. 5. Step 5.

Update $i=i+1$, if $i<u$, go to Step 2; else, output $\bm{M}_{\textrm{best}}$ = $\bm{M}$.

4 Simulation results

In this section, the first example illustrates that the SESE algorithm has good properties. In our second example, for the design with large runs and factors, we give some comparative studies, which show the efficient two-part algorithm with desirable performance. In these examples, we select the combined space-filling measurement (6). For simplicity, we only consider any column of FSLHDs with all in (1) being 1/2 when updating (6) in our proposed algorithm.

4.1 Example 1

As depicted in Figure 5(a), we randomly generate an initial design FSLHD($4,8,12$; 3,2) with optimal univariate uniformity. It is clear that the space-filling property is poor for the whole design and for each slice of the design. Based on the combined space-filling measurement $\phi_{\textrm{CSM}}$ ($t=50$ ) in (6), we improve the space-filling property of the design by the SESE algorithm ($P=20$). The initial design with $\phi_{\textrm{CSM}}=14.4740$ is showed in Figure 5(a). After operating the SESE algorithm, the resulting design with $\phi_{\textrm{CSM}}=5.7958$ in Figure 5(b) has good space-filling property over the experiment region.

For comparison, we randomly generate FSLHDs by the method in Section 2 for 100000 times and calculate the corresponding values of $\phi_{\textrm{CSM}}$. The resulting FSLHDs with good space-filling properties account for a small portion of 100000 FSLHDs. The smallest value of $\phi_{\textrm{CSM}}$ from the 100000 FSLHDs is 6.8387, while the value of $\phi_{\textrm{CSM}}$ in Figure 5(b) is 5.7958. The values between 6.8387 and 8 of $\phi_{\textrm{CSM}}$ account for 0.22 percent of all $\phi_{\textrm{CSM}}$ values from the 100000 FSLHDs. It can be seen that the SESE algorithm is useful to improve the space-filling property of the whole design and each slice of the initial design.

4.2 Example 2

To show the good performance of the two-part algorithm for design with large runs and factors, we compare its performance with the SESE algorithm. We repeat each algorithm for 100 times with a random initial design FSLHD in Table 1 (SD, standard deviation and AT, average time ). In the SESE algorithm, we set stopping rules $P=30$ for FSLHD$(15,30;2,2)$ and $P=40$ for FSLHD$(5,10,15,30;4,6)$. Conclusion can be obtained from Table 1 as follows:

1. (i)

The average time of the operation shows that the two-part algorithm has higher efficiency than the SESE algorithm. 2. (ii)

For FSLHD$(5,10,15,30;4,6)$, since $n_{i}^{q}>>n$ with $i=1\cdots 4$, the $\phi_{\textrm{CSM}}$ values of the resulting FSLHD from Part-I algorithm are desirable when compared with those values from two-part algorithm. However, the results of Part-I algorithm for FSLHD$(15,30;2,2)$ are not good enough. Therefore, if $q$ is large and $n_{i}^{q}>>n$, we need not to run the Part-II algorithm. 3. (iii)

Based on the $\phi_{\textrm{CSM}}$ values of the resulting FSLHDs, we can see that the $\phi_{\textrm{CSM}}$ values are close to each other. It can be concluded that both the two-part algorithm and the SESE algorithm are stable and do not heavily rely on the initial design.

By comparison, the resulting designs is better after using SESE algorithm. However, for generating space-filing FSLHDs with large runs and factors as well as considering the cost of time, the two-part algorithm is preferable.

5 Discussion the methods for evaluating the combined space-filling measurement

Recall that $\bm{D}$ is the design matrix of an ${\rm FSLHD}(n_{1},\ldots,n_{u};u,q)$. Since we generate a neighbor of design by exchanging two elements in one column of $\bm{D}$, we do not need to recalculate all the inter-site distances when we update $\phi_{t}(\bm{D})$ or $\phi_{t}(\bm{D}^{(i)})$. The calculative efficiency of optimality criteria for the LHD has been discussed in Jin et al. (2016). Here, based on above three exchange procedures for the FSLHD, we give updating expressions of $\phi_{\textrm{CSM}}(\bm{D})$ using the previous $\phi_{t}(\bm{D})$ and $\phi_{t}(\bm{D}^{(i)})$ for our proposed algorithm.

For the design matrix $\bm{D}=(x_{ij})_{n\times q}$ with $n$ design points $\{\bm{x}_{1},\cdots,\bm{x}_{n}\}$, we exchange $x_{rk}$ and $x_{sk}$ in the $k$th column of the design. Let $d(\cdot,\cdot)$ be the inter-site distance before exchanging. Let $v\neq r,s$, $1\leq v\leq n$, as defined in (3), the new related inter-site distance of the two design points $\bm{x}_{r}$ and $\bm{x}_{s}$ should be updated:

[TABLE]

[TABLE]

where $h(r,s,k,v)=|x_{sk}-x_{vk}|^{m}-|x_{rk}-x_{vk}|^{m}$ and the other inter-site distances are unchanged. We give a new $\phi_{\textrm{CSM}}^{\prime}(\bm{D})$ based on the previous $\phi_{t}(\bm{D})$ and $\phi_{t}(\bm{D}^{(i)})$ as follows:

[TABLE]

For three different procedures, the values of $\phi_{t}^{\prime}(\bm{D})$ and $\phi_{t}^{\prime}(\bm{D}^{(i)})$ are determined as follow:

(i) Within-slice exchange procedure. For $e\in\{1,\cdots,u\}$, suppose that the design points $\bm{x}_{r}$ and $\bm{x}_{s}$ are in the $e$th slice. Let $n_{0}=0$. Then we have $r,s\in J_{e}=\{\sum_{l=0}^{{e-1}}n_{l}+1,\cdots,\sum_{l=0}^{{e}}n_{e}\}$ and

[TABLE]

[TABLE]

(ii) Different-slice exchange procedure. For $e,e^{\prime}\in\{1,\cdots,u\}$, suppose that the design points $\bm{x}_{r}$ are in the $e$th slice and $\bm{x}_{s}$ in the $e^{\prime}$th slice. Let $n_{0}=0$. Then we have $r\in J_{e}=\{\sum_{l=0}^{{e-1}}n_{l}+1,\cdots,\sum_{l=0}^{{e}}n_{e}\}$, $s\in J_{e^{\prime}}=\{\sum_{l=0}^{{e^{\prime}-1}}n_{l}+1,\cdots,\sum_{l=0}^{{e^{\prime}}}n_{e^{\prime}}\}$ and

[TABLE]

[TABLE]

(iii) Out-slice exchange procedure. For $e\in\{1,\cdots,u\}$, suppose that the element $\bm{x}_{rk}$ are in the $e$th slice and the element $x_{sk}^{\prime}\in(0,1)$ are in the out slice. Let $n_{0}=0$. Then we have $r\in J_{e}=\{\sum_{l=0}^{{e-1}}n_{l}+1,\cdots,\sum_{l=0}^{{e}}n_{e}\}$, $h(r,s,k,v)=|x_{sk}^{\prime}-x_{vk}|^{m}-|x_{rk}-x_{vk}|^{m}$, $v\neq r$ and

[TABLE]

[TABLE]

Through the above description of the updating formulas, we can improve the efficiency of re-evaluating $\phi_{\textrm{CSM}}(\bm{D})$ for our proposed algorithm.

6 Conclusions

In this article, we propose a method to construct SLHDs with arbitrary run sizes. Based on such designs, we give an SESE algorithm to search the optimal FSLHDs. Moreover, we provide an efficient two-part algorithm to improve the optimization efficiency in generating the space-filling FSLHDs with large runs and factors. We believe that FSLHDs with optimal univariate uniformity and good space-filling properties are more widely used in computer experiments. Orthogonality is also an appealing feature for SLHDs. Orthogonal SLHDs are constructed in Yang et al. (2013); Huang et al. (2014); Cao and Liu (2015), however, orthogonality does not ensure a good space-filling property. In the future, we will study the construction of an orthogonal-maximin SLHD with slices of arbitrary run sizes. Such a design have both orthogonality and space-filling property.