Latin Hypercube Sampling Overview

• Latin Hypercube Sampling (LHS) is a method that creates space-filling designs by exactly stratifying each input dimension into equal-probability intervals.
• It is widely applied in uncertainty quantification, computer experiments, and numerical integration to efficiently reduce variance in estimating function means.
• Advanced adaptations of LHS, such as LHSD and constraint-based methods, extend its use to models with dependent inputs and complex sampling requirements.

Latin hypercube sampling (LHS) is a randomized space-filling design for generating $N$ points in $[0,1]^d$ with the property that each coordinate projection forms an exact one-dimensional stratification. By construction, in each marginal, the $N$ samples occupy each of the $N$ equal-probability intervals exactly once. LHS is widely employed in uncertainty quantification, computer experiments, numerical integration, sensitivity analysis, sequential design, and high-dimensional simulation contexts. Its primary appeal is the strong one-dimensional stratification, efficient variance reduction for mean estimation of smooth or low-dimensional functions, and algorithmic simplicity.

1. Construction and Formal Definition

Let $N,d\in\mathbb N$. LHS is constructed by, for each coordinate $j=1,\ldots,d$, independently drawing a random permutation $T_j$ of $\{1,\ldots,N\}$ and independent uniforms $u_{n,j}\sim \mathrm{Uniform}[0,1)$ for $n=1,\ldots,N$. The $n$th point is

$$X_n = (X_{n,1}, \ldots, X_{n,d}), \quad X_{n,j} = \frac{T_j(n)-1 + u_{n,j}}{N}.$$

Each marginal projection is a permutation of the $N$ strata $[0,1/N), [1/N,2/N), ..., [(N-1)/N,1)$. This construction ensures exact one-per-stratum placement along all $d$ axes. The design can be discretized for factor-type parameters by replacing $[0,1)$ intervals with categorical levels (Doerr et al., 2017, Donovan et al., 2015).

2. Discrepancy and Negative Dependence Properties

The uniformity of LHS is quantified by the star discrepancy:

$$D^*_N(P)=\sup_{x\in[0,1]^d}\left|\frac{1}{N}\left|P\cap[0,x)\right| - \prod_{i=1}^d x_i\right|.$$

Probabilistic bounds establish that, with high probability,

$D^*_N(X) = \Theta(\sqrt{d/N})$

matching the rate for i.i.d. (Monte Carlo) sampling; explicit constants are $2.6442$ for LHS and $2.5287$ for MC. Lower and upper bounds coincide up to these factors (Doerr et al., 2017, Gnewuch et al., 2021). LHS points are not independent, but indicators of sample inclusion for axis-aligned boxes satisfy $\gamma$-negative dependence with $\gamma\leq e^d$, providing Chernoff-Hoeffding–type exponential concentration for empirical measures—albeit with exponential-in-$d$ constants (Doerr et al., 2021).

3. Space-Filling Metrics and Optimization

Space-filling designs are evaluated via pairwise separation (maximin criterion), discrepancy (star, $L^2$, centered), and minimum spanning tree statistics (Damblin et al., 2013).

• Maximin LHS: Among random LHS, select the design maximizing minimal Euclidean distance.
• Discrepancy optimized LHS: Select LHS with minimal centered or wrap-around discrepancy.
• Algorithmic approaches: Simulated annealing and stochastic evolutionary algorithms iteratively swap columns or entries, ensuring marginal stratification and improving global metrics.

Well-optimized LHS designs can achieve superior two-dimensional and higher-order projection uniformity, but trade-offs exist between global distance and marginal coverage (Damblin et al., 2013).

4. Statistical Properties: Estimation and Variance Reduction

For integral estimation $\mu = \mathbb E[f(X)]$, LHS yields unbiased estimators $\hat\mu_{LHS} = \frac1N \sum_{i=1}^N f(X_i)$. Variance reduction is explained by the Hoeffding/ANOVA decomposition:

$f(x) = f_0 + \sum_{i=1}^d f_i(x_i) + r(x),$

with $f_i$ zero-mean. Under LHS, variance contributions from all main effects are eliminated, leaving only interaction terms and $O(1/N)$ remainders, so:

• If $f$ is additive: $\operatorname{Var}(\hat\mu_{LHS}) = O(N^{-2})$.
• For general $f$, $$\operatorname{Var}(\hat\mu_{LHS}) \leq \operatorname{Var}(\hat\mu_{MC})$$ with $O(N^{-1})$ scaling (Kucherenko et al., 2015, Hakimi, 10 Feb 2025).

A central limit theorem (CLT) applies:

$$\sqrt{N}(\hat\mu_{LHS}-\mu) \Rightarrow \mathcal N(0, \sigma_{LHS}^2)$$

with reduced asymptotic variance compared to i.i.d. sampling (Hakimi, 10 Feb 2025). This also extends to Z-estimators and quasi-score statistics. For sequential stratification, combining LHS in each stratum with optimal weighting yields variance decay rates as strong as $O(N^{-2})$ in some adaptive procedures (Krumscheid et al., 2023).

5. Extensions to Dependence, Constraints, and Sequential Design

Standard LHS stratifies each margin but does not preserve input dependencies or operate natively under general constraints. Key research advances include:

• LHS with dependence (LHSD): Retains marginal stratification but reconstructs rank-based dependence (copula) among coordinates, providing variance reduction whenever $f$ is monotone and the copula has positive dependence (Aistleitner et al., 2013).
• Quantization-based LHS (Q-LHS): Applies Voronoi vector quantization to blockwise-dependent variables, then samples per cell, preserving full input dependence for space-filling designs. Unbiasedness and variance advantages hold, especially for block-structured models (Lambert et al., 16 May 2024).
• Constraint handling (CASTRO): For mixture or synthesis constraints (e.g., $\sum x_k=1$, combinatorial feasibility), CASTRO decomposes into low-dimensional LHS subproblems and recombines via Cartesian products, followed by candidate pruning to maximize coverage and uniformity. Distance-based selection can incorporate prior experimental data (Schenk et al., 23 Jul 2024).
• Sequential and adaptive LHS: Hierarchical stratification with local LHS and adaptive partitioning achieves global variance reduction consistent with optimal weighting arguments (Krumscheid et al., 2023). Recent methods for expansion ("LHS in LHS") merge new LHS points in void regions, quantified by the LHS-degree metric, providing near-LHS uniformity for extensibility (Boschini et al., 29 Aug 2025).

6. Coverage, Orthogonality, and Comparison with Other Designs

Coverage of $d$-dimensional discrete grids (strata per axis) after $k$ LHS trials is given by

$P(k,n) = 1 - (1-n^{-(d-1)})^k,$

approximated by $1-\exp(-k n^{-(d-1)})$ for large $n$; for $t$-dimensional projections:

$P_t(k,n) = 1 - (1-n^{-(t-1)})^k,$

independent of ambient $d$. This guides trial counts for desired subspace coverage in uncertainty quantification and model population design (Donovan et al., 2015, Burrage et al., 2015).

Orthogonal Sampling (OS): Guarantees exact uniform coverage of $t$-dimensional projections at the sub-block level. For full design coverage, LHS and OS are combinatorially equivalent when $n=p^d$. OS requires construction of suitable orthogonal arrays, often more restrictive but preferred when precise subspace uniformity is critical (Donovan et al., 2015, Burrage et al., 2015).

7. Practical Applications and Limitations

LHS is universally adopted in computer experiments, global sensitivity analysis (HSIC/statistical screening), surrogate modeling, uncertainty propagation in engineering, and efficient adversarial query strategies (Wang et al., 2022). In practice,

• Basic LHS is preferred for $N$ small to moderate and $d$ moderate.
• Optimized LHS (maximin/discrepancy minimizing) is used for high accuracy in model fitting; computational costs grow rapidly with $N,d$.
• For problem structures involving constraints or dependencies, state-of-the-art variants (e.g., Q-LHS, CASTRO, LHSD) are essential for maintaining statistical properties.

Limitations include the inability to achieve low star discrepancy rates beyond $\Theta(\sqrt{d/N})$ in high dimensions, difficulties in sequential expansion without loss of stratification, and challenges in handling strong input dependence or combinatorial constraints without specialized algorithms (Doerr et al., 2017, Boschini et al., 29 Aug 2025, Lambert et al., 16 May 2024, Schenk et al., 23 Jul 2024).

---

References:

• (Doerr et al., 2017): Probabilistic Lower Bounds for the Discrepancy of Latin Hypercube Samples
• (Gnewuch et al., 2021): Discrepancy Bounds for a Class of Negatively Dependent Random Points Including Latin Hypercube Samples
• (Doerr et al., 2021): On Negative Dependence Properties of Latin Hypercube Samples and Scrambled Nets
• (Donovan et al., 2015): Estimates of the coverage of parameter space by Latin Hypercube and Orthogonal sampling
• (Damblin et al., 2013): Numerical studies of space filling designs: optimization of Latin Hypercube Samples and subprojection properties
• (Lambert et al., 16 May 2024): Quantization-based LHS for dependent inputs: application to sensitivity analysis of environmental models
• (Boschini et al., 29 Aug 2025): LHS in LHS: A new expansion strategy for Latin hypercube sampling in simulation design
• (Aistleitner et al., 2013): A central limit theorem for Latin hypercube sampling with dependence and application to exotic basket option pricing
• (Krumscheid et al., 2023): Sequential Estimation using Hierarchically Stratified Domains with Latin Hypercube Sampling
• (Hakimi, 10 Feb 2025): Robust estimation with latin hypercube sampling: a central limit theorem for Z-estimators
• (Kucherenko et al., 2015): Exploring multi-dimensional spaces: a Comparison of Latin Hypercube and Quasi Monte Carlo Sampling Techniques
• (Shields et al., 2015): The generalization of Latin hypercube sampling
• (Park et al., 2020): Variance Reduction for Sequential Sampling in Stochastic Programming
• (Wang et al., 2022): Query-Efficient Adversarial Attack Based on Latin Hypercube Sampling
• (Schenk et al., 23 Jul 2024): A Novel Constrained Sampling Method for Efficient Exploration in Materials and Chemical Mixture Design
• (Burrage et al., 2015): Populations of models, Experimental Designs and coverage of parameter space by Latin Hypercube and Orthogonal Sampling