Computer experiments are widely used in scientific research and industrial applications, where complex computer codes, often high-fidelity simulators, are used to generate data instead of physical experimentation [9, 35]. Because their outputs are deterministic and free of random errors, replication is unnecessary [2, 3, 22]. Latin hypercube designs (LHDs, [30]), characterized by their avoidance of replication in each dimension and their uniform one-dimensional projections, may be regarded the most popular design class for computer experiments [9, 48]. Depending on practical objectives, different forms of optimal LHDs are used, including space-filling, maximum projection, and orthogonal designs. Although numerous construction methods have been developed, obtaining high-quality LHDs for moderate to large design sizes remains challenging [9, 21, 48, 50]. One challenge is that theoretical results, such as algebraic constructions guaranteeing properties like maximin distance or orthogonality, only exist for certain design sizes, which limits their general applicability. For sizes where such guarantees are unavailable, search algorithms are typically employed; however, their performance varies with the choice of algorithm, the extent of the search space explored, and the computational resources allocated, meaning global optimality cannot be ensured. More fundamentally, optimal LHD construction is a discrete optimization problem, and although exhaustive enumeration would identify the true optimum, it becomes infeasible for all design sizes due to the exponential growth of the permutation space.

An LHD with n runs and k factors is represented by an \(n \times k\) matrix in which each column is a permutation of integers \(1, \ldots , n\). Throughout this paper, n refers to the run size and k refers to the factor size. Space-filling LHDs aim to distribute design points as uniformly as possible across the input domain, thereby reducing unsampled regions in all dimensions. Various criteria have been proposed to quantify this property, including maximin and minimax distance measures [19, 31], different forms of discrepancy [9, 10, 16], and entropy-based measures [36]. Because the number of possible LHDs increases exponentially as \((n!)^{k}\) for a given n and k, an exhaustive enumeration is impractical for all design sizes. The current literature contains search algorithms [2, 6, 14, 18, 22, 25, 26, 31, 51] and algebraic constructions [44, 47, 52] designed to produce space-filling LHDs.

Search algorithms are commonly used to construct space-filling LHDs with flexible run sizes. For example, [31] proposed a simulated annealing (SA) algorithm that helps avoid being trapped in local optima and improves the likelihood of finding globally optimal designs. Building on their work and that of [25, 41] developed an SA-based approach for searching space-filling orthogonal array–based LHDs (OALHDs). [22] introduced a multi-objective criterion and an adapted SA algorithm that simultaneously accounts for orthogonality and space-filling properties. [2] extended the framework of sliced Latin hypercube designs (SLHDs; [33]) and proposed a two-stage SA algorithm. Beyond SA-based algorithms, several other search algorithms have been developed to identify optimal designs. [51] proposed the columnwise pairwise (CP) algorithm for generating efficient symmetric LHDs. [18] introduced the enhanced stochastic evolutionary (ESE) algorithm, which combines an exchange procedure (inner loop) with a threshold-determination mechanism (outer loop). [26] employed a genetic algorithm (GA) with a strategy that directly updates toward the global best. [14] explored an iterated local search heuristic that alternates between local improvement and perturbation. [6] developed a particle swarm optimization (PSO) algorithm in which the search process gradually reduces the Hamming distance between each particle and the global or personal best through element exchanges. In general, search algorithms can be computationally intensive, especially for constructing large designs. For certain design sizes, algebraic constructions provide an attractive alternative, often yielding optimal LHDs at very low computational cost. For instance, [47] used Costas arrays to construct space-filling saturated LHDs or Latin squares with run sizes \(n=p-2\), \(p-1\) or p, where p is any prime or prime power. [44] applied the Williams transformation [45] to linearly permuted good lattice point sets to obtain maximin LHDs when the number of factors does not exceed the count of positive integers that are coprime to the run size.

While traditional space-filling designs emphasize coverage of the full dimensional input space, it is often equally important to ensure space-filling properties in intermediate dimensions. To address this need, [21] introduced maximum projection (MaxPro) designs, which enhance space-filling performance across all sub-spaces of dimension 2 through \(k-1\). Compared with classical maximin distance LHDs, MaxPro LHDs typically yield better separation in these lower-dimensional views. However, generating high-quality MaxPro LHDs remains challenging, particularly for large design sizes. [21] proposed an SA-based algorithm for their construction. In this work, we explore a GA framework that is capable of producing MaxPro LHDs with substantially improved performance, and illustrative examples are provided in Section 5.

Unlike space-filling LHDs, which aim to reduce similarity among design points (rows), orthogonal LHDs (OLHDs) focus on minimizing dependence among factors (columns). OLHDs are particularly appealing because they achieve zero column-wise correlation, thereby facilitating more reliable estimation of factor effects. Algebraic constructions have been developed for specific design sizes. For example, [50] proposed an algebraic construction of OLHDs with run sizes \(n=2^m+1\) and factor sizes \(k=2m-2\), where m is any integer no less than 2. [7] further extended this approach to accommodate more factors. [37] developed a construction method based on factorial designs and group rotations for designs with \(n=2^{2^m}\) and \(k=2^{m}t\), where m is any positive integer and t is the number of rotation groups. [39] extended their earlier work [38] to construct OLHDs with \(n=r2^{c+1}\) or \(n=r2^{c+1}+1\) and \(k=2^c\), where c and r are positive integers. [49] proposed using generalized orthogonal designs to construct OLHDs and nearly orthogonal LHDs (NOLHDs) with \(n=2^{r+1}\) or \(n=2^{r+1}+1\) and \(k=2^r\), where r is any positive integer. [12] introduced a construction method based on orthogonal arrays and their full fold-overs, yielding designs with \(n=2ak\) runs and k factors, where k is the dimension of the orthogonal matrix and a is any positive integer. [3] applied the Williams transformation [45] to construct OLHDs under a second-order cosine model for odd prime n and \(k\le n -1\). In addition, several new types of optimal LHDs can be derived by leveraging properties of orthogonal arrays. For example, [41] proposed to construct OALHDs by deterministically replacing elements in orthogonal arrays, resulting in designs with improved space-filling properties and low column-wise correlations. [27] proposed using OLHDs or NOLHDs as starting designs and coupling them with orthogonal arrays to obtain OALHDs. Their approach requires fewer runs to support large numbers of factors, producing designs with \(n^2\) runs and 2fp factors, where n and p are design sizes of the OLHDs or NOLHDs, and 2f is the number of columns in the coupled orthogonal array. However, many design sizes fall outside the scope of these theoretical results, and in such cases OLHDs can be obtained through search algorithms. The optimization algorithms discussed above can be adapted to the target orthogonality.

In this paper, we provide a comprehensive review and comparison of current popular algebraic constructions for optimal LHDs and widely adopted search algorithms for generating high-quality LHDs. By examining their theoretical properties and empirical performance across a variety of criteria, we present a unified perspective on their respective strengths and limitations. The comparisons offered here are intended to guide practitioners in selecting design strategies suited to different objectives and constraints. In conducting these evaluations, we also identify several designs that outperform existing ones; these improved designs are reported in Section 5.

The remainder of the paper is organized as follows. Section 2 introduces different optimality criteria for LHDs. Section 3 describes some widely used search algorithms and their implementations. Section 4 presents several popular algebraic constructions. Section 5 reports numerical results and comparisons. Section 6 concludes with a discussion and potential directions for future research.

A variety of criteria have been proposed to quantify designs’ space-filling performances [9, 10, 16, 19, 36]. Among these, the maximin distance criterion [19] has become one of the most widely used. It aims to distribute design points throughout the experimental region so that the smallest pairwise distance between points is as large as possible, thereby promoting uniform coverage. Let \({\textbf {X}}\) denote an LHD matrix throughout this paper. For two runs \(x_i\) and \(x_j\) of \({\textbf {X}}\), the \(L_q\)-distance is defined as \(d_q(x_i, x_j) = \left\{ \sum _{m=1}^{k} \vert x_{im}-x_{jm}\vert ^q \right\} ^{1/q}\), where q is an integer, \(1 \le i \le n\), \(1 \le j \le n\), and \(i \ne j\). In this paper, we focus on the cases \(q=1\) and \(q=2\), corresponding to the Manhattan (\(L_1\)) and Euclidean (\(L_2\)) distances. The overall \(L_q\)-distance of the design \({\textbf {X}}\) is defined as \(d_q({\textbf {X}}) = \text {min} \{d_q(x_i, x_j), 1 \le i<j \le n \}\). A design \({\textbf {X}}\) is said to be a maximin \(L_q\)-distance design if it achieves the largest \(d_q({\textbf {X}})\) among all LHDs of the same size. If multiple designs share the same largest \(d_q({\textbf {X}})\), the maximin distance design sequentially maximizes the next minimum inter-site distances. To facilitate numerical optimization of the maximin distance criterion, Morris and Mitchell [31] and Jin et al. [18] proposed to minimize a scalar value:

where p is a tuning parameter. This \(\phi _{p}\) criterion in Equation (1) approaches the maximin distance criterion as \(p \rightarrow \infty \). In practical applications, choosing \(p=15\) typically suffices [31].

Maximin distance LHDs primarily emphasize the space-filling properties in the full k-dimensional design space, but they do not necessarily ensure good separation when designs are projected onto lower-dimensional subspaces. To address this limitation, [21] introduced the maximum projection (MaxPro) criterion, which explicitly evaluates space-filling performance in all possible projections. An LHD \({\textbf {X}}\) is a MaxPro LHD if it minimizes

Equation (2) penalizes pairs of design points that are close in any coordinate, thereby encouraging separation in every marginal projection. Therefore, MaxPro LHDs promote uniformity not only in the full-dimensional space but also across all lower-dimensional subspaces. It should be noted that the original MaxPro criterion was defined for designs scaled to the unit hypercube \([0,1]^{k}\) [21], whereas the designs considered here use integer-valued levels. The two representations are equivalent up to a simple linear transformation. For instance, scaling can be reversed by applying the transformation \({\textbf {X}}_{Integer} = {\textbf {X}}_{Scaled}*n-0.5\), which maps unit-hypercube points to integer-level designs. This transformation is implemented in Section 5 to ensure a fair comparison of the design performance.

Orthogonal and nearly orthogonal designs, which seek to reduce the dependence among factors, are widely employed in experiments [11, 37, 40]. The degree of orthogonality in a design is commonly quantified using correlation-based measures. Two frequently used criteria are the average absolute correlation and the maximum absolute correlation [11], denoted by ave(|q|) and max|q|, respectively:

where \(q_{ij}\) represents the correlation between the ith and jth columns of the design matrix \({\textbf {X}}\). A design is considered orthogonal if both ave(|q|) and max|q| equal zero, a condition that cannot be achieved for many combinations of run sizes and factor dimensions. In practical applications, designs with small values of either criterion are generally favored.

Simulated annealing (SA, [24]) is a stochastic optimization technique inspired by the physical annealing process in metallurgy, in which a material is heated to a high temperature and then slowly cooled to reduce structural defects. As a single solution based algorithm, SA begins from an initial candidate solution and iteratively explores the search space, accepting or rejecting new solutions according to a probabilistic rule until convergence is reached. The algorithm has been widely applied to optimization problems involving continuous variables. To accommodate the discrete nature of LHDs, the classical SA framework has been adapted by researchers to enable effective search over the permutation-based design space. In this article, we focus on minimizing the optimality criteria outlined in Section 2, meaning only minimization optimization problems are considered.

Morris and Mitchell [31] proposed a modified simulated annealing algorithm in which elements of an LHD are randomly exchanged to search for potential improvements. Specifically, the algorithm begins with a randomly generated LHD, from which a column is selected at random. Then, two random elements within this column are swapped to produce a new candidate design. For example, consider an LHD with four runs in which a selected column contains the elements (1, 2, 3, 4). Exchanging the first two elements yields the column (2, 1, 3, 4), resulting in a new LHD. If this exchange improves the design under a specified optimality criterion, it is accepted; otherwise, it is accepted with probability \(\hbox {exp}[-(\Phi ({\textbf {X}}_{new})-\Phi ({\textbf {X}}))/T]\), where \(\Phi \) denotes the optimality criterion, \({\textbf {X}}\) is the current design, \({\textbf {X}}_{new}\) is the design obtained after the exchange, and T is the current temperature. This exchange procedure is applied iteratively to progressively improve the design. When no improvement is observed after a prescribed number of attempts, the temperature T is reduced according to the annealing schedule, thus decreasing the probability of accepting inferior designs and promoting convergence. The algorithm ends when convergence is achieved or a predefined computational budget is exhausted, at which point the best LHD encountered is returned. The pseudo-code of this SA is outlined in Algorithm 1.

Leary et al. [25] extended the simulated annealing algorithm of Morris and Mitchell [31] to construct orthogonal array–based Latin hypercube designs (OALHDs). Tang [41] demonstrated that OALHDs generally exhibit superior space-filling properties compared to randomly generated LHDs. The algorithm proposed by Leary et al. [25] begins with a randomly selected OALHD, from which a column is chosen at random. Then, two elements within this column that correspond to the same entry in the underlying orthogonal array (OA) are exchanged. For instance, consider an OALHD with nine runs in which a selected column contains elements (1, 2, 3, 4, 5, 6, 7, 8, 9), and the corresponding OA has entries \(({ 1,1,1},{ 2,2,2},{ 3,3,3})\). In this case, elements (1, 2, 3) share the same original OA entry \({ 1}\), elements (4, 5, 6) share the entry \({ 2}\), and elements (7, 8, 9) share the entry \({ 3}\). The remaining steps of the algorithm follow those of the simulated annealing procedure described in [31]. It is important to note that the existence of OALHDs depends on the availability of the corresponding underlying orthogonal arrays.

Joseph and Hung [22] proposed another variant of simulated annealing for constructing orthogonal–maximin LHDs, which simultaneously account for orthogonality and maximin distance criteria. The algorithm begins by generating a random LHD and then identifies the column exhibiting the largest average pairwise correlation with all other columns. Specifically, the \(l^{*}\)th column is selected, where \(l^{*} = \mathop {\textrm{argmax}}\limits \rho _{l}^2\) with \(\rho _{l}^2=\frac{1}{k-1}\sum _{j\ne l}\rho _{lj}^2\), and \(\rho _{lj}\) denotes the correlation between the lth and jth columns. Next, the algorithm selects the row with the largest total row-wise distance from all other rows. That is, the \(i^{*}\)th row is chosen such that \(i^{*}=\mathop {\textrm{argmax}}\limits (\sum _{j\ne i}d(x_i, x_j)^{-p})^{1/p}\), where \(d(x_i, x_j)\) represents the distance between the ith and jth rows. The element located at the intersection of the \(i^{*}\)th row and the selected column is then swapped with a randomly chosen element of the same column. The remaining steps of the procedure follow the simulated annealing framework described in [31]. Although this algorithm is capable of producing LHDs that simultaneously minimize the \(\phi _{p}\) and ave(|q|) criteria, it is computationally intensive, as all \(\rho _{l}^2\) values and row-wise distances must be recomputed at each iteration. Table 1 summarizes the similarities and differences among the three simulated annealing algorithms discussed above.

[2] and [33] proposed methods for constructing space-filling sliced Latin hypercube designs (SLHDs). An SLHD with n runs can be partitioned into t slices, \({\textbf {X}}=\bigcup _{i=1}^{t} {\textbf {X}}_{i}\), where each slice \({\textbf {X}}_{i}\) is itself an LHD with \(m=n/t\) runs. They introduced a simulated annealing algorithm consisting of two stages. The first stage optimizes the design from the slice perspective and begins with t randomly generated slices, each being an \(m \times k\) LHD. If duplicated rows are detected across slices, the algorithm randomly selects one duplicated row, then chooses another row within the same slice, and swaps the two elements in these rows for a randomly selected column. This procedure is repeated until no duplicate rows remain. Subsequently, elements in a randomly chosen column of a randomly selected slice are exchanged. If this exchange results in duplicated rows, the previous procedure is invoked again; otherwise, the exchange is retained when it leads to an improvement. The second stage refines the design at the element level by replacing the t occurrences of the level l (\(l=1, \ldots , m\)) with a random permutation of \(\{ (l-1)t+1,\ldots , lt \}\), followed by swapping two randomly selected elements within that permutation. When \(t=1\), the algorithm reduces to a procedure for constructing classical space-filling LHDs.

Particle swarm optimization (PSO, [23]) is a metaheuristic optimization algorithm inspired by collective behaviors observed in nature, such as coordinated flight patterns of bird flocks. PSO is a population-based algorithm in which a set of particles, representing candidate solutions, explores the search space by updating their positions and velocities according to predefined mathematical rules, with the goal of iteratively improving a specified optimality criterion. The movement of each particle is influenced by both its own best-known position (personal best) and the best positions discovered by the swarm as a whole (global best). Through this cooperative information-sharing mechanism, the swarm is expected to converge toward high-quality solutions.

Recent studies [5, 6, 46] have adapted the classical PSO framework to the construction of high-quality experimental designs, a problem characterized by a discrete search space. In particular, Chen et al. [6] proposed a PSO-based algorithm, known as LaPSO, to identify maximin distance LHDs. This approach modifies the standard PSO update rules for particle velocities and positions. Specifically, it seeks to reduce the Hamming distance between each particle and the global best (or its corresponding personal best) by deterministically exchanging elements. Here, the Hamming distance between two LHD particles is defined as the number of positions at which their corresponding elements differ.

LaPSO begins by generating m random LHDs, where m denotes the number of particles and serves as a tuning parameter. Initially, the personal best (PB) of each particle is set to the particle itself, while the global best (GB) is defined as the particle LHD that achieves the best value under a specified optimality criterion \(\Phi \). The algorithm then applies an exchange procedure to each column of every particle LHD. Specifically, for a given column, an element in the current particle (denoted as \(e_r\))is randomly selected, and the corresponding element located at the same position in the PB (denoted as \(e_p\)) is identified. The elements \(e_r\) and \(e_p\) are then swapped within the current column, thereby reducing the Hamming distance between the current particle and its PB. An analogous exchange procedure is performed to reduce the Hamming distance between the current particle and the GB. These two operations emulate particle movement toward the PB and the GB, respectively.

Either of these steps may be skipped or repeated, but they cannot be omitted simultaneously. Although no strict upper bound is imposed on the number of repetitions, performing a large number of exchanges, i.e. 2n, will make the current column nearly identical to that of the PB or the GB. After all particle LHDs have been updated, a random number z is drawn from the standard uniform distribution for each column of each particle. If \(z < p_0\), where \(p_0\) is a user-defined tuning parameter, two randomly selected elements in the current column are exchanged. This random perturbation helps prevent premature convergence to local optima. Once the algorithm converges or the computational budget is exhausted, the GB is returned as the final solution. Unlike the exchange mechanisms used in SA based algorithms, LaPSO explicitly links each particle to its personal best, the global best, or both. However, since each exchange operation can reduce the Hamming distance by at most one, LaPSO may become computationally demanding when constructing large high-quality designs. The pseudo-code of LaPSO is outlined in Algorithm 2.

Genetic algorithms (GAs) are nature-inspired metaheuristic optimization algorithms that emulate the principle of natural selection proposed by Charles Darwin [13, 17]. A GA is a population-based algorithm that typically consists of stages of \({ selection}\), \({ crossover}\), \({ mutation}\) and \({ fitness}\) \({ evaluation}\). In accordance with evolutionary terminology, candidate solutions are referred to as chromosomes, which together constitute a population. Based on their performance under a specified criterion, called the fitness, a subset of chromosomes is selected as parents for generating a new population, a process known as selection. The selected parents are then recombined through crossover and mutation operators to produce offspring that are expected to exhibit improved performance.

[26] proposed a GA for constructing maximin distance LHDs. The initial population consists of m randomly generated LHDs, where m is an even, user-defined population size. During the \({ selection}\) stage, the entire population is evaluated according to the specified optimality criterion, and the top half of the population is retained as \({ survivors}\). Among these, the current global best is designated as \({ best}\) \({ survivor}\) (BS). The BS and the remaining survivors are then used to generate a new population.

Specifically, the BS is assigned as the first individual, and the second through the (m/2)th individuals are created via the following \({ crossover}\) procedure. For each survivor other than the BS, a column is randomly selected from that survivor and used to replace the column with the same index in the BS, thereby generating a new individual. For the remaining half of the population, the BS is designated as the \((m/2+1)\)th individual, and the \((m/2+2)\)th through the mth individuals are generated through an alternative \({ crossover}\) operation. In this case, for each survivor except the BS, a column is randomly selected from the survivor and replaced by the column with the same index in the BS to form a new individual.

After constructing the new population of m LHDs, a \({ mutation}\) step is applied to reduce the risk of premature convergence. During \({ mutation}\), for each column of every new LHD except the first design, a random number z is drawn from the standard uniform distribution. If \(z < p_{mut}\), where \(p_{mut}\) is a user-defined tuning parameter, two randomly selected elements within that column are exchanged. Next, the entire population is then re-evaluated to identify the survivors for the next iteration. Once the algorithm converges or the computational budget is exhausted, the BS is returned as the final design.

Unlike LaPSO, which relies on element-wise exchanges, this GA operates by replacing entire columns in both the current global best design and candidate designs. Consequently, the Hamming distance for the replaced columns becomes zero immediately with respect to the global best. Compared with other search algorithms, this approach typically requires less CPU time, particularly for large design sizes. The pseudo-code of GA is outlined in Algorithm 3.

For certain design sizes, algebraic constructions are available, and theoretical results have been developed to guarantee the efficiency of the resulting designs. Such constructions require little to no search and are therefore especially attractive for designs with large sizes. In this section, we review several practically useful algebraic constructions for maximin distance LHDs and orthogonal LHDs.

[44] proposed to generate maximin distance LHDs via good lattice point (GLP) sets [52] and Williams transformation [45]. In practice, their method can lead to space-filling designs with relatively flexible sizes, where the run size n is flexible but the factor size k must be no greater than the number of positive integers that are co-prime to n. They proved that the resulting designs of sizes \(n \times (n-1)\) (with n being any odd prime) and \(n \times n\) (with \(2n+1\) or \(n+1\) being odd prime) are optimal under the maximin \(L_1\) distance criterion.

The construction method in [44] can be summarized into three steps. First, generate an \(n \times k\) GLP design D whose element is \(x_{ij}=i \times h_{j}\) (mod n) with \(i=1, \ldots , n\), \(j=1, \ldots , k\) and \(h=(h_1, \ldots ,h_k)\) being a set of distinct positive integers that are coprime to n. Second, for any \(b\in \left\{ 0,\ldots ,n-1\right\} \), generate \(D_b=D+b \text { (mod } n\text {)}\) and \(E_b=W(D_b)\), where \(W: \mathcal {Z}_n \rightarrow \mathcal {Z}_n\) is the Williams transformation [45]:

Third, find the best \(b^{*} \in \{ 0,\ldots ,n-1\}\) such that \(E_{b^{*}}\) has the smallest \(\phi _p\) value.

As shown in [52] and [44], the designs D and E defined above are \(n \times k\) LHDs whose last rows contain the same elements. By removing the last rows and reordering the levels, the remaining designs can be transformed into LHDs, referred to as leave-one-out designs \(\tilde{D}\) and \(\tilde{E}\). After performing the three steps outlined above, the resulting designs \(E_{b^{*}}\) and \(\tilde{E}_{b^{*}}\) exhibit good space-filling properties. The construction method proposed by [44] is particularly attractive for generating large maximin distance LHDs; however, it may be too complex for practitioners to implement directly. In the LHD package in R, the function FastMmLHD() implements this method with a user-friendly interface.

[41] proposed to construct orthogonal array-based LHDs (OALHDs) from existing orthogonal arrays (OAs). Suppose \({\textbf {A}}\)=OA(n, m, s, r) is an orthogonal array with n rows, m columns, s levels (\(n > s \ge 2\)) and r strength. If each \(n \times r\) sub-matrix of \({\textbf {A}}\) includes all possible \(1\times r\) row vectors with the same frequency \(\lambda \), then \(\lambda \) is called the index of the array where the run size \(n=\lambda s^r\). An \(n \times m\) LHD can be considered as an OA(n, m, n, 1) with \(\lambda = 1\). For every column of an OA, [41] proposed to replace the \(ns^{-1}\) (equal to \(\lambda s^{r-1}\)) positions of entry k by a permutation of numbers: \((k-1)ns^{-1}+1, (k-1)ns^{-1}+2, \ldots , (k-1)ns^{-1}+ns^{-1}=kns^{-1}\), where \(k=1,\ldots ,s\), then the new generated design matrix would be an \(n \times m\) LHD. [41] showed that OALHDs can have better space-filling properties than general ones.

Orthogonal LHDs (OLHDs) have zero pairwise correlations between any two columns and are therefore widely used in practice. There is a rich literature on the construction of OLHDs for various design sizes. In this section, we summarize several currently popular construction methods [3, 7, 27, 39, 41, 50]. The corresponding design sizes are summarized in Table 2.

[50] proposed a construction for OLHDs with run sizes \(n=2^m+1\) and factor sizes \(k=2m-2\) where m is any integer greater than 2. It involves the constructions of three matrices \({\textbf {M}}\), \({\textbf {S}}\), and \({\textbf {T}}\). The first column in \({\textbf {M}}\), denoted as \({\textbf {e}}\), is a random permutation of \((1, 2,\ldots , 2^{m-1})\). The 2nd to the mth columns in \({\textbf {M}}\) are calculated as \({\textbf {A}}_{L} {\textbf {e}}\), where \(L=1,2,\ldots ,m-1\). Each \({\textbf {A}}_{L}\) is defined as:

where \({\textbf {I}}\) is the \(2 \times 2\) identity matrix, \({\textbf {R}} = {\textbf {1}} {\textbf {1}}^T -\) \({\textbf {I}}\) and \(\otimes \) is the Kronecker product. The last \(m-2\) columns in \({\textbf {M}}\) are calculated as \({\textbf {A}}_{i} {\textbf {A}}_{m-1} {\textbf {e}}\), where \(i=1,2,\ldots ,m-2\). Similarly, to construct the matrix \({\textbf {S}}\), set its first column, denoted as \({\textbf {j}}\), to be the \(2^{m-1}\times 1\) vector of \(+1\)’s. The 2nd to the mth columns in \({\textbf {S}}\), denoted as \({\textbf {a}}_{K}\) with \(K=1,2,\ldots ,m-1\) are defined as \({\textbf {a}}_{K}={\textbf {B}}_{1} \otimes {\textbf {B}}_{2} \otimes \ldots \otimes {\textbf {B}}_{m-1}\), where all the \({\textbf {B}}\)’s are \([1,1]^T\) except for \({\textbf {B}}_{m-k}=[-1,1]^T\). The last \(m-2\) columns in \({\textbf {S}}\) are calculated as \({\textbf {a}}_{1} {\textbf {a}}_{j}\), where \(j=2,\ldots ,m-1\). Matrix \({\textbf {T}}\) is the element-wise product of \({\textbf {M}}\) and \({\textbf {S}}\), and the matrix \({\textbf {X}}\) is constructed by \({\textbf {X}}=[{\textbf {T}}^{T},{\textbf {0}},{\textbf {-T}}^{T}]^T\), where \({\textbf {0}}\) is a column vector containing k zeros and \({\textbf {-T}} = -1 \times {\textbf {T}}\). For any choice of \({\textbf {e}}\) in \({\textbf {M}}\), the resulting design \({\textbf {X}}\) is an OLHD.

[7] extended [50]’s method to construct OLHDs with run size \(n=2^m+1\) and factor size \(k=m+ {m-1 \atopwithdelims ()2}\), where m is any integer greater than 2. Their construction adopted the same matrix \({\textbf {T}}\) in [50] but used different matrices \({\textbf {M}}\) and \({\textbf {S}}\). In \({\textbf {M}}\), for the last \(m-2\) columns, instead of using \({\textbf {A}}_{i} {\textbf {A}}_{m-1} {\textbf {e}}\) in [50], they adopted \({\textbf {A}}_{i} {\textbf {A}}_{j} {\textbf {e}}\) with \(i=1,\ldots ,m-2\) and \(j=i+1,\ldots ,m-1\), which accommodates more factors. Similarly, for the matrix \({\textbf {S}}\), for the last \(m-2\) columns, they adopted \({\textbf {a}}_{i} {\textbf {a}}_{j}\) with \(i=1,\ldots ,m-2\) and \(j=i+1,\ldots ,m-1\). Given the same number of runs, the method in [7] is capable of accommodating more factors. Note that the choice of \({\textbf {e}}\) in constructing \({\textbf {M}}\) is important. If \({\textbf {e}}=[1, 2,\ldots , 2^{m-1}]^T\), then \({\textbf {X}}\) is guaranteed to be orthogonal, while a poor choice of \({\textbf {e}}\) may not guarantee \({\textbf {X}}\) to be orthogonal.

[39] extended their earlier work [38] to construct OLHDs with \(n=r2^{c+1}+1\) or \(n=r2^{c+1}\) and \(k=2^c\), where r and c are positive integers. Their method for constructing OLHD with \(n=r2^{c+1}+1\) and \(k=2^c\) consists of constructing three matrices: \(S_c\), \(T_c\) and \(A_{r2^c\times 2^c}\). The matrix \(S_c\) is defined as

\(S_c=\begin{bmatrix} 1 & 1 \\ 1 & -1 \end{bmatrix}\) for \(c=1\), \(S_c=\begin{bmatrix} S_{c-1} & -S_{c-1}^{*} \\ S_{c-1} & S_{c-1}^{*} \end{bmatrix}\) for \(c>1\),

where the operator \(*\) works on any matrix with an even number of rows by multiplying the top half entries by \(-1\). The matrix \(T_c\) is defined as

The matrix \(A_{r2^c\times 2^c}\) is defined as \(A_{r2^c\times 2^c}=[(T_c^{1})^{T},\ldots ,(T_c^{r})^{T}]^{T}\), where \(T_c^{i}=T_c+(i-1)2^cS_c\) for \(i=1,\ldots ,r\). The design matrix \({\textbf {X}}\) is given by \({\textbf {X}}=[A_{r2^c\times 2^c}^{T},{\textbf {0}},-A_{r2^c\times 2^c}^{T}]^T\), where \({\textbf {0}}\) is a column vector containing \(2^c\) zeros. The construction method for OLHD with \(n=r2^{c+1}\) and \(k=2^c\) also relies on the \(S_c\) and \(T_c\). Let \(H_c=T_c-S_c/2\), \(B_{r2^c\times 2^c}\) is defined as \(B_{r2^c\times 2^c}=[(H_c^{1})^{T},\ldots ,(H_c^{r})^{T}]^{T}\), where \(H_c^{i}=H_c+(i-1)2^cS_c\) for \(i=1,\ldots ,r\). The design matrix \({\textbf {X}}\) is constructed as \({\textbf {X}}=[B_{r2^c\times 2^c}^{T},-B_{r2^c\times 2^c}^{T}]^T\). The proposed method is flexible in terms of run sizes as they can be either even or odd. It accommodates more factors than that in [50] but the factor sizes are restricted to be powers of two.

[27] constructed OLHDs and NOLHDs with \(n^2\) runs and 2fp factors by coupling OAs. In practice, it starts with an OLHD(n, p) or NOLHD(n, p), which will be coupled with an OA(\(n^2,2f,n,2\)). Let \(b_{ij}\) denote the elements of OLHD(n, p), where \(i=1,\ldots ,n\) and \(j=1,\ldots ,p\), and \({\textbf {A}}\)=OA(\(n^2,2f,n,2\)). For \(j=1,\ldots ,p\), construct an \(n^2 \times (2f)\) matrix \(A_j\) from \({\textbf {A}}\) by replacing its levels \(1, \ldots , n\) with \(b_{ij},\ldots ,b_{nj}\). Partition each \(A_j\) matrix into f pieces, i.e. \(A_j=[A_{j1},\ldots ,A_{jf}]\), where each of \(A_{j1},\ldots ,A_{jf}\) contains exactly two columns. Let \(M_{j}=[A_{j1}V,\ldots ,A_{jf}V]\), where

and the design matrix \({\textbf {X}}=[M_{1},\ldots ,M_{p}]\). For example, an OLHD(11, 7), coupled with an OA(121,12,11,2), would yield an OLHD(121, 84). Similarly, a NOLHD(169, 168) can be generated from an NOLHD(13, 12) coupled with an OA(169,14,13,2). One advantage of their construction method is that it requires fewer runs to accommodate a large number of factors. The number of factors here must be even numbers, and the run size is restricted by the availability of the OAs.

[3] proposed a method to construct OLHDs with the run size n being odd primes via the Williams transformation [45]. When the factor size \(k \le (n-1)/2\), let Y denote an \(n \times k\) matrix and the elements of Y are defined as

where \(y_{ij} \in \{0, 1, \ldots , n-1 \}\), \(i=1, 2, \ldots , n\), \(j=1, 2, \ldots , k\), and \(g_1,g_2,\ldots ,g_k\) are distinct elements in the set \(\{1, 2, \ldots , (n-1)/2\}\). The design matrix would be \({\textbf {X}}=W(Y)\), where the Williams transformation W is defined by Equation (4). When \((n-1)/2+1 \le k \le (n-1)\), let \(r=k-(n-1)/2\) and the design matrix would be partitioned as \({\textbf {X}}=({\textbf {X}}_0,{\textbf {X}}_1)\), where \({\textbf {X}}_0\) is the \(n \times (n-1)/2\) design matrix generated by the same procedure as above and \({\textbf {X}}_1\) is an \(n \times r\) design matrix having the following elements before applying the Williams transformation:

\(x_{ij} \equiv ig_j\) mod n,

where \(x_{ij} \in \{0, 1, \ldots , n-1 \}\), \(i=1, 2, \ldots , n\), \(j=1, 2, \ldots , r\), and \(g_1,g_2,\ldots ,g_r\) are distinct elements in the set \(\{1, 2, \ldots , (n-1)/2\}\). Note that n can be any odd number in theory, but it is often assumed to be a prime number so that the generated designs are LHDs.

In this section, we aim to provide guidance for practitioners on how to choose appropriate method(s) for generating optimal or near-optimal LHDs. Simulations are conducted to compare different search algorithms and algebraic constructions. To ensure fairness, all algorithms and algebraic constructions were implemented in the same software environment, namely R [34]. We also make all implementations publicly available through the R package LHD ( [42], see also [43]), which enables practitioners to directly apply the methods discussed in this paper.

For each design size considered in the simulations, each search algorithm is run 20 times, and we record both the best design obtained and the average CPU time across the 20 runs. In all reported tables, a value of “0” indicates that the CPU time is less than one second. For all algorithms, the maximum number of iterations is set to 500, except in the convergence analysis shown in Figure 2, where we extended the maximum number of iterations to 2000 to better illustrate the convergence behavior of different algorithms. In addition to the best-found results reported in the tables, Figures 1 and 3 present boxplots of the objective values across repeated runs. These boxplots summarize the distribution of performance, including the median, variability, and potential outliers, and thus provide insight into the stability of each algorithm.

For the SA based algorithms (SA [31], SA2008 [22], and SLHD [2]), the temperature decrease rate is set to \(10 \%\), following [31] for generally good results.

For LaPSO, let \({ SameNumP}\) denote the number of trials used to reduce the Hamming distance with respect to the personal best for a given column within a particle, and let \({ SameNumG}\) denote the corresponding quantity with respect to the global best. We set \({ SameNumP} = 0\) and \({ SameNumG} = n/4\), as recommended by [6]. It is worth noting that setting \({ SameNumP} = n/2\) and \({ SameNumG} = 0\) could yield better designs, but requires substantially more iterations (five times as many) [6]. The tuning parameter \(p_{0}\), described in Section 3.2, is set to \(1/(k-1)\), also following [6]. We set the number of particles to 100.

For GA, the population size m is set to 100, matching the number of particles used in LaPSO. The mutation probability \(p_{\text {mut}}\), described in Section 3.3, is set to \(1/(k-1)\), consistent with the choice of \(p_{0}\) in LaPSO. We aim to have equivalent settings to ensure the fairness in comparison. All parameter settings are held fixed across all design sizes. This avoids instance-specific tuning and ensures that the comparison reflects the intrinsic performance of the algorithms rather than differences in parameter alternation.

First, we compare five search algorithms: SA [31], SA2008 [22], SLHD [2], LaPSO [6], and GA [26] together with the algebraic method FastMm [44] for constructing maximin distance LHDs under both the \(L_1\) and \(L_2\) distance criteria. All of these methods and their implementation details are described in Sections 3 and 4, respectively. We discuss the results under the \(L_1\) distance in Table 3, and those under the \(L_2\) distance in Tables 4 and 5. Note that maximin \(L_1\) distance designs can differ substantially from maximin \(L_2\) distance designs when the design sizes are small; however, for large design sizes, the resulting designs tend to be similar.

Table 3 reports the minimum \(\phi _{p}\) values [31] under the \(L_1\) distance (i.e., the Manhattan distance, corresponding to \(q=1\) in Equation (1)), along with the average CPU time for each case. Designs generated by the FastMm method are theoretically proven to be optimal when their sizes are \(n \times (n-1)\) with n being any odd prime, and \(n \times n\) with either \(2n+1\) or \(n+1\) being prime [44]. In most cases, these designs also achieve the smallest \(\phi _{p}\) values in Table 3. Consequently, all designs constructed using the FastMm method corresponding to these theoretically covered sizes are optimal and can be obtained instantly. This includes the \(12 \times 12\) case, for which \(n + 1 = 13\) is prime.

We note, however, that for the \(12 \times 12\) case, the \(\phi _{p}\) values reported by LaPSO and GA (0.0256) are slightly smaller than the value reported for FastMm (0.0258) in Table 3. This apparent discrepancy arises because the theoretical optimality of FastMm is established with respect to the \(L_1\) distance criterion, \(d_1({\textbf {X}})\), defined at the beginning of Section 2, whereas Table 3 reports the \(\phi _{p}\) criterion. Although \(\phi _{p}\) is commonly used as a surrogate measure for constructing maximin designs, it depends on the entire set of pairwise distances rather than only the minimum distance. Consequently, for a finite value of p, minimizing \(\phi _{p}\) does not necessarily induce the same ranking as maximizing \(d_1({\textbf {X}})\). For example, a design may have a slightly smaller minimum distance but a more favorable overall distance distribution, resulting in a smaller \(\phi _{p}\) value. We verified that, for the \(12 \times 12\) case, the design generated by FastMm has \(d_1({\textbf {X}}) = 50\), whereas the best designs found by GA and LaPSO have \(d_1({\textbf {X}}) = 49\). Therefore, the FastMm design remains optimal with respect to the maximin \(L_1\) distance criterion, consistent with the theoretical result, because it achieves a larger \(L_1\) distance.

Among the five search algorithms, GA achieves the smallest \(\phi _{p}\) values while also requiring the lowest average CPU time, and in some cases identifies the same optimal designs as FastMm. LaPSO performs comparably to GA in terms of \(\phi _{p}\) values but requires substantially longer CPU time. Overall, for the special design sizes covered by [44], the algebraic construction FastMm is recommended. When a search algorithm is required, GA is the preferred choice.

Since the \(L_2\) distance (i.e., the Euclidean distance, corresponding to \(q=2\) in Equation (1)) is the most popular choice in practice, we focus on maximin \(L_2\) distance designs in the following discussion. Table 4 presents the results for designs with run sizes ranging from 4 to 14. When \(k = 2\) with \(n = 4\) or \(n = 5\), all methods achieve the same minimum \(\phi _{p}\) value, while FastMm requires the least CPU time. For \(k = 2\) with \(n = 6, 7,\) or 8, all search algorithms yield the same minimum \(\phi _{p}\), and GA has the lowest CPU time. When \(k = 2\) with \(n = 10\) or \(n = 12\), both LaPSO and GA attain the same minimum \(\phi _{p}\). For the \(14 \times 2\) design, GA achieves the smallest minimum \(\phi _{p}\). When \(k = 3\) with \(n \le 7\) and \(k = 4\) with \(n \le 6\), nearly all search algorithms produce the same minimum \(\phi _{p}\), with GA again requiring the least CPU time. In all remaining cases, GA and LaPSO achieve similar minimum \(\phi _{p}\) values and outperform the other search algorithms; however, GA consistently requires the least computational time.

In Table 5, we consider designs following the rule of thumb run size (\(n = 10k\)) commonly used in computer experiments [4, 15, 20, 29]. Table 5 shows that as both n and k increase, the GA tends to require the least CPU time among the five search algorithms. In most cases, GA and LaPSO outperform the other methods in terms of the \(\phi _{p}\) values, with GA often yielding better results than LaPSO. LaPSO generally incurs the highest CPU time. Therefore, we recommend GA when computational resources are limited.

We further present boxplots of the \(\phi _{p}\) values for selected designs with rule of thumb run sizes in Figure 1. For the case of \(n = 20\) and \(k = 2\) (the top-left panel of Figure 1), GA yields the smallest median \(\phi _{p}\) value, while LaPSO appears less stable than the other algorithms. In the remaining three panels, GA and LaPSO outperform the other methods. In most cases, the relative performance patterns observed in terms of the best-found values are consistent with those of the median values, suggesting that the conclusions drawn from the tables are robust. Moreover, algorithms with tighter interquartile ranges, such as GA, tend to exhibit more stable performance across runs. As the design sizes increase, the advantage of using GA becomes increasingly obvious.

Figure 2 illustrates the convergence behavior of each algorithm for different design sizes. For the case of \(n = 20\) and \(k = 2\) (the top-left panel of Figure 2), SA converges around the 1300th iteration, whereas the remaining algorithms all converge before the 500th iteration. For the \(40 \times 4\) and \(60 \times 6\) cases (the top-right and bottom-left panels of Figure 2), SLHD continues to show improvements until approximately the 2000th iteration, while the other algorithms stop improving around the 1200th iteration. Overall, GA converges faster than all other methods while achieving strong performance, with LaPSO ranking second. Note that GA generates new candidate designs by exchanging entire columns, whereas the other methods rely on element-wise exchanges. Element-wise operations are effective when n and k are relatively small; however, as the size of designs increases, the number of elements grows exponentially, and element-wise exchange strategies may become less efficient.

We note that the computational meaning of an “iteration” differs across algorithms. For population-based algorithms (i.e., GA and LaPSO), one iteration involves updating multiple candidate designs, whereas for SA-based algorithms, each iteration updates a single design. Accordingly, the convergence curves should be interpreted qualitatively as indicators of stabilization behavior rather than as measures of computational efficiency.

[21] adopted a simulated annealing algorithm to identify MaxPro LHDs, which is implemented in the R package MaxPro [1] (denoted as “MaxPro” in the tables). Here, we compare this approach with the LaPSO and GA algorithms implemented in the R package LHD [42]. The objective is to identify MaxPro LHDs by minimizing the objective function \(\psi \) defined in Equation (2). For each design size, we run each algorithm 20 times and report the best results (i.e., the minimum \(\psi \) values) along with the average CPU time in Tables 6 and 7. For the algorithm from the MaxPro package, we set the number of iterations to 500 and keep all other input arguments at their default settings. For both LaPSO and GA, all parameter settings follow those specified at the beginning of Section 5.

Note that the MaxPro package was developed in C++ using Rcpp framework, while the LHD package was developed entirely in the standard R environment. Consequently, MaxPro achieves substantially lower CPU times in all cases. Therefore, the reported CPU times are provided for information and reference purposes only.

In Table 6, when the design sizes are small (i.e., \(k = 2\) with \(n \le 8\), \(k = 3\) with \(n \le 6\), and \(k = 4\) with \(n \le 5\)), all three algorithms yield the same minimum \(\psi \) values, while the MaxPro method has the lowest CPU time. For the remaining design sizes (except for the case of \(12 \times 3\)), the GA produces the best results. The LaPSO also outperforms the MaxPro method, with \(\psi \) values close to those obtained by the GA. In Table 7, we present additional results for designs with rule of thumb run sizes. For relatively large design sizes, the GA consistently identifies designs with the smallest \(\psi \) values. The MaxPro method requires the least CPU time, whereas the LaPSO generally outperforms MaxPro in terms of design quality but incurs the highest CPU time among the three methods.

Figure 3 presents boxplots of the \(\psi \) values for several cases with rule of thumb design sizes. In the upper left panel of Figure 3 (\(k = 2\) and \(n = 20\)), the MaxPro method appears less stable than the others, as indicated by its more dispersed boxplot. The GA identifies the best-found design and exhibits the greatest stability. In the remaining three panels of Figure 3, the GA again outperforms the other methods. As the design size increases, the advantage of the GA becomes more obvious. Generally speaking, GA is recommended for constructing MaxPro LHDs of moderate to large sizes, offering reasonable CPU times. When computational time is the primary concern, the MaxPro method is preferred.

In this subsection, we aim to identify orthogonal LHDs (OLHDs) and nearly orthogonal LHDs (NOLHDs) by minimizing the ave(|q|) or max|q| criteria defined in Equation (3). Tables 8 to 10 present results for selected designs with run sizes up to \(2^8\). In these tables, “CM” indicates whether an algebraic construction method is available for the given design size. If “CM” is “Yes,” the corresponding OLHD can be obtained using the construction method; otherwise, search algorithms are required. For a summary of available algebraic constructions across different design sizes, please refer to Table 2 in Section 4.2. For all the algorithms considered in this subsection, the parameter settings follow those specified at the beginning of Section 5.

In Table 8, algebraic constructions are available for two cases (\(8 \times 4\) and \(9 \times 4\)). SA2008 can also identify these OLHDs, with both ave(|q|) and max|q| equal to 0. For the remaining cases, SA2008 is capable of identifying NOLHDs with very small ave(|q|) and max|q| values. In general, SA2008 outperforms all other search algorithms, except for the \(10 \times 6\) and \(12 \times 6\) cases, in which GA performs better. In Table 9, we further compare SA2008 and GA for larger design sizes. Note that the SA, the SLHD and the LaPSO perform clearly worse than the SA2008 and the GA in our study, and we do not include them in Tables 9 to 10 for conciseness. When \(n \ge 16\) and \(k \ge 6\), GA provides the best results, achieving very small ave(|q|) and max|q| values. Moreover, its CPU time is substantially lower than that of SA2008.

In Table 10, we present additional cases with even larger run sizes. Only a few of these cases have available algebraic constructions. The designs generated by GA achieve max|q| values below 0.025 in all cases while requiring reasonable CPU time. Moreover, in Table 11, we consider cases with rule of thumb design sizes, comparing all five algorithms implemented in the R package LHD. For the \(20 \times 2\) design, all five algorithms successfully identify OLHDs, with GA and SA achieving the lowest CPU times. For cases with \(k \ge 3\), GA achieves the smallest max|q| values (below 0.0024) and requires the least CPU time in all instances. In general, SA2008 performs best for small designs (i.e., \(n < 16\)), while GA is recommended for moderate and large designs. This conclusion aligns with the discussion in Section 5.1: element-wise operations are effective when n and k are relatively small, while GA generates new candidate designs by exchanging entire columns, which is less affected as design sizes increase.

In this paper, we focus on constructing three commonly used types of optimal LHDs: maximin distance LHDs, maximum projection LHDs, and (nearly) orthogonal LHDs. We have summarized and compared several widely used search algorithms, including SA [31], SA2008 [22], SLHD [2], LaPSO [6], and GA [26], along with various algebraic constructions [3, 7, 27, 39, 41, 44, 50]. Our goal is to provide guidance for practitioners in selecting appropriate experimental designs. Algebraic constructions are preferred when available, particularly for large designs, whereas search algorithms are employed to generate near-optimal LHDs with flexible sizes.

From the numerical studies presented in Section 5, GA emerges as the most reliable approach to generate near-optimal LHDs, offering a good balance between performance and computational time. GA often provides competitive or superior results for moderate and large design sizes. LaPSO also performs well, but it generally requires substantially more computational time. For constructing MaxPro LHDs, the R package MaxPro [1] is computationally fast, but its performance is often inferior to that of GA implemented in the R package LHD [42]. We provide three tables summarizing the recommended methods for different selected design sizes in Tables 12, 13, 14. We note that the relative performance of these methods may depend on the design size, tuning parameters, and computational budget, and thus the recommendations should be interpreted as practical guidelines rather than definitive rankings.

The search algorithms discussed in this paper can also be adapted to generate other types of efficient experimental designs. An interesting future research is to extend the scope to include additional design types, such as fractional factorial designs [8], supersaturated designs [28], and order-of-addition designs [32].