Kriging Metamodeling & Uncertainty Calibration

Overview

In simulation-based reliability assessment of Non-Destructive Evaluation (NDE) inspections, Gaussian Process (Kriging) metamodeling provides a powerful non-parametric surrogate for mapping complex physics across multi-dimensional input domains, as formalized in Malkiel et al. (2026).

digiqual includes a complete Kriging suite featuring:

Anisotropic Covariance Function Selection (Rank 2): Automated evaluation and selection across Matérn 3/2, Matérn 5/2, Gaussian (RBF), and Rational Quadratic kernels with per-dimension anisotropic length scale vectors \(\boldsymbol{\theta} = [\theta_1, \dots, \theta_d]\).
Matrix-Based Leave-One-Out Cross-Validation (\(MSE_{\text{LOO}}\)): Fast \(O(N^3)\) matrix inverse \(B = S^{-1}\) for exact LOO predictions \(\mu_{\hat{Y}_{-i}}\) and variances \(\sigma_{\hat{Y}_{-i}}^2\).
Standardized LOO Residual Outlier Calibration (Rank 4): Automated detection of extreme observation errors \(|e_i| > 3\) and prior variance scaling (\(\hat{\sigma}^2 \leftarrow \gamma \cdot \hat{\sigma}^2\)) to guarantee conservative 95% worst-case bounds.

1. Covariance Kernel Selection & Anisotropic Fitting

Unlike classical isotropic models that assume equal sensitivity across all directions, digiqual optimizes anisotropic length scale parameters \(\boldsymbol{\theta} = [\theta_1, \theta_2, \dots, \theta_d]\) for each input dimension.

Candidate Kernel Families

During model training (fit_all_robust_mean_models), digiqual evaluates four physical covariance functions:

Matérn 3/2: \[ R(\mathbf{x}, \mathbf{x}') = \left(1 + \sqrt{3} d\right) \exp\left(-\sqrt{3} d\right) \]
Matérn 5/2: \[R(\mathbf{x}, \mathbf{x}') = \left(1 + \sqrt{5} d + \frac{5}{3} d^2\right) \exp\left(-\sqrt{5} d\right)\]
Gaussian (RBF): \[ R(\mathbf{x}, \mathbf{x}') = \exp\left(-\frac{1}{2} d^2\right) \]
Rational Quadratic: \[ R(\mathbf{x}, \mathbf{x}') = \left(1 + \frac{d^2}{2 \alpha \ell^2}\right)^{-\alpha} \]

where the anisotropic distance \(d\) is given by: \[ d = \sqrt{\sum_{j=1}^d \left( \frac{x_j - x_j'}{\theta_j} \right)^2 } \]

Selection Metric (\(MSE_{\text{LOO}}\))

Each kernel family is evaluated via Leave-One-Out Cross-Validation. The kernel candidate achieving the lowest normalized LOO Mean Squared Error is automatically selected:

\[ MSE_{\text{LOO}} = \frac{1}{m} \sum_{i=1}^m \frac{\left(y_i - \mu_{\hat{Y}_{-i}}\right)^2}{\text{var}(y)} \]

2. LOO Residuals & Outlier Calibration

To ensure that Kriging surrogate uncertainty bounds are mathematically valid for safety-critical inspection qualification, digiqual validates prediction errors using standardized LOO residuals.

Augmented Covariance Inversion

Construct the augmented covariance matrix \(S\): \[ S = \begin{bmatrix} \hat{\sigma}^2 R + \alpha I & \mathbf{1} \\ \mathbf{1}^T & 0 \end{bmatrix} \]

Inverting \(S\) yields matrix \(B = S^{-1}\). The Leave-One-Out prediction mean \(\mu_{\hat{Y}_{-i}}\) and standard deviation \(\sigma_{\hat{Y}_{-i}}\) for training observation \(i\) are calculated analytically:

\[ \mu_{\hat{Y}_{-i}} = -\sum_{j \neq i} \frac{B_{ij}}{B_{ii}} y_j, \qquad \sigma_{\hat{Y}_{-i}} = \sqrt{\frac{1}{B_{ii}}} \]

Standardized Residuals (\(e_i\))

The standardized LOO residual vector measures how many standard deviations each prediction deviates from the observed signal:

\[ e_i = \frac{y_i - \mu_{\hat{Y}_{-i}}}{\sigma_{\hat{Y}_{-i}}} \]

Variance Inflation Factor (\(\gamma\))

Standardized residuals are expected to follow a Standard Normal distribution \(\mathcal{N}(0, 1)\). Residuals outside \([-3, 3]\) indicate localized surrogate mismatch or extreme scatter.

If \(\max |e_i| > 3\), digiqual calculates a conservative variance scaling factor: \[ \gamma = \max\left(1.0, \frac{\max |e_i|}{3.0}\right) \]

During noise estimation in fit_variance_model, raw residuals are scaled by \(\sqrt{\gamma}\), inflating local prediction variance \(\hat{\sigma}^2(x)\) by factor \(\gamma\). This guarantees mathematically reliable 95% worst-case confidence bounds even in the presence of extreme outliers.

3. Python Usage Example

Below is a complete Python script demonstrating Kriging surrogate model fitting, kernel optimization, residual computation, diagnostic plotting, and PoD generation:

import numpy as np
import matplotlib.pyplot as plt
from digiqual.pod import (
    fit_all_robust_mean_models,
    compute_kriging_loo_residuals,
    fit_variance_model,
    infer_best_distribution,
    compute_pod_curve,
)
from digiqual.plotting import plot_kriging_diagnostics

# 1. Generate Synthetic 2D NDE Data (Flaw Length & Probe Angle)
np.random.seed(42)
N = 100
length = np.random.uniform(0.5, 5.0, N)
angle = np.random.uniform(-15.0, 15.0, N)
X = np.column_stack([length, angle])

# Signal response with non-smooth physical behavior
y = 2.0 * length - 0.3 * (length**2) + 0.05 * angle + np.random.normal(0, 0.3, N)
# Inject an outlier to test calibration
y[10] += 3.5

# 2. Fit Surrogate Pool & Evaluate Anisotropic Kernels
models, cv_scores, winner_key = fit_all_robust_mean_models(X, y)
gpr = models[("Kriging", None)]

print(f"Selected Kriging Kernel: {gpr.best_kernel_name_}")
print(f"Outlier Calibration Factor (gamma): {gpr.outlier_scale_factor_:.3f}")

# 3. Compute LOO Residuals
loo_means, loo_stds, std_res, gamma = compute_kriging_loo_residuals(gpr, X, y)

# 4. Generate Diagnostic Plot
ax = plot_kriging_diagnostics(
    std_residuals=std_res,
    outlier_scale_factor=gamma,
    best_kernel_name=gpr.best_kernel_name_,
)
plt.show()

Selected Kriging Kernel: RBF (Gaussian)
Outlier Calibration Factor (gamma): 1.912

Summary of API Functions

Function Name	Description
`pod.fit_all_robust_mean_models`	Fits polynomial pool and evaluates anisotropic Kriging kernel candidates via \(MSE_{\text{LOO}}\).
`pod.compute_kriging_loo_residuals`	Inverts augmented matrix \(S\) to compute exact LOO predictions, standardized residuals \(e_i\), and scale factor \(\gamma\).
`plotting.plot_kriging_diagnostics`	Generates dual-panel diagnostic plot visualizing standardized residual distribution vs \(\mathcal{N}(0, 1)\) and \([-3, 3]\) bounds.