Replication or exploration? Sequential design for stochastic simulation experiments

1 Introduction

• Stochastic simulations can require larger experiments and more sophisticated GP models to separate signal from noise and track input-dependent variance.The paper notes that available tools often add modeling and computational complexity.Source: §1 · pp. 1–2
• The unresolved design problem is how to allocate unique locations and replication to obtain the best overall GP fit, rather than using a uniform replicate count.Prior replication criteria were limited to finite candidate designs.Source: §1 · p. 2
• A sequential IMSPE criterion weighs exploring now against replicating now by accounting for the options each choice creates later.It also exploits fast GP updates and a closed-form IMSPE expression and derivatives for numerical optimization.Source: §1 · p. 2
• The lookahead criterion is not full lookahead, and longer horizons tend to favor more replication.The horizon can be tuned online; adaptive choices may prioritize modeling cost or out-of-sample RMSE.Source: §1 · pp. 2–3

2 Background and proof of concept

• The section’s proof of concept is that replication can benefit both predictive accuracy and computational cost in GP surrogate modeling.Source: §2 · p. 3

2.1 Gaussian process regression with replication

• The background assumes known kernel hyperparameters and noise variance for its initial setup, before addressing the more realistic unknown-parameter case later.Source: §2.1 · p. 3
• Replication reduces GP computation from O(N^3) to O(n^3) without approximation, which matters when noisy, low-signal simulations require many observations.Here N is the total observation count and n is the number of unique input locations.Source: §2.1 · p. 4

2.2 Sequential design for GPs

• One-shot IMSPE design is difficult in an N-by-d search space, so the paper adds one run at a time and optimizes the next location.The sequential IMSPE minimizer can be an existing location, making the next run a replicate.Source: §2.2 · pp. 4–5
• In Figure 1, the blue noise regime favors replicating at x2, whereas the green regime favors exploring near x = 0.32.The figure contrasts the two noise-variance functions and their resulting IMSPE objectives.Source: Figure 1 · p. 5, §2.2 · p. 5
• When IMSPE surfaces are multimodal, a discrete search over replicates can find lower IMSPE than some local minima from continuous optimization.Source: §2.2 · p. 6
• Learning the noise process during sequential sampling permits non-uniform sampling that can adapt to input-dependent noise.The passage reports that assuming constant noise for design can lead to inferior prediction.Source: §2.2 · p. 6

3.1 IMSPE closed-formed expressions

• The unique-n IMSPE formulation supports sequential replicate-versus-new-site comparisons, with O(n^2) calculations and O(n) replicate updates.It also yields an explicit condition for preferring replication: sufficiently high noise can make replication optimal, and it can also be optimal in homoskedastic settings, often at input-space edges depending on the design and replicate counts.Source: §3 · p. 6, §3.1 · pp. 7–9

3 IMSPE through the lens of replication

• The general IMSPE development provides a closed-form expression, including a unique-n representation suited to calculations under replication.The representation uses kernel-product integrals and can be advantageous even when closed forms are unavailable, because numerical integration inside the kernel matrix may be easier than integrating a highly multimodal predictive variance.Source: §3 · p. 6, §3.1 · p. 7

3.2 Gradient expressions

• Closed-form gradients support sequential IMSPE optimization without further matrix decompositions, at a computational cost of O(n^2).The authors present these as the first derivative expressions of this kind of which they are aware.Source: §3 · p. 6, §3.2 · pp. 9–10

3.3 Looking ahead over replication

• Rollout evaluates h+1 paths of explore/replicate choices and selects the first action on the path with the lowest horizon IMSPE.Longer horizons weakly tilt designs toward replication, but exploration can still be preferred; in the h=3 example, the best path replicates three times before exploring.Source: §3.3 · pp. 11–12, Figure 2 · p. 11

4.1 Heteroskedastic modeling

• The heteroskedastic model uses replication and learns a latent log-variance GP jointly with the mean GP, under the premise that little is known about noise.That premise may not fit mechanisms such as Monte Carlo simulations with known convergence rates.Source: §4 · p. 12, §4.1 · p. 13

4.2 Sequential heteroskedastic modeling

• Sequential updates keep GP hyperparameters fixed during lookahead, but after an observation they update all hyperparameters, using warm starts instead of restarting costly O(n^3) optimization.Previous values are reused where available, with the newest latent requiring initialization.Source: §4.2 · pp. 13–14

4.3 Defining the horizon

• The Target heuristic adjusts the horizon to maintain a chosen ratio of unique locations to total design size, increasing it to encourage replication and decreasing it to encourage exploration.The ratio target is arbitrary; the paper chooses it because its magnitude is easy to interpret.Source: §4.3 · p. 15
• Horizons up to 20 add negligible search cost relative to GP updates, while softly favoring replication without ruling out exploration.The empirical work reports that horizon targeting manages computational costs.Source: §4.3 · p. 15
• For accuracy, an adaptive horizon links the target ratio to signal-to-noise considerations and compares estimated replication allocations with replicates already collected.The allocation heuristic reflects noise levels but does not directly account for GP-estimated correlations.Source: §4.3 · p. 16
• The allocation recipe addresses replication only: design locations must be entered a priori, so it is not directly implementable as a sequential design.The authors describe generating a candidate design first as one possible workaround.Source: §4.3 · p. 16
• In the illustrative 1D example, low-noise regions warrant fewer replicates, while noisier regions warrant more; batch allocation is contrasted with greedy sequential allocation.The greedy comparison reflects that sequential design decisions cannot be undone.Source: §4.3 · p. 16

5.1 Illustrative one-dimensional example

• With h = 4 and Adapt, the best RMSE used 60 unique sites (ρ = 0.12), with some points replicated up to 30 times.Increasing the horizon shifted design density in high-variance regions toward replication without harming performance.Source: §5.1 · p. 17, Figure 4 · p. 18

5.2 Synthetic simulation experiment

• Across 1,000 Monte Carlo replicates, sequential hetGP had the best out-of-sample RMSE; random-replication hetGP won proper scores, while sequential hetGP ranked second on score and log-noise RMSE.The sequential approach's value was much diminished with a homoskedastic GP.Source: Figure 5 · p. 19, §5.2 · p. 19

5.3 Susceptible-Infected-Recovered (SIR) epidemic model

• On SIR, RMSEs were barely distinguishable, but h = 4 and Target achieved the best scores; replication helped where signal-to-noise was low.The underlying mean function was very smooth, while the variance surface was more challenging.Source: Figure 6 · p. 20, §5.3 · pp. 20–21
• Longer SIR horizons replaced clusters of singleton sites with replicates, retaining essentially the same statistical efficiency while reducing computational overhead by more than 3×.The computational cost of looking ahead was negligible next to the savings from using fewer unique locations.Source: §5.3 · pp. 21–22

5.4 Inventory management

• On ATO, shorter lookahead performed better; Adapt matched h = −1, whereas Target at ρ = 0.2 explored too little.Adapt took an average of 682 minutes to reach N = 2000, while Target took 183 minutes.Source: §5.4 · pp. 22–23, Figure 9 · p. 23

6 Conclusion and perspectives

• The sequential design chooses whether to explore or replicate at the next sample, a simpler task than computationally intractable fully batched allocation.The procedure remains myopic; a tractable lookahead scheme emphasizes replication, and tuning its horizon trades inference thriftiness against out-of-sample accuracy.Source: §6 · p. 23
• IMSPE is used for its popularity and because it yields closed-form optimization derivatives, horizon updates, and simplifications for considering replicates.Other design criteria target quantities such as entropy or information about unknown model parameters, and may lead to similar strategies.Source: §6 · pp. 23–24
• The implementation and empirical work use hetGP, while reproducing the larger examples is less straightforward because they require more computing and/or R–MATLAB linking.Supplementary code reproduces the smaller examples; code for the other examples is available upon request.Source: §6 · p. 24
• Heteroskedastic processes favor more replication in high-noise than low-noise regions; under homoskedasticity, the inferred replication level is generally global, with possible increases at input-space edges.The results also apply in the homoskedastic setting, though with less pronounced effects.Source: §6 · p. 24
• Across three example sets, lookahead-based IMSPE yields accurate prediction, and its slight replication bias produces more computationally thrifty predictors without compromising accuracy.The examples span low to high noise, low to moderate input dimension, and synthetic to real simulation experiments.Source: §6 · p. 24

A Detailed gradient expressions

• Appendix A provides expressions for the IMSPE gradient discussed in Section 3.2.The expressions are presented as supplementary detail to that discussion.Source: A Detailed gradient expressions · p. 29
• The derivation notes that W_n does not depend on x̃.The appendix also points to Eq. (11) in presenting the expressions.Source: A Detailed gradient expressions · p. 29
• Expressions for c1 and c2 with particular kernels are located in Appendix B.Source: A Detailed gradient expressions · p. 29

B Expressions for common kernels

• The appendix considers Gaussian and Matérn kernels with α = 5/2, 3/2, and 1/2; the last is the exponential kernel.The kernels are stationary, expressed as k(x, x′) = νc(x − x′).Source: B Expressions for common kernels · p. 29
• For separable kernels over D = [0, 1]^d, the appendix uses separability to obtain corresponding expressions.Source: B Expressions for common kernels · p. 29
• The appendix provides univariate parameterizations and corresponding expressions for w, d, c1, and c2.Source: B Expressions for common kernels · p. 29

B.1 Gaussian kernel

• The Gaussian-kernel expressions include the error function, erf.Source: §B.1 · p. 30
• The subsection gives component expressions for the vector c1 indexed by 1 ≤ i ≤ n.Source: §B.1 · p. 30
• The kernel used in Figure 1 has hyperparameters ν = 1 and θ = 0.01.Source: §B.1 · p. 30

• It gives a Matérn-kernel parameterization for specified α values and expressions for c2 and derivatives with respect to x.Source: §B.2 · pp. 30, 34
• Closed-form derivatives of w(x_i, x_j) are derived by assuming x_i ≤ x_j, removing absolute values, and dividing the integral into three regions.The regions are 0 → x_i, x_i → x_j, and x_j → 1; symbolic solvers are used for the most tedious components.Source: §B.2 · pp. 30–31
• The appendix provides derivatives with respect to x_i and x_j for both orderings, obtaining the x_i > x_j case by swapping the inputs.Source: §B.2 · pp. 32–34