LHS in LHS: A new expansion strategy for Latin hypercube sampling in simulation design

TL;DR

This paper introduces LHS in LHS, an expansion algorithm for Latin Hypercube Sampling that allows adding samples to existing designs while maintaining their properties, along with a new metric to measure design quality.

Contribution

It proposes a novel expansion strategy for LHS that preserves its properties and introduces a new metric for assessing LHS quality, implemented in a Python package.

Findings

The algorithm effectively adds samples to existing LHS designs.

The new metric quantifies deviation from ideal LHS properties.

Implementation is available as an open-source Python package.

Abstract

Latin Hypercube Sampling (LHS) is a prominent tool in simulation design, with a variety of applications in high-dimensional and computationally expensive problems. LHS allows for various optimization strategies, most notably to ensure space-filling properties. However, LHS is a single-stage algorithm that requires a priori knowledge of the targeted sample size. In this work, we present LHS in LHS, a new expansion algorithm for LHS that enables the addition of new samples to an existing LHS-distributed set while (approximately) preserving its properties. In summary, the algorithm identifies regions of the parameter space that are far from the initial set, draws a new LHS within those regions, and then merges it with the original samples. As a by-product, we introduce a new metric, the LHS degree, which quantifies the deviation of a given design from an LHS distribution. Our public implementation is distributed via the Python package expandLHS.

Figures

Equations

$$\mathcal{D}(S)=\frac{1}{NP}\sum _{j=0}^{P-1}\sum _{l=0}^{N-1}\min \left[\sum _{i=0}^{N-1}{I}_{\left[\frac{l}{N},\frac{l+1}{N}\right)}(S_{ij}),1\right],$$

$$I_{\left[a,b\right)}(x)=\{\begin{array}{ll}1& x\in \left[a,b\right)\\ 0& x\notin \left[a,b\right)\end{array}.$$

$$\mathcal{D}\left(S\right)=1⟺S=\mathrm{LHS}(P,N),$$

$$\mathcal{D}=1-\frac{1}{6(1+M/N{)}^3},$$

Full text

Keywords: Simulation design, Latin hypercube sampling, Space-filling

 

Code metadata

[TABLE]

††journal: Software X

\affiliation

[a]organization=Dipartimento di Fisica “G. Occhialini”, Universitá degli Studi di Milano-Bicocca, addressline=Piazza della Scienza 3, city=Milano, postcode=20126, country=Italy \affiliation[b]organization=INFN, sezione di Milano-Bicocca, addressline=Piazza della Scienza 3, city=Milano, postcode=20126, country=Italy \affiliation[c]organization=Dipartimento di Informatica, Sistemistica e Comunicazione, Universitá degli Studi di Milano-Bicocca, addressline=Viale Sarca 336, city=Milano, postcode=20126, country=Italy

1 Motivation and significance

The continuous growth of available data is a common trend across disciplines, often enabling the study of increasingly complex phenomena. A useful strategy to address these challenges is to perform simulated experiments, where one analyzes the outcomes of a model and/or the predictions of a theory while controlling the effects of input parameters. Simulations are widely used in various fields, ranging from mathematics to physics, medicine, economics, and education [1]. More often than not, simulating complex phenomena carries a high computational cost, making it critical to determine where to place simulations in the parameter space.

To this end, simulation design is a branch of statistics that aims to optimize computational experiments [2, 3]. Simulation design assists in developing algorithmic strategies that explore the parameter space effectively, ensuring an optimal distribution of samples (where optimality needs to be appropriately defined), while mitigating computational costs. A set of samples, where simulations will be placed, can therefore be preferred over another by enforcing specific properties such as space-filling, one-dimensional projection, and low correlation. There is a vast literature on simulation design, with a variety of implementations and model-independent algorithms based on different notions of optimality [2]. In this paper, we focus on one such approach: Latin hypercube sampling (LHS)[4].

LHS was first introduced in the 1970s [5, 6] and was further developed to optimize projection properties, improve space-filling design, and avoid spurious correlations [7, 8, 9, 10, 11, 12, 13, 14]. The concept of Latin square comes from combinatorial mathematics — an $N\times N$ square with $N$ different symbols appearing only once per column and row. A Latin hypercube generalizes this property to a $P$-dimensional hypercube, where each dimension is binned into $N$ disjoint intervals $[i/N,(i+1)/N)$ where $i=0,1,...,N-1$ with marginal probability $1/N$. LHS distributes points to tile the space while preserving the one-dimensional projection property—i.e., ensuring that each marginalized one-dimensional bin contains exactly one sample. Multiple designs with this property are possible for a given $N$, and one can further optimize the final result by imposing additional criteria, such as orthogonality and distance optimization [4]. However, these are typically single-stage sampling algorithms, requiring all samples to be drawn simultaneously or, at most, in a finite number of discrete steps. This implies that the number of targeted samples must be known a priori, and once the set is drawn, additional samples generally cannot be added while preserving the desired projection property.

This is a crucial limitation in many practical applications. In most realistic scenarios, simulation design does not happen all at once but rather gradually and iteratively, with subsequent batches that need to be properly initialized. Suppose one initially plans for $N$ simulations based on the available computational budget and distributes them using LHS. Later, realizing that this is insufficient for the targeted application, one might secure additional computing time for $M$ more simulations. How should these $M$ additional simulations be integrated into the initial LHS-distributed set of $N$ simulations so that the combined set of $N+M$ samples still retains desirable space-filling properties?

In this paper, we propose a new model-free expansion algorithm for LHS, accompanied by an open-source package for the Python programming language, which we dubbed expandLHS. Compared to other proposed LHS expansion strategies [15, 16, 17, 18, 19, 20], our approach allocates new samples by rebinning the original hypercube to preserve the projection property. Building on this idea, we can further leverage the versatility of Latin hypercubes by essentially designing an LHS within the existing LHS. As a by-product, this paper introduces the notion of “LHS degree” to quantify how closely a given set of samples resembles an LHS with the same number of points.

2 Software description

2.1 Algorithm

Let us denote by ${\rm LHS}(P,N)$ an initial set of $N$ points in $P$ dimensions, defined within the hypercube $[0,1)^{P}$ and satisfying the one-dimensional projection property of LHS. That is, there is one and only one sample per interval if each dimension is uniformly binned into $N$ disjoint intervals. We wish to distribute $M$ new samples to obtain an extended set ${\rm eLHS}(P,N+M)$.

The algorithm we propose consists of three main steps, as illustrated in Fig. 1:

1. (a)

Regridding. Each dimension is binned into $N+M$ equal-width, disjoint intervals. The initial samples ${\rm LHS}(P,N)$ now fall onto a different grid, which creates $Q\geq M$ empty bins and may lead to overlaps of samples in other bins. 2. (b)

Void extraction. We consider all empty intervals and select $M$ of them. 3. (c)

“LHS in LHS”: In this $M$-dimensional subgrid, we distribute $M$ samples using, once more, an LHS strategy. These are added to the initial LHS to obtain the targeted extended set ${\rm eLHS}(P,N+M)$.

The probability of obtaining multiple samples in the same interval after regridding depends on both the initial set ${\rm LHS}(P,N)$ and the value of $M$. When this happens, the final extended set ${\rm eLHS}(P,N+M)$ partially loses the one-projection property. In fact, if $Q>M$ in a certain dimension, there will necessarily be $Q-M$ bins with multiple points and, consequently, empty bins after the expansion.

To quantify this loss, we introduce a new metric $\mathcal{D}$, which we refer to as the “LHS degree” of a sample set. Consider a generic sample set $S$ with $N$ elements in $P$ dimensions; this is composed of real numbers $S_{ij}$ with $i=0,\dots,N-1$ and $j=0,\dots,P-1$. We define

[TABLE]

where $I$ is the indicator function defined as

[TABLE]

In particular, one has $0<\mathcal{D}(S)\leq 1$ and

[TABLE]

The degree $\mathcal{D}$ can thus be interpreted as the fractional closeness to an LHS set. Essentially, this procedure checks all the $i$-intervals for each dimension $j$. If they are populated by a sample $S_{ij}$, a weight $1/(NP)$ is assigned; otherwise, the weight is 0. The metric is the sum of all the weights.

After the expansion, one has $\mathcal{D}({\rm eLHS})\leq 1$. Indeed, the sums in Eq. (1) trace the presence of at least one sample per interval. An interval containing overlapping samples still contributes a weight of 1 to the weighted sum described above, due to the use of the $\min$ function. For each such interval, there exists a corresponding empty interval in the same dimension that receives a weight of zero. Consequently, the degree $\mathcal{D}({\rm eLHS})$ is reduced. Figure 1 shows an example with $N=7$, $M=3$, and $P=2$. In this case, the final degree is $\mathcal{D}=0.95$. A “perfect” expansion with $\mathcal{D}({\rm eLHS})=1$ corresponds to the case where the new $N+M$ set of points is also LHS-distributed. In general, this is not always possible.

The selection of $M$ empty intervals is not unique whenever $Q<M$. Among the possible choices, we select the one that maximizes either the centered discrepancy or the geometric discrepancy [2, 21]. The former measures the uniformity of the sample set as a proxy for its space-filling properties, while the latter tracks the minimum Euclidean distance between samples. In practice, we generate multiple expansions and select the one that minimizes (maximizes) the discrepancy (geometric discrepancy), up to a certain tolerance level. This becomes more relevant in a larger number of dimensions $P$ and as the number of new samples $M$ approaches the initial sample size $N$.

2.2 Implementation

Our software, expandLHS, is implemented in Python 3 and leverages Numpy’s arrays and Numba’s just-in-time compilation for fast computations. expandLHS is distributed via the Python Package Index and can be installed via:

Ψpip install expandLHS

Dependencies include Numpy [22], Scipy [23], and Numba (version $\geq 0.57.0$) [24]. If not present, these libraries will be installed/updated along with the package.

The expandLHS module is implemented as a single Python class. This is imported within a Python console, script, or Jupyter Notebook using e.g.

Ψfrom expandLHS import ExpandLHS

Our software is in its v1.1 release. The source code is distributed under the git version control system and the MIT permissive license at github.com/m-boschini/expandLHS. See Fig. 2 for a flowchart describing the overall implementation.

The expandLHS class has additional functionalities implemented. The method ExpandLHS.degree computes the degree $\mathcal{D}$ defined in Eq. (1). The method ExpandLHS.optimal_expansion iterates over possible values of $M$ to identify the expansion strategy that maximizes the degree of the resulting eLHS. Finally, the inner LHS in step (c) of the algorithm above can be optimized to achieve a user-defined threshold in the chosen metric. At present, our package implements two different metrics from the scipy.stats.qmc submodule: discrepancy and geometric_discrepancy.

3 Illustrative examples

3.1 LHS expansion

The code snippet below introduces the basic usage of expandLHS. We extend an initial LHS set with $N=20$ points in $P=2$ dimensions, adding $M=18$ new samples.

⬇

from scipy.stats.qmc import LatinHypercube

from expandLHS import ExpandLHS

initial LHS set

P = 2 # hypercube dimension

N = 20 # initial sample size

lhs_sampler = LatinHypercube(P)

lhs_set = lhs_sampler.random(N)

LHS expansion

M = 18 # expansion size

eLHS = ExpandLHS(lhs_set)

elhs_set = eLHS(M)

elhs_opt = eLHS(M, optimize=’discrepancy’)

Outputs are illustrated in Fig. 3, where we show the initial set, an unoptimized expanded set, and the optimized expanded set. Optimizing using the discrepancy criterion increases the uniformity of the expanded sample set. For this case, we report a discrepancy estimate of $1.0\times 10^{-3}$ for the unoptimized expansion and $4.7\times 10^{-4}$ for the optimized expansion. This particular case happens to be a perfect expansion, $\mathcal{D}({\rm eLHS})=1$.

3.2 How many new samples?

In the following example, we simulate different expansions for the same initial set by adding a number of new points $M$ between 4 and 12, ranking our attempts according to the resulting degree $\mathcal{D}$. Note that this metric depends only on the initial set and the value of $M$, and not on the added samples.

⬇

optimal expansion size with M in [4, 12]

verbose = False returns the optimal expansion

verbose = True returns all the estimates

eLHS.optimal_expansion((4,12), verbose=True)

M degree

(0, 1.0), # no expansion

(12, 1.0), # best option

(9, 0.9828),

(7, 0.9815),

(10, 0.9667),

(6, 0.9615),

(11, 0.9516),

(5, 0.92),

(4, 0.9167),

(8, 0.9107) # worst option

We obtain a perfect expansion for $M=12$, while the worst case is that with $M=8$. This functionality offers some interesting use cases. In the example above, suppose one has budgeted computational time for $M=8$ new simulations. It turns out a minor increase to $M=9$ ensures substantially better coverage of the parameter space, with $\mathcal{D}$ increasing from $\sim 0.91$ to $\sim 0.98$.

3.3 General behavior of eLHS

We now present some general properties of our expansion algorithm with a particular focus on the degree metric introduced in Sec. 2. Figure 4 shows the degree $\mathcal{D}$ as a function of the expansion size $M$, for initial LHS sets with different numbers of samples $N$ and parameter-space dimensions $P$. For stability, results are averaged over multiple realizations.

Lower values of $M$ correspond to a higher likelihood of overlaps, leading to a lower degree and thus a less space-filling eLHS. This effect diminishes as $M$ increases; when $M>N$, the expansion dominates over the initial set, and the degree approaches 1. Figure 4 shows “spikes” with $\mathcal{D}=1$ taking place at multiples of the initial set size, i,e. $M=kN$ with $k\in\mathbb{N}$. This is not surprising: in those expansions, the existing boundaries of each bin remain unchanged and each initial interval is divided into $k$ subintervals. This implies that overlaps never occur, resulting in a perfect expansion. Hints of this behavior were already presented in Refs. [18, 25], even though the authors only considered the case with $M=N$. Less trivially, Fig. 4 shows that expansions with sizes $M=\left(k+\frac{1}{2}\right)N$ exhibit higher degrees than adjacent values, although they do not reach $\mathcal{D}=1$.

We find that the LHS degree is largely independent of the number of dimensions $P$ and presents a self-similar behavior with $N$. In other words, $\mathcal{D}$ depends on the ratio $M/N$ and not on any of these two parameters independently. In particular, the degree is well predicted by the following phenomenological expression111More precisely, a least-square fit to the curves of Fig. 4 with ansatz $\mathcal{D}=1+a(b+M/N)^{c}$ returns $a=-0.167$, $b=1.01$, and $c=-2.99$.

[TABLE]

which, however, does not capture the high-$\mathcal{D}$ spikes described above.

Finally, the middle panel of Fig. 4 compares our results against repeated unitary expansions, that is, expanding the LHS $M$ times with one new sample at a time, instead of performing a direct one-step expansion of size $M$ as considered here. As expected, a unitary expansion strategy is highly suboptimal, as it leads to an overall increase in the likelihood of overlaps. Furthermore, it drastically increases the computational time by a factor $\propto M$.

3.4 Computational performance

Figure 5 illustrates the computational cost of expandLHS as a function of $M$ and $P$. Just-in-time compilation via Numba significantly speeds up the computation, particularly when the expansion is performed without optimization. Unsurprisingly, the execution time increases when optimization is required. In particular, optimization degrades the algorithm’s scaling with the dimensionality of the hypercube.

Both of our optimization schemes show a similar trend as the number of samples grows, with discrepancy performing slightly better than geometric discrepancy. This scaling behavior can be explained as follows. The discrepancy focuses on uniformly filling the hypercube and quantifies the difference between the input distribution of points and the expected uniform coverage. Its evaluation requires comparing the fraction of samples in a hypercube subset with the fraction of volume it occupies, which depends heavily on $P$ and mildly on $M$. On the other hand, the geometric discrepancy estimates the minimum distance between any pair of points and is thus mainly affected by the number of samples $N+M$.

4 Impact

LHS plays a major role in simulation design, as evidenced by the vast literature on the topic (for a review see Ref. [4]). The additional flexibility provided by our expansion algorithms broadens the use cases of this powerful tool. One interesting application is in training machine learning models. A model initially trained on a set of $N$ samples may fail to achieve the required accuracy. A common solution in machine learning is to increase the training set size or extend the coverage of the input space. An expansion algorithm for LHS can assist in selecting an appropriate distribution for new training samples while preserving existing information. A study on this topic is presented in Ref. [26]. To this end, expanded LHS might complement the usage of learning curves in machine learning [27]. Quoting the specific field of research of some of us, gravitational-wave astronomy, two activities to which we plan to apply expandLHS are the development of numerical-relativity surrogate models [28, 29] and the interpolation of stellar-physics simulations [30, 31].

Our expandLHS algorithm shares some similarities with the method presented in Ref. [20]. The key differences are that our approach assumes fewer constraints on the final samples and directly targets the required number of additional samples $M$ while allowing for expansions that are not perfect Latin Hypercube sets. Reference [18] also presented a related algorithm, which they dubbed Progressive LHS, supported by a public Python implementation. Their approach samples an LHS while also ensuring additional properties derived from Sliced Latin Hypercube sampling [32]. In particular, they construct a sequence of LHS slices whose progressive union still forms a Latin hypercube. The complete set obtained by the union of all slices constitutes a Latin hypercube that maximizes space-filling properties. Crucially, this remains a one-stage sampling procedure, meaning that the size of the final Latin hypercube and the number of slices must be specified a priori. Moreover, the algorithm permits sampling only those configurations in which the total number of samples is an integer multiple of the number of slices. Our implementation is different; the regridding and “LHS in LHS” operations we propose are arguably closer to the original procedure for distributing points via LHS and, by relaxing the requirement of a perfect Latin hypercube, we allow for expansions of any size.

5 Conclusions

In this paper, we presented a new algorithm for expanding an existing LHS while optimally preserving its space-filling properties. Our procedure leverages the flexibility of LHS by embedding a new LHS within the initial one. Loss of optimality arises from the potential occurrence of overlapping samples in the regridded space, as quantified by the LHS degree $\mathcal{D}$.

Our new algorithm has broad applications in simulation design and machine learning, assisting in scenarios where a carefully planned initial simulation must later be expanded to meet targeted requirements. However, the expansion strategy introduced in Sec. 2 is specialized in sampling the uniform hypercube. Other LHS-like algorithms have been proposed to sample constrained spaces, handle non-uniform boundaries, and incorporate weighted samples (e.g. [33, 34, 35]). Generalizing our implementation to these non-standard LHS techniques is left to future work.

Our implementation, expandLHS, is available as an open-source package for the Python programming language, utilizing fast array manipulation and just-in-time compilation.

Acknowledgements

M.B., D.G., A.C. and M.F. are supported by ERC Starting Grant No. 945155–GWmining, Cariplo Foundation Grant No. 2021-0555, MUR PRIN Grant No. 2022-Z9X4XS, MUR Grant “Progetto Dipartimenti di Eccellenza 2023-2027” (BiCoQ), and the ICSC National Research Centre funded by NextGenerationEU. D.G. is supported by MSCA Fellowship No. 101064542–StochRewind, MSCA Fellowship No. 101149270–ProtoBH, and MUR Young Researchers Grant No. SOE2024-0000125. Computational work was performed at CINECA with allocations through INFN and Bicocca.