Space-Filling Regularization

• Space-Filling Regularization is a design principle that enforces nearly uniform coverage across input spaces, latent states, or loss fields to prevent clustering and gaps.
• It employs both soft (objective-based) and hard (constraint-based) formulations, optimizing metrics such as maximin distances, discrepancy, and Arrwwid numbers.
• Techniques include surrogate adjustment, privacy sets, PDE-based loss diffusion, and learned space-filling curves to maintain structural uniformity.

Space-filling regularization denotes a family of strategies that make samples, latent trajectories, loss fields, or recursive traversals cover a domain in a controlled and nearly uniform way. In the current literature, the phrase is explicit in nonlinear state-space identification, where it penalizes collapse of the state trajectory in extended input/state space (Klein et al., 10 Jul 2025). Closely related formulations appear as hard feasibility constraints in experimental design, as conditioning-based control of Kriging kernels when the sample set is not space-filling, as elliptic regularization of the loss field over the data domain, and as multiscale orderings induced by space-filling curves (Benková et al., 2015, Peri, 2023, Hasan et al., 4 Mar 2025, Wang et al., 2022). Across these settings, the common objective is to suppress clustering, holes, unstable interpolation, or poor locality by imposing coverage, dispersion, or connected traversal structure.

1. Conceptual scope

In computer experiments, a space-filling design is classically understood as a finite point set in a bounded region, typically $[0,1]^k$, that fills the region “as uniformly as possible” (Lin et al., 2022). The original motivation is deterministic simulation: repeated runs at the same input are unnecessary, so the design objective shifts from replication and randomization toward coverage, dispersion, and low redundancy. The same intuition extends beyond design points. In state-space identification, the regularized object is the model-induced trajectory in the extended variable $\tilde u=[x,u]^T$, and the goal is to prevent optimization from collapsing that trajectory into a narrow subset of the available region (Klein et al., 10 Jul 2025). In elliptic loss regularization, the regularized object is the loss landscape $u(X,y)$ over a domain $\mathcal D\subset \mathbb R^{d+k}$, so that supervision propagates through unsampled interior regions rather than remaining concentrated on a discrete empirical support (Hasan et al., 4 Mar 2025).

A useful taxonomy separates “soft” and “hard” formulations. In the terminology of privacy-set designs, soft methods encode space-fillingness in the objective function $\Phi(\xi)$, whereas hard methods encode it in the feasible set $\Gamma\subseteq \Xi$ itself (Benková et al., 2015). This distinction is broad enough to include maximin and discrepancy optimization, privacy-set constraints, split-plot constructions under randomization restrictions, conditioning-based surrogate regularization, and PDE-based loss regularization. The shared structure is regularizing pressure toward uniform coverage, but the mathematical carrier of that pressure differs: objective, constraint, kernel, trajectory penalty, or recursive ordering.

The phrase therefore does not denote a single canonical algorithm. It names a design principle that recurs in several technical guises: coverage of the input region, coverage of the output manifold, coverage of the latent state space, coverage of the loss field, or locality-preserving linearization of a higher-dimensional domain.

2. Design-theoretic formulations

The classical design-theoretic core is the space-filling design literature for computer experiments. Latin hypercube designs provide exact one-dimensional stratification: in each coordinate, the $n$ points occupy the $n$ equal intervals exactly once. McKay, Beckman, and Conover showed that if $f$ is monotonic in each input variable, then $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$, and Stein refined the asymptotics by showing that LHDs remove main-effect variation from Monte Carlo error (Lin et al., 2022). However, one-dimensional balance alone is insufficient. The same chapter organizes the main space-filling criteria into distance-based criteria such as maximin, minimax, Audze–Eglais, and Morris–Mitchell’s $\tilde u=[x,u]^T$0; correlation-based criteria such as $\tilde u=[x,u]^T$1 and $\tilde u=[x,u]^T$2; and discrepancy criteria such as star discrepancy, $\tilde u=[x,u]^T$3 discrepancy, and centered or symmetric $\tilde u=[x,u]^T$4 discrepancies (Lin et al., 2022). In this setting, regularization is a geometric prior toward marginal balance, pairwise repulsion, projection quality, and global uniformity.

Hard-constrained formulations make the regularization explicit at the feasibility level. Privacy sets assign to each point $\tilde u=[x,u]^T$5 a forbidden neighborhood $\tilde u=[x,u]^T$6 with $\tilde u=[x,u]^T$7; a design $\tilde u=[x,u]^T$8 is permissible if $\tilde u=[x,u]^T$9 and no distinct $u(X,y)$0 violate $u(X,y)$1 (Benková et al., 2015). This framework subsumes exact no-replication designs with $u(X,y)$2, Latin-hypercube-like coordinate exclusions, time-separation constraints, and Bridge designs with

$u(X,y)$3

The design problem becomes

$u(X,y)$4

so space-fillingness acts as a hard regularizer on admissible geometry (Benková et al., 2015).

Restricted randomization introduces a further layer. In split-plot design, whole-plot factors are hard to change and subplot factors are easy to change. The split-plot Fast Flexible Filling construction adapts Lekivetz and Jones by first clustering a large random cloud in the full space to obtain $u(X,y)$5 representative runs, and then reclustering those runs in the whole-plot factor subspace to produce $u(X,y)$6 whole plots under Ward’s criterion

$u(X,y)$7

The result is a design-side regularization toward broad coverage of the joint space, the whole-plot subspace, and the subplot subspace while respecting split-plot structure (Muehlenstaedt et al., 2020).

A more radical modification arises for nonlinear simulators. If a deterministic model is a map $u(X,y)$8, a design that is uniform in parameter space need not be uniform on the output manifold $u(X,y)$9. The manifold-aware design literature therefore replaces parameter-space uniformity by output-space uniformity with respect to Hausdorff measure, using the Jacobian factor

$\mathcal D\subset \mathbb R^{d+k}$0

and target density proportional to $\mathcal D\subset \mathbb R^{d+k}$1 (Rhee et al., 2017). The derivative-free algorithm estimates local output density through $\mathcal D\subset \mathbb R^{d+k}$2-nearest-neighbor radii and reweights sparse and dense output regions accordingly. In this formulation, space-filling regularization shifts from the input coordinates to the scientifically relevant image manifold.

The same tension between structural interpretability and geometric coverage appears in screening. Space-Filling One-Factor-At-A-Time designs remain within the OFAT class needed for screening, but optimize a GP-based space-filling criterion

$\mathcal D\subset \mathbb R^{d+k}$3

under the multiplicative inverse multiquadric kernel. The design preserves one-factor-at-a-time contrasts for total Sobol’ screening while improving global coverage relative to MOFAT (Yu et al., 1 Jun 2026). Here the regularization is not an additive penalty but optimization within a constrained design family.

3. Surrogate, latent-state, and loss-field regularization

A direct model-side formulation appears in Kriging on irregularly spaced data. Ordinary Kriging assumes interpolation weights derived from a self-correlation matrix $\mathcal D\subset \mathbb R^{d+k}$4, but adaptive optimization, digital twins, response-surface modeling from operational data, and irregular time series often generate non-space-filling samples: dense local clusters and neglected regions. The reported pathologies are clustered samples, neglected or sparsely sampled regions, loss of interpolation quality, spurious noise or high-frequency oscillations, and ill-conditioning of the matrix inversion (Peri, 2023). The proposed remedy does not add a nugget, penalty term, or validation-based criterion. Instead it fixes $\mathcal D\subset \mathbb R^{d+k}$5 in the exponential-product kernel

$\mathcal D\subset \mathbb R^{d+k}$6

and chooses the kernel parameters by solving

$\mathcal D\subset \mathbb R^{d+k}$7

The practical procedure perturbs an initial $\mathcal D\subset \mathbb R^{d+k}$8, evaluates condition numbers, selects the best candidate, and then applies a Direct Search method of the NEWUOA / Powell type for local refinement (Peri, 2023). This is a surrogate-side compensation for a fixed, non-space-filling design: the sample geometry is not changed, but the interpolation operator is regularized through conditioning.

The state-space formulation is more explicit. In Local Model State Space Networks, the model is

$\mathcal D\subset \mathbb R^{d+k}$9

$\Phi(\xi)$0

with LOLIMOT organizing local affine models over the extended space $\Phi(\xi)$1 (Klein et al., 10 Jul 2025). The baseline fit

$\Phi(\xi)$2

does not constrain the geometry of the latent trajectory, so optimization can deform that trajectory severely. The paper introduces three indicators of state-space coverage—Convex Hull Volume, mean minimum distance to a grid $\Phi(\xi)$3, and Kullback–Leibler divergence—and uses the grid-based quantity

$\Phi(\xi)$4

as the regularizer. Two penalties are proposed: $\Phi(\xi)$5 and

$\Phi(\xi)$6

On the Bouc–Wen Hysteretic System, the target-based version with $\Phi(\xi)$7 and $\Phi(\xi)$8 prevents instability on test data and allows deeper LOLIMOT refinement (Klein et al., 10 Jul 2025). In this line of work, space-filling regularization acts on the distribution of latent states rather than on the static design.

Elliptic loss regularization moves the regularized object one level further outward, from states to the loss field itself. The method defines a loss landscape $\Phi(\xi)$9 over a domain $\Gamma\subseteq \Xi$0 and imposes the elliptic condition

$\Gamma\subseteq \Xi$1

with boundary condition $\Gamma\subseteq \Xi$2 on $\Gamma\subseteq \Xi$3 (Hasan et al., 4 Mar 2025). The central theoretical statement is a maximum principle: interior loss values are bounded between the minimum and maximum boundary losses. Practically, the PDE is approximated through Brownian-bridge trajectories between data points, so training minimizes loss along stochastic paths through the interior of the domain rather than only at empirical samples. This makes the method a PDE-based space-filling regularizer for the loss field: supervision is diffused through the domain, and sparse interior regions are controlled through elliptic structure rather than left unconstrained.

4. Space-filling curves, scalarization, and locality

A distinct branch of the literature regularizes by imposing a recursive ordering on space. The pandimensional recurrence framework constructs continuous maps

$\Gamma\subseteq \Xi$4

and inverse scalarizations

$\Gamma\subseteq \Xi$5

from serpentine Hamiltonian paths on $\Gamma\subseteq \Xi$6, unifying Peano, Hilbert, and higher-dimensional generalizations (Jaffer, 2014). The integer recurrences $\Gamma\subseteq \Xi$7 and $\Gamma\subseteq \Xi$8, the limit representation

$\Gamma\subseteq \Xi$9

and the explicit inverse recurrences provide a multiscale mechanism for scalarizing a hypercube while preserving recursive structure. The paper analyzes isotropy and dimension reduction and reports that Hilbert curves perform somewhat better than Peano curves and their isotropic variants for dimension reduction (Jaffer, 2014). In regularization terms, this line furnishes deterministic, locality-aware 1D parametrizations of higher-dimensional domains.

The Fibonacci space-filling curve provides a non-dyadic alternative. It is built from the Cartesian product of the 1D Fibonacci substitution with itself, producing a 2D substitution $n$0 on four rectangle prototiles, then refining it to a decorated substitution $n$1 on 24 decorated tiles so that local traversals concatenate into a single curve (Ozkaraca, 2024). The resulting limit map is a continuous surjection

$n$2

with connectedness in the sense that consecutive rectangles share an edge. The construction induces an order on the $n$3 rectangles of each level-$n$4 partition, and the interval partition is chosen so that interval length equals rectangle area. This suggests a particularly structured form of regularization: the square is linearized by a connected, multiscale, area-adapted ordering.

Neural Space-filling Curves replace fixed orderings by learned ones. The method covers an $n$5 image graph by $n$6 circuits, constructs a dual graph $n$7, predicts edge weights $n$8, computes an MST on $n$9, and merges the circuits into a Hamiltonian circuit and then a path $n$0 (Wang et al., 2022). The learned order is optimized for downstream objectives such as lag-$n$1 autocorrelation

$n$2

or LZW code length. The method therefore regularizes the 2D-to-1D conversion by topology, Hamiltonian validity, and dataset-level sharing of the traversal.

The locality of such orderings can be measured by fragmentation criteria. The Arrwwid number of a recursive tiling or space-filling curve is the smallest $n$3 such that any ball $n$4 can be covered by at most $n$5 tiles or curve fragments with total volume at most $n$6 for some constant $n$7 (Haverkort, 2010). This yields sharp lower and upper bounds: in two dimensions, recursive tilings and curves can achieve Arrwwid number $n$8, whereas Hilbert, Morton/Lebesgue, and Peano have Arrwwid number $n$9; in three dimensions, any uniform cube tiling has Arrwwid number $f$0, and any curve based on a uniform cube tiling has Arrwwid number at least $f$1 (Haverkort, 2010). The criterion formalizes a regularization objective of low worst-case fragmentation.

The geometry of the underlying tiles matters. Recursive reptilings of acute triangles are too rigid to support face-continuous space-filling curves: no such curve can be based on a reptiling of an acute triangle (Gottschau et al., 2016). The obstruction is combinatorial and geometric: acute triangles admit only very restricted fan and cap structures, and those restrictions prevent the conforming Hamiltonian paths required for recursive face continuity. This establishes that space-filling regularization by recursive ordering is domain-dependent rather than purely formal.

5. Objective functions and algorithmic mechanisms

Across the literature, space-filling regularization is implemented by a small number of recurrent objective forms: repulsion and covering metrics, conditioning criteria, distribution-matching criteria, PDE constraints, and traversal-fragmentation metrics.

The supplementary material for sliced mixture design gives a recent distribution-matching formulation based on energy distance. For a design $f$2 and target distribution $f$3 on the constrained mixture region $f$4, the energy distance is

$f$5

The same supplementary material states a sliced criterion

$f$6

which regularizes both the full design and each slice, together with one-shot and sequential MM procedures whose limit points are stationary solutions of the finite-sample optimization problem (Xiong et al., 26 Sep 2025).

A concise comparison of representative objective forms is given below.

Object regularized	Criterion or mechanism	Representative source
Design points in $f$7	maximin, minimax, $f$8, discrepancy, orthogonality	(Lin et al., 2022)
Feasible design geometry	privacy sets $f$9 and constrained $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$0	(Benková et al., 2015)
Kriging kernel matrix	$$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$1	(Peri, 2023)
Latent state trajectory	$$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$2, CHV, KLD; penalties on $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$3	(Klein et al., 10 Jul 2025)
Loss field on $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$4	elliptic PDE $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$5	(Hasan et al., 4 Mar 2025)
Sliced mixture design	energy distance $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$6	(Xiong et al., 26 Sep 2025)
Screening design geometry	$$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$7 under MIM GP kernel	(Yu et al., 1 Jun 2026)
SFC fragmentation	Arrwwid number	(Haverkort, 2010)

Algorithmically, these criteria lead to very different optimizers. Space-filling design chapters emphasize simulated annealing, columnwise-pairwise exchange, genetic algorithms, branch-and-bound, iterated local search, particle swarm optimization, and threshold accepting for distance or discrepancy optimization (Lin et al., 2022). Privacy-set design uses the Privacy Sets Algorithm, which combines greedy augmentation and mutation by temporary privacy-set violations followed by repair (Benková et al., 2015). Split-plot FFF uses two stages of Ward clustering (Muehlenstaedt et al., 2020). Kriging regularization uses candidate perturbations followed by NEWUOA-type direct search (Peri, 2023). LMSSN space-filling regularization is embedded inside quasi-Newton training with a BFGS Hessian approximation (Klein et al., 10 Jul 2025). Elliptic loss regularization uses Brownian bridges and Euler–Maruyama simulation (Hasan et al., 4 Mar 2025). Learned scan orders use CNN/GNN weight generators, evaluators, MST extraction, and alternating optimization (Wang et al., 2022). The unifying theme is geometric control, but the operational mathematics ranges from combinatorial design search to PDE surrogates and graph neural networks.

6. Trade-offs, limitations, and open questions

A first limitation is epistemic rather than computational: improved coverage does not create information where none exists. In condition-number-based Kriging regularization, local error maps show that regularization reduces interpolation error over most of the domain but can increase local error in poorly sampled corners; the method “cannot improve the content of information, but only reshape interpolation more regularly” (Peri, 2023). The same warning appears in manifold-aware design: a parameter-space-uniform design may still be highly nonuniform on the output manifold, so “space-filling” is only meaningful relative to the space one chooses to regularize (Rhee et al., 2017).

A second limitation is that stronger regularization usually induces a fit–coverage trade-off. In nonlinear state-space models, penalizing $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$8 or enforcing a target $$\operatorname{Var}(\hat\mu_{lhs})\le \operatorname{Var}(\hat\mu_{srs})$$9 improves robustness and interpretability, but the output-error loss $\tilde u=[x,u]^T$00 no longer reaches the same low values as in the unregularized fit (Klein et al., 10 Jul 2025). Elliptic loss regularization similarly controls interior behavior inside the chosen domain $\tilde u=[x,u]^T$01, often taken as the convex hull of the data, but it is not a guarantee for arbitrary far-out-of-support points (Hasan et al., 4 Mar 2025). Space-filling regularization is therefore strongest as an interior or reachable-region prior, not as unrestricted OOD certification.

A third limitation concerns metric mismatch. The Arrwwid number is a precise worst-case fragmentation metric, but low Arrwwid number does not necessarily imply better average-case performance, and the original paper explicitly notes that preliminary experiments showed no large performance gap among several 2D curves despite Arrwwid differences (Haverkort, 2010). Likewise, the Fibonacci construction proves continuity, surjectivity, connectedness, and a multiscale ordering, but does not provide distortion bounds, Hölder exponents, or bi-Lipschitz locality guarantees (Ozkaraca, 2024). The regularization one obtains from a space-filling curve is therefore structural rather than metrically complete.

A fourth limitation is feasibility. Hard-constrained methods depend on admissible set geometry. Privacy-set algorithms assume that maximal permissible designs have size $\tilde u=[x,u]^T$02, an assumption that can fail under strong constraints or in constrained regions (Benková et al., 2015). Split-plot space-filling designs must simultaneously satisfy coverage and restricted randomization (Muehlenstaedt et al., 2020). Acute triangles admit reptilings, but not the cap/fan structure needed for face-continuous recursive curves (Gottschau et al., 2016). Space-filling regularization is thus always constrained by the combinatorics of the ambient domain.

The current frontier is not a single best criterion but a family of partially compatible ones. The literature already spans discrepancy, repulsion, condition number, energy distance, trajectory occupancy, PDE regularity, GP worst-case uncertainty, and fragmentation. This suggests that future work will continue to be problem-specific: surrogate modeling on irregular designs, latent-state robustness, sliced and constrained experiments, and learned spatial orderings each require different regularized objects, different operators, and different notions of what it means to “fill space.”