surrogate-modeling

name	surrogate-modeling
description	Surrogate modeling (metamodeling) — Gaussian process/kriging (DACE), polynomial response surface (RSM), radial basis functions (RBF), support vector regression, neural network surrogates, design of experiments (LHS, CCD, Sobol), cross-validation (LOOCV, k-fold), PRESS statistic, RMSE/R²/Q², surrogate-based optimization (EGO, expected improvement), multi-fidelity surrogates, and variance-based sensitivity (Sobol indices).
metadata	{"priority":7,"promptSignals":{"phrases":["surrogate model","metamodel","Gaussian process regression","kriging model","response surface method","surrogate-based optimization"],"minScore":3}}

Surrogate Modeling (Metamodeling) — Complete Skill

Purpose and Framework

Surrogate model (metamodel): cheap-to-evaluate mathematical approximation of expensive simulation/experiment Replace: FEA (hours per run), CFD (days per run), physical test Enable: optimization, uncertainty quantification, sensitivity analysis, real-time prediction

Workflow:

1. Design of experiments (DoE) → sampling plan x₁...x_N in design space X ⊆ ℝᵈ
2. Run expensive model at each sample: y_i = f(x_i); collect training data {X, y}
3. Fit surrogate ŷ(x) to training data
4. Validate: cross-validation metrics (RMSE, Q², PRESS)
5. Use surrogate: optimization, uncertainty propagation, sensitivity

Design of Experiments (DoE)

Space-Filling Designs

Latin Hypercube Sampling (LHS): Divide each dimension [0,1] into N equal intervals; place exactly one sample per interval in each dimension; random permutation LHS ensures projection: no two samples share row/column in any 2D projection Correlation-minimized LHS (maximin): maximize minimum distance between any two points N_recommended = 10d to 20d (d = input dimension) for kriging; minimum N = (d+1)(d+2)/2 for quadratic RSM

Sobol quasi-random sequences: Low-discrepancy sequences fill space more uniformly than random; van der Corput base-2 sequences Better space coverage than random for N < 1,000; widely used for UQ and sensitivity

Full factorial: N = L^d (L levels, d factors); impractical for large d Central composite design (CCD): 2^d factorial + axial points at ±α = 2^(d/4) + center points; d+1 ≤ N ≤ 2d+2d+1

Optimal designs (D-optimal, I-optimal): D-optimal: minimize det(X^T X)^(-1) → minimizes parameter estimation variance I-optimal: minimize average prediction variance over design space Computed by coordinate exchange algorithm (JMP, Design-Expert)

Polynomial Response Surface (RSM)

Regression Model

Linear RSM: ŷ(x) = β₀ + Σᵢ βᵢxᵢ [p = d+1 coefficients; n ≥ p samples]

Quadratic RSM (most common): ŷ(x) = β₀ + Σᵢ βᵢxᵢ + Σᵢ βᵢᵢxᵢ² + Σᵢ<ⱼ βᵢⱼxᵢxⱼ [p = (d+1)(d+2)/2 coefficients]

Least squares fit: β = (X^T X)^(-1) X^T y [matrix X contains polynomial basis functions evaluated at design points] Variance of β: Cov(β) = σ² (X^T X)^(-1) [σ² = MSE from residuals]

Prediction variance: var(ŷ(x)) = σ² × x^T(X^T X)^(-1)x [increases away from training points]

Applicability: Best for: smooth, well-behaved responses; global approximation; simple interpretation Weakness: poor for highly nonlinear responses; polynomial order must be specified

Kriging (Gaussian Process Regression)

DACE Model

Kriging predictor: ŷ(x) = β₀ + Σᵢ αᵢ × R(x, xᵢ) [β₀ = global trend; R = correlation function; αᵢ = weights]

Or in Gaussian process notation: ŷ(x) = μ + r^T(x) R⁻¹ (y - μ1) [universal kriging; R = N×N correlation matrix between training points]

Prediction mean and variance: μ̂(x) = μ + r^T R⁻¹ (y - μ1) σ̂²(x) = σ²_process × (1 - r^T R⁻¹ r + (1 - 1^T R⁻¹ r)² / (1^T R⁻¹ 1))

σ̂²(x) = 0 at training points (interpolation); increases away from training points (UQ)

Correlation Functions

Gaussian (squared exponential) — most common: R(xᵢ, xⱼ) = exp(−Σₖ θₖ × |xᵢₖ − xⱼₖ|²) [smooth, infinitely differentiable] θₖ: length-scale parameters (larger θₖ → faster correlation decay in dimension k)

Matérn 5/2: R(r) = (1 + √5 r/ℓ + 5r²/(3ℓ²)) × exp(−√5 r/ℓ) [less smooth; better for physical responses with kinks] Preferred over Gaussian in many engineering applications (Matérn ν=5/2 → twice differentiable)

Power exponential: R(xᵢ, xⱼ) = exp(−Σₖ θₖ × |xᵢₖ − xⱼₖ|^pₖ) [p = 2 → Gaussian; p → 0 → Ornstein-Uhlenbeck]

Kriging hyperparameter estimation: Maximum likelihood estimation (MLE): maximize L(θ, σ², μ) w.r.t. θ, σ², μ Θ* = argmax [−N/2 × log(σ̂²) − 1/2 × log(det(R))] Solved by numerical optimization (Nelder-Mead, L-BFGS-B); start from multiple initial θ values

Software: DACE toolbox (Matlab); scikit-learn GaussianProcessRegressor; SMT (Python Surrogate Modeling Toolbox); PyKrige

Radial Basis Functions (RBF)

RBF approximation: ŷ(x) = Σᵢ wᵢ × φ(||x − xᵢ||) [wᵢ = weights; φ = radial basis function]

Common basis functions:

Type	φ(r)	Properties
Gaussian	exp(−(εr)²)	smooth; shape parameter ε
Multiquadric	√(1 + (εr)²)	globally supported; Hardy
Inverse multiquadric	1/√(1 + (εr)²)	interpolating
Thin plate spline	r² ln r	no shape parameter needed
Cubic	r³	polynomial precision 1

Weight solution: Φ w = y [Φᵢⱼ = φ(||xᵢ − xⱼ||); matrix equation; typically augmented with polynomial for well-posedness]

Shape parameter ε selection: cross-validation or Rippa's algorithm (leave-one-out analytically)

Neural Network Surrogates

Feedforward network for regression: ŷ = W_L × σ(W_{L-1} × ... × σ(W₁ × x + b₁) ... + b_{L-1}) + b_L

Typical architecture: d → 32 → 32 → 1 (2 hidden layers; tanh or ReLU activation) Training: Adam optimizer; batch size 32–128; early stopping on validation loss Regularization: L2 weight decay (λ = 1e-4 to 1e-3); dropout (p = 0.1–0.3)

Bayesian neural network (BNN): Variational inference: approximate posterior q(w|θ) ≈ p(w|D); predictions with uncertainty More expensive to train; MC dropout as cheap approximation (Gal & Ghahramani)

When to use NN vs. kriging: d < 20: kriging usually better (fewer hyperparameters; principled UQ) d > 50: NN often better (curse of dimensionality less severe with NN depth) N > 10,000: NN required (sparse GP or inducing points for large-N kriging)

Validation Metrics

Cross-Validation

Leave-One-Out Cross-Validation (LOOCV): Fit model on N-1 points; predict left-out point; repeat N times PRESS (Predicted Residual Sum of Squares): PRESS = Σᵢ (yᵢ − ŷ₋ᵢ(xᵢ))² [ŷ₋ᵢ = prediction without point i in training]

LOOCV analytically (ordinary least squares): ŷ₋ᵢ = ŷ(xᵢ) − hᵢᵢ × eᵢ / (1 − hᵢᵢ) [hᵢᵢ = i-th diagonal of hat matrix H = X(X^TX)^{-1}X^T] PRESS computed without refitting N times — computationally efficient for RSM

k-fold CV: Split N samples into k folds; train on k-1, test on 1; rotate; k = 5 or 10 typical More stable than LOOCV for small N

Goodness-of-Fit Metrics

Root Mean Square Error (RMSE): RMSE = √(PRESS/N_test) or √(Σ(yᵢ−ŷᵢ)²/N_test) [same units as y; lower is better]

Q² (predictive R²) — preferred for validation: Q² = 1 − PRESS / Σ(yᵢ − ȳ)² [based on PRESS; analogous to R² but for prediction] Q² > 0.9: excellent; Q² > 0.8: acceptable; Q² < 0.7: poor or overfitted

R² (coefficient of determination — training only): R² = 1 − RSS/TSS = 1 − Σ(yᵢ−ŷᵢ)²/Σ(yᵢ−ȳ)² [high R² on training ≠ good predictive model] Always use Q² for surrogate validation, not R²

Maximum absolute error: MAX_ERR = max|yᵢ − ŷᵢ| [important for optimization — worst-case prediction error]

Example validation: Training N = 50 runs; validation N_test = 10 held-out runs; kriging model y_test = [1.23, 2.45, 0.87, ...]; ŷ_test = [1.19, 2.52, 0.91, ...] RMSE = √(Σ(error²)/10) = 0.042; Q² = 0.971 → excellent kriging fit

Surrogate-Based Optimization (SBO)

Efficient Global Optimization (EGO)

Expected Improvement (EI) infill criterion: EI(x) = (y_best − μ̂(x)) × Φ((y_best − μ̂(x))/σ̂(x)) + σ̂(x) × φ((y_best − μ̂(x))/σ̂(x)) [Φ = standard normal CDF; φ = PDF; y_best = current best; μ̂, σ̂ = kriging prediction and uncertainty] Balances exploitation (low μ̂) and exploration (high σ̂)

EGO algorithm (Jones et al. 1998):

1. Build initial kriging from DoE (N₀ = 10d runs)
2. Maximize EI over design space (cheap optimizer — COBYLA, GA)
3. Run expensive model at x* = argmax(EI)
4. Add (x*, y*) to training set; refit kriging
5. Repeat until EI_max < tolerance or budget exhausted

Convergence: EI → 0 as kriging uncertainty collapses; typical budget: 50–200 function evaluations for d ≤ 10

Parallel EGO: q-EI for batch of q parallel evaluations; computed via MC integration KB (Kriging Believer): compute EI assuming pending evaluations return ŷ

Other Infill Criteria

Probability of Improvement: PI(x) = Φ((y_best − μ̂(x))/σ̂(x)) [less explorative than EI; converges faster near optimum]

Lower Confidence Bound (LCB): LCB(x) = μ̂(x) − κ × σ̂(x) [κ = 1–3; larger κ → more exploration; used in Bayesian optimization]

MAXIMIN: maximize σ̂(x) only → pure exploration; for space filling

Variance-Based Sensitivity (Sobol Indices)

Sobol decomposition: V(Y) = Σᵢ Vᵢ + Σᵢ<ⱼ Vᵢⱼ + ... + V₁₂...ₐ [ANOVA-like decomposition; V(Y) = total output variance]

First-order Sobol index: Sᵢ = Vᵢ/V(Y) = V[E(Y|Xᵢ)]/V(Y) [fraction of variance due to Xᵢ alone]

Total Sobol index: STᵢ = 1 − V[E(Y|X_{~i})]/V(Y) [total effect including interactions; STᵢ ≥ Sᵢ always]

Saltelli estimator (Monte Carlo): N_MC = N × (d+2) runs total [N = 1,000–10,000; d = number of inputs] With surrogate: N_samples → ∞ essentially free; most practical approach

Interpretation: Sᵢ ≈ STᵢ: additive model; no interactions STᵢ − Sᵢ large: Xᵢ has important interactions with other inputs STᵢ < 0.05: Xᵢ can be fixed for optimization without significant loss

Multi-Fidelity Surrogates

Co-kriging (Kennedy-O'Hagan): High-fidelity: f_H(x); expensive (N_H runs) Low-fidelity: f_L(x); cheap (N_L >> N_H runs) f_H(x) = ρ × f_L(x) + δ(x) [ρ = scale factor; δ = discrepancy process; both modeled as GPs]

Autoregressive form: Posterior: p(f_H | D_H, D_L) analytically from joint GP

Savings: total N_H can be 5–10× smaller; use low-fidelity (coarse mesh, ROM) to guide high-fidelity sampling

Space mapping: align low-fidelity to high-fidelity by parameter mapping P: x_L → x_H; iterate

Standards and References

Reference	Scope
Jones, Schonlau, Welch (1998) J. Global Opt.	EGO original paper
DACE toolbox (Lophaven, Nielsen, Søndergaard)	Kriging implementation
Sacks et al. (1989) Statistical Science	Design and analysis of computer experiments
SMT: Surrogate Modeling Toolbox (Python)	Open-source; LHS, kriging, RBF, NN
Saltelli (2002) Computer Physics Comm.	Efficient Sobol sensitivity estimators
ISO 21748	Measurement uncertainty guide (relevant for surrogate UQ)

Output

Provide: problem description (output y; input variables x₁...xd with bounds; expensive model: FEA/CFD/test; cost per run; total budget N_budget), DoE (design type: LHS/CCD/Sobol; N samples; rationale for N relative to d; scaled to unit hypercube), surrogate type selection (kriging: smooth response, d < 20; RBF: large dataset, multivariate; NN: d > 50, large N; RSM: linear/quadratic known; reasoning), training (correlation function or basis; hyperparameters estimated by MLE or CV; training time), validation metrics (RMSE = [value]; Q² = [value]; MAX_ERR = [value]; N_validation holdout; verdict: excellent/acceptable/poor), sensitivity analysis (Sobol Sᵢ and STᵢ for each input; rank by importance; which inputs can be fixed), optimization (if SBO: EGO with EI infill; N_init + N_iterations; x_optimal from surrogate; verify with true model at x_optimal; actual vs. predicted y at x_optimal [%error]), multi-fidelity (if applicable: N_H = [value]; N_L = [value]; ρ = [value]; savings vs. single-fidelity), and applicable reference (DACE toolbox; Jones 1998 EGO; SMT toolbox for implementation).