Surrogate Models

This notebook demonstrates how to train, evaluate, and serialize surrogate models using the Standard Evaluator library. Surrogate models are fast mathematical approximations of expensive evaluators, enabling rapid prediction of response values without re-running the underlying analysis. You will learn how to use the PolynomialModel, SMT-based models (such as RBF), serialize and reconstruct trained models, and configure model options.

Installation

The SMT-based surrogate models require the optional smt dependency. Install it with:

pip install standard-evaluator[smt]

The PolynomialModel is available with the base installation and does not require any extra dependencies.

import numpy as np
import pandas as pd
import standard_evaluator as se

matlabengine not available.

Model Hierarchy

The surrogate model classes follow this inheritance hierarchy:

  • SurrogateModel (base class, inherits from NumpyEvaluator)

    • PolynomialModel — Polynomial response surface model (base dependencies only)

    • AbstractSmtModel — Base class for all SMT-based models (requires smt)

      • RadialBasisFunctionModel (RBF) — Radial basis function interpolation

      • InverseDistanceWeightingModel (IDW) — Inverse distance weighting interpolation

      • GradientEnhancedNeuralNetworksModel (GENN) — Gradient-enhanced neural network

      • LeastSquaresApproximationModel (LS) — Least-squares polynomial approximation

      • RegularizedMinimalEnergyTensorProductBSplines (RMTS) — Regularized tensor-product B-splines

      • SecondOrderPolynomialApproximationModel (SOPA) — Second-order polynomial approximation

All surrogate models share a common API:

  • Training: Pass a pd.DataFrame of sites (inputs + true responses) to the constructor

  • Prediction: Call the model with a DataFrame of new input sites (modifies in-place)

  • Serialization: Use to_dict() and from_dict() for round-trip persistence

  • Options: Configure model behavior via a Pydantic options model

Training a PolynomialModel

We start by training a PolynomialModel on sample data from the Sphere test evaluator. First we generate training sites by sampling within the variable bounds, then we evaluate the true responses.

# Create a test evaluator (Sphere function: f = x0^2 + x1^2)
evaluator = se.evaluators.test.Sphere(num_independent=2, num_dependent=1)

# Get variable bounds from the opt_problem
opt_problem = evaluator.opt_problem
variables = opt_problem.variables

# Generate 20 random training sites within bounds
np.random.seed(42)
n_train = 20
train_data = {}
for var in variables:
    lb, ub = var.bounds
    train_data[var.name] = np.random.uniform(lb, ub, n_train)

train_sites = pd.DataFrame(train_data)

# Evaluate true responses at training sites
evaluator(train_sites)
print(f"Training sites shape: {train_sites.shape}")
train_sites.head()

Training sites shape: (20, 3)
x0 x1 f
0 -2.509198 2.237058 11.300501
1 9.014286 -7.210123 133.243225
2 4.639879 -4.157107 38.810014
3 1.973170 -2.672763 11.037061
4 -6.879627 -0.878600 48.101209

Now we train a degree-2 polynomial model on these sites and predict at new test points.

from standard_evaluator.surrogate_models import PolynomialModel, PolynomialModelOptions

# Train the polynomial model (degree 2 fits a quadratic — perfect for Sphere)
poly_model = PolynomialModel(
    sites=train_sites,
    options=PolynomialModelOptions(degree=2),
    opt_problem=opt_problem,
)

# Generate 5 new test sites for prediction
n_test = 5
test_data = {}
for var in variables:
    lb, ub = var.bounds
    test_data[var.name] = np.random.uniform(lb, ub, n_test)

test_sites = pd.DataFrame(test_data)

# Predict using the polynomial model (modifies DataFrame in place)
poly_model(test_sites)
predicted = test_sites["f"].values.copy()

# Compute true values for comparison
true_sites = test_sites[evaluator.inputs].copy()
evaluator(true_sites)
true_values = true_sites["f"].values

# Display predicted vs true
comparison = pd.DataFrame({
    "x0": test_sites["x0"],
    "x1": test_sites["x1"],
    "Predicted": predicted,
    "True": true_values,
    "Error": np.abs(predicted - true_values),
})
print("Polynomial Model Predictions vs True Values:")
comparison

Polynomial Model Predictions vs True Values:
x0 x1 Predicted True Error
0 -7.559235 3.250446 67.707436 67.707436 7.105427e-14
1 -0.096462 -3.765778 14.190392 14.190392 1.421085e-14
2 -9.312230 0.401360 86.878710 86.878710 9.947598e-14
3 8.186408 0.934206 67.890017 67.890017 8.526513e-14
4 -4.824400 -6.302911 63.001525 63.001525 3.552714e-14

The degree-2 polynomial model fits the Sphere function exactly (since the Sphere function is itself a degree-2 polynomial), resulting in near-zero prediction errors.

Serialization and Reconstruction

Surrogate models can be serialized to a dictionary using to_dict() and reconstructed using SurrogateModel.from_dict(). This enables saving trained models to disk and loading them later without retraining.

from standard_evaluator.surrogate_models import SurrogateModel

# Serialize the trained model
model_dict = poly_model.to_dict()
print(f"Serialized model type: {model_dict['type']}")
print(f"Dictionary keys: {list(model_dict.keys())}")

# Reconstruct from dictionary
reconstructed_model = SurrogateModel.from_dict(model_dict)

# Verify predictions match
verify_sites = pd.DataFrame(test_data)  # fresh copy of test inputs
reconstructed_model(verify_sites)
reconstructed_predictions = verify_sites["f"].values

# Compare original and reconstructed predictions
max_diff = np.max(np.abs(predicted - reconstructed_predictions))
print(f"\nMax difference between original and reconstructed predictions: {max_diff:.2e}")
print("Predictions match:", np.allclose(predicted, reconstructed_predictions))

Serialized model type: PolynomialModel
Dictionary keys: ['type', 'info', 'opt_problem', 'version', 'name']

Max difference between original and reconstructed predictions: 0.00e+00
Predictions match: True

The reconstructed model produces identical predictions to the original, confirming that the serialization round-trip preserves the model state completely.

Training an RBF Model (SMT)

The RadialBasisFunctionModel uses radial basis function interpolation from the SMT library. Let’s train it on a more complex function and assess its accuracy.

from standard_evaluator.surrogate_models import RadialBasisFunctionModel

# Use the Rosenbrock test evaluator (a harder, non-polynomial function)
rosenbrock = se.evaluators.test.Rosenbrock(num_independent=2, num_dependent=1)
rb_opt_problem = rosenbrock.opt_problem
rb_variables = rb_opt_problem.variables

# Generate 30 training sites (Latin hypercube-like random sampling)
np.random.seed(123)
n_rb_train = 30
rb_train_data = {}
for var in rb_variables:
    lb, ub = var.bounds
    rb_train_data[var.name] = np.random.uniform(lb, ub, n_rb_train)

rb_train_sites = pd.DataFrame(rb_train_data)
rosenbrock(rb_train_sites)

# Train the RBF model
rbf_model = RadialBasisFunctionModel(
    sites=rb_train_sites,
    opt_problem=rb_opt_problem,
)

# Generate 15 test sites
n_rb_test = 15
rb_test_data = {}
for var in rb_variables:
    lb, ub = var.bounds
    rb_test_data[var.name] = np.random.uniform(lb, ub, n_rb_test)

rb_test_sites = pd.DataFrame(rb_test_data)

# Predict with RBF model
rbf_model(rb_test_sites)
rbf_predicted = rb_test_sites["f"].values.copy()

# Get true values
rb_true_sites = rb_test_sites[rosenbrock.inputs].copy()
rosenbrock(rb_true_sites)
rbf_true = rb_true_sites["f"].values

# Compute RMSE
rmse = np.sqrt(np.mean((rbf_predicted - rbf_true) ** 2))
# Compute R² score
ss_res = np.sum((rbf_true - rbf_predicted) ** 2)
ss_tot = np.sum((rbf_true - np.mean(rbf_true)) ** 2)
r_squared = 1 - ss_res / ss_tot

print(f"RBF Model Accuracy on Rosenbrock function:")
print(f"  RMSE: {rmse:.4f}")
print(f"  R²:   {r_squared:.4f}")

___________________________________________________________________________
   
                                    RBF
___________________________________________________________________________
   
 Problem size
   
      # training points.        : 30
   
___________________________________________________________________________
   
 Training
   
   Training ...
      Initializing linear solver ...
         Performing LU fact. (30 x 30 mtx) ...
         Performing LU fact. (30 x 30 mtx) - done. Time (sec):  0.0003672
      Initializing linear solver - done. Time (sec):  0.0004151
      Solving linear system (col. 0) ...
         Back solving (30 x 30 mtx) ...
         Back solving (30 x 30 mtx) - done. Time (sec):  0.0001214
      Solving linear system (col. 0) - done. Time (sec):  0.0001593
   Training - done. Time (sec):  0.0010283
___________________________________________________________________________
   
 Evaluation
   
      # eval points. : 15
   
   Predicting ...
   Predicting - done. Time (sec):  0.0000215
   
   Prediction time/pt. (sec) :  0.0000014
   
RBF Model Accuracy on Rosenbrock function:
  RMSE: 1.3163
  R²:   0.9656

The RBF model provides good interpolation accuracy. Since RBF is an interpolation method, it passes exactly through training points and approximates well between them.

Model Options

Each surrogate model class defines its configurable options via a Pydantic model. Options control model behavior such as polynomial degree, regularization parameters, or basis function settings.

  • PolynomialModelOptions: degree (int, default=1), coefficient_ordering (enum)

  • RadialBasisFunctionModelOptions: d0 (float, default=1.0), poly_degree (int, default=-1), reg (float, default=1e-10)

You can pass an options instance to the model constructor to override defaults. Let’s see the effect of changing the polynomial degree.

# Train a degree-1 (linear) polynomial model on the same Sphere data
poly_linear = PolynomialModel(
    sites=train_sites,
    options=PolynomialModelOptions(degree=1),  # Linear fit (non-default)
    opt_problem=opt_problem,
)

# Predict at test sites
linear_test = pd.DataFrame(test_data)
poly_linear(linear_test)
linear_predicted = linear_test["f"].values

# Compare linear vs quadratic fit
options_comparison = pd.DataFrame({
    "x0": test_sites["x0"],
    "x1": test_sites["x1"],
    "True": true_values,
    "Degree 1 (linear)": linear_predicted,
    "Degree 2 (quadratic)": predicted,
})
print("Effect of polynomial degree on predictions:")
options_comparison

Effect of polynomial degree on predictions:
x0 x1 True Degree 1 (linear) Degree 2 (quadratic)
0 -7.559235 3.250446 67.707436 65.141901 67.707436
1 -0.096462 -3.765778 14.190392 72.438339 14.190392
2 -9.312230 0.401360 86.878710 67.555390 86.878710
3 8.186408 0.934206 67.890017 69.076106 67.890017
4 -4.824400 -6.302911 63.001525 74.223847 63.001525

The degree-1 model cannot capture the quadratic Sphere function accurately, while the degree-2 model fits it exactly. This demonstrates how the degree option directly affects model capability and accuracy.