Kriging (Gaussian Process Regression) Surrogate Model

Theoretical Foundations of Kriging

The Surrogate-Model Idea

Gaussian Process Formulation

Treat the response \( f(\mathbf{x}) \) as a Gaussian process with mean \( \mu(\mathbf{x}) \) and covariance kernel \( k(\mathbf{x}, \mathbf{x}') \). Given training points \( X \) and observations \( \mathbf{y} \), the prediction at \( \mathbf{x}_* \) is closed-form:

$$ \hat{f}(\mathbf{x}_*) = \mu + \mathbf{k}_*^T K^{-1} (\mathbf{y} - \mu\mathbf{1}), \qquad \hat{\sigma}^2(\mathbf{x}_*) = k(\mathbf{x}_*, \mathbf{x}_*) - \mathbf{k}_*^T K^{-1} \mathbf{k}_* $$

where \( K \) is the training covariance matrix and \( \mathbf{k}_* \) the train-test covariances. The mean interpolates the training data exactly; the variance is zero at training points and grows away from them.

Kernel Choice and the Smoothness Assumption

Per-dimension length scales \( \ell_i \) (ARD) express "how far you move in that direction before the response changes" — after training, their magnitudes read directly as sensitivity information (large \( \ell_i \) = inert factor).

Training and Validation

Hyperparameter Estimation and Numerical Care

Length scales and variances are set by maximizing the log marginal likelihood — a multimodal optimization, so multi-start (about 10 restarts) is mandatory. The covariance matrix goes ill-conditioned when training points cluster, breaking the Cholesky factorization; the fix is a nugget added to the diagonal. The nugget is not just numerics — it is observation-noise variance, and CAE has real numerical noise (remeshing, convergence truncation), so estimating the nugget as a parameter is the safe default.

Initial Sampling Design

Latin hypercube sampling (maximin) is standard, with the rule of thumb \( n = 10d \) — as a starting point, on the assumption that active learning (below) will add more. Set ranges slightly wider than the domain you'll use downstream, because Kriging extrapolates terribly: keep the usage region covered by the training region.

Validation — LOO Cross-Validation and Q²

Quality is checked with leave-one-out cross-validation (closed-form and fast for Kriging):

$$ Q^2 = 1 - \frac{\sum_i (y_i - \hat{y}_{-i})^2}{\sum_i (y_i - \bar{y})^2} $$

Guidelines: \( Q^2 \ge 0.9 \) generally, ~0.8 acceptable for trend-finding optimization, ≥0.95 for reliability tail probabilities. Also check standardized residuals \( (y_i - \hat{y}_{-i})/\hat{\sigma}_{-i} \) stay within ±3 — otherwise the predictive variance itself is untrustworthy (revisit noise and kernel).

Active Learning — Adding Samples Where It Matters

Kriging's superpower is sequential sampling driven by predictive variance, with the acquisition rule chosen per goal:

• Global accuracy — add at maximum predictive variance (simple; drifts toward corners)
• Optimization — Expected Improvement (EI; Bayesian optimization / EGO)
• Reliability — points where the limit-state classification \( \hat{f}(\mathbf{x}) = 0 \) is most ambiguous (U-function, AK-MCS)

The loop "initial LHS → validate → goal-driven additions → convergence check" beats one-shot fixed sampling by a wide margin at equal budget.

Applying It in Practice

Standard Workflow

1. Screen — above ~10 factors, reduce with Morris first (Kriging struggles in high dimensions)
2. Normalize inputs — each factor to [0,1]; stabilizes length-scale estimation across units
3. Initial DOE — LHS with \( 10d \) points; flag and log failed runs
4. Fit and validate — Matérn 5/2 + ARD + estimated nugget as default; LOO Q² and standardized residuals
5. Sequential additions — goal-driven acquisitions to budget; track Q² over rounds
6. Production use — Monte Carlo, Sobol' indices, optimization on the surrogate; re-verify final candidates with real CAE runs

Failed Runs and Discontinuous Responses

When parts of the design space break the analysis (buckling divergence, meshing failure), silently dropping those points lets the surrogate interpolate smoothly across a region it knows nothing about — and then predict the optimum inside the infeasible zone. Countermeasures: learn feasibility with a separate classifier (GP classification) and compose, or restrict ranges if the physical boundary is known. For discontinuous responses (contact on/off, buckling mode switches), a single global Kriging struggles — use domain splitting (cluster + local models) or drop to Matérn 3/2.

From Scalars to Field Outputs

To surrogate full fields (stress maps, temperature fields), first compress with POD to a few modal coefficients and fit one Kriging per coefficient (POD+Kriging). Choose mode count by ~99% cumulative energy and report reconstruction error separately from surrogate error.

Tools and CAE Integration

Library and Tool Comparison

Minimal Example (scikit-learn)

import numpy as np
from sklearn.gaussian_process import GaussianProcessRegressor
from sklearn.gaussian_process.kernels import Matern, WhiteKernel, ConstantKernel

kernel = (ConstantKernel() * Matern(length_scale=np.ones(d), nu=2.5)
          + WhiteKernel(noise_level=1e-6))          # nugget = numerical noise
gp = GaussianProcessRegressor(kernel=kernel, normalize_y=True,
                              n_restarts_optimizer=10)  # multi-start
gp.fit(X_train, y_train)                            # X_train: CAE results at LHS points
y_pred, y_std = gp.predict(X_new, return_std=True)  # prediction and its std

Make it a habit to print gp.kernel_ after fitting — the length scales double-check both physical plausibility (do influential factors have small scales?) and hyperparameter health (no boundary sticking).

Research Frontiers

Multi-Fidelity Kriging

Co-Kriging / MFK fuses cheap low-fidelity data (coarse meshes) with a few expensive high-fidelity runs — a natural fit for CAE. The low-fidelity model captures the global shape; a handful of high-fidelity samples corrects it, routinely reaching target accuracy at a fraction of the high-fidelity-only cost. Hierarchies beyond mesh density work too: 2-D vs. 3-D, linear vs. nonlinear.

Gradient-Enhanced Kriging (GEK)

Where adjoint solvers provide cheap gradients (CFD adjoints, structural design sensitivities), embedding gradient observations in the covariance structure yields \( 1+d \) pieces of information per run — dramatic sample-efficiency gains in higher dimensions. Gradient noise is handled with its own nugget.

High Dimensions and Deep-Learning Hybrids

Practical Kriging tops out around \( d \approx 20 \); beyond that, active-subspace identification and rotation methods lead. Deep kernel learning (NN feature extraction feeding a GP head) pursues discontinuity-tolerant flexibility while keeping GP variance — with the adoption criterion, as ever, being whether the predictive variance stays calibrated (check the standardized-residual distribution).

Troubleshooting

Symptoms, Causes, and Fixes

How Far to Trust the Surrogate

Related: Morris screening, Sparse PCE with LAR, Bayesian calibration.