GitHub

GEKPLS Surrogate Tutorial

Gradient Enhanced Kriging with Partial Least Squares Method (GEKPLS) is a surrogate modeling technique that brings down computation time and returns improved accuracy for high-dimensional problems. The Julia implementation of GEKPLS is adapted from the Python version by SMT which is based on this paper.

Surrogates.GEKPLS — Type
GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, theta;
       nugget = 10.0 * eps(), noise = 0.0)
GEKPLS(x, y, grads, lb, ub; n_comp = 2, delta_x = 1.0e-4, extra_points = 2,
       theta = fill(1.0e-2, n_comp), nugget = 10.0 * eps(), noise = 0.0)

Fit a gradient-enhanced Kriging model after projecting the input dimensions with partial least squares (PLS). This model is intended for problems where function values and input gradients are available at every training point.

Fields

  • x: training points in their caller-provided representation.
  • y: training responses.
  • x_matrix: matrix form of the training points.
  • y_matrix: column-matrix form of the responses.
  • grads: gradient matrix associated with x_matrix.
  • xl: two-column matrix of lower and upper bounds.
  • delta: Taylor-expansion step used to create gradient-enhanced points.
  • extra_points: number of additional Taylor points per training point.
  • num_components: number of retained PLS components.
  • beta: generalized least-squares trend coefficients.
  • gamma: correlation residual coefficients.
  • theta: correlation scales in the reduced PLS space.
  • reduced_likelihood_function_value: fitted reduced likelihood value.
  • X_offset: input centering values.
  • X_scale: input scaling values.
  • X_after_std: standardized gradient-enhanced training matrix.
  • pls_mean: mean PLS projection matrix.
  • y_mean: response centering value.
  • y_std: response scaling value.
  • nugget: starting jitter added to the correlation diagonal.
  • noise: observation-noise term added alongside the nugget.

Arguments

  • x: vector of training points.
  • y: scalar response at each point in x.
  • grads: gradient at each point in x, in the same input-coordinate order.
  • n_comp::Integer: number of PLS components to retain.
  • delta_x: first-order Taylor-expansion step.
  • lb: lower bound for each input coordinate.
  • ub: upper bound for each input coordinate.
  • extra_points::Integer: number of gradient-enhanced points used by PLS.
  • theta: correlation scales, one per PLS component. Used as given unless optimize_theta is set, in which case it is the starting point of the fit. The value matters a great deal — on the welded-beam benchmark theta = 1 predicts about thirteen times more accurately than the conventional 0.01.

Keywords

  • nugget: starting jitter added to the correlation diagonal. It is raised by factors of ten only as far as the Cholesky factorization requires, so a well-conditioned problem pays the smallest jitter that works. Oversizing it is not free: the fixed 1e6 * eps() this replaced costs more than a factor of two in RMSE on the welded-beam problems.
  • noise: observation-noise term added alongside the nugget; increasing it lowers the reduced likelihood value.
  • optimize_theta: whether to fit theta by maximizing the reduced likelihood. Defaults to false, unlike Kriging and GEK: every evaluation factorizes an nt × nt matrix with nt = n(1 + extra_points), so the search is far more expensive here than the fit itself. Worth setting on a modest design — on welded beam it lowered RMSE by a factor of 2.75 — and worth avoiding on a large one.
  • n_start: Latin-hypercube starts for that search.

Returns

A callable GEKPLS. Calling it with one point returns a scalar prediction. Training points outside [lb, ub] are rejected with an ArgumentError.

Example

using Surrogates, Zygote

f(x) = sum(abs2, x)
lb = [-1.0, -1.0]
ub = [1.0, 1.0]
x = sample(20, lb, ub, SobolSample())
y = f.(x)
grads = Zygote.gradient.(f, x)
surrogate = GEKPLS(x, y, grads, 1, 1.0e-4, lb, ub, 1, [0.01])
surrogate((0.25, 0.5))
source

The following are the inputs when building a GEKPLS surrogate:

  1. x - The vector containing the training points
  2. y - The vector containing the training outputs associated with each of the training points
  3. grads - The gradients at each of the input X training points
  4. n_comp - Number of components to retain for the partial least squares regression (PLS)
  5. delta_x - The step size to use for the first order Taylor approximation
  6. lb - The lower bound for the training points
  7. ub - The upper bound for the training points
  8. extra_points - The number of additional points to use for the PLS
  9. theta - The correlation scales, one per PLS component

theta must have exactly n_comp entries, one per PLS component. By default it is used as given — unlike KPLS, GEKPLS does not fit it unless asked — and its magnitude matters a great deal. On the welded-beam problem below, theta = 1 predicts about thirteen times more accurately than the conventional 0.01 starting point. reduced_likelihood_function_value ranks candidate scales, so it can be used to choose between them.

The following keyword arguments are also accepted:

  • optimize_theta: fit theta by maximizing the reduced likelihood, taking the supplied value as the starting point. Defaults to false, unlike Kriging and GEK: every evaluation factorizes an nt × nt matrix with nt = n(1 + extra_points), so the search costs far more here than the fit itself. Worth setting on a modest design and worth avoiding on a large one.
  • n_start: Latin-hypercube starts for that search.
  • nugget: the starting jitter added to the correlation diagonal, raised by factors of ten only as far as the Cholesky factorization requires.
  • noise: an observation-noise term added alongside the nugget.

Basic GEKPLS Usage

The following example illustrates how to use GEKPLS:

using Surrogates
using Zygote

function water_flow(x)
    r_w = x[1]
    r = x[2]
    T_u = x[3]
    H_u = x[4]
    T_l = x[5]
    H_l = x[6]
    L = x[7]
    K_w = x[8]
    log_val = log(r / r_w)
    return (2 * pi * T_u * (H_u - H_l)) /
           (log_val * (1 + (2 * L * T_u / (log_val * r_w^2 * K_w)) + T_u / T_l))
end

n = 1000
lb = [0.05, 100, 63070, 990, 63.1, 700, 1120, 9855]
ub = [0.15, 50000, 115600, 1110, 116, 820, 1680, 12045]
x = sample(n, lb, ub, SobolSample())
grads = gradient.(water_flow, x)
y = water_flow.(x)
n_test = 100
x_test = sample(n_test, lb, ub, GoldenSample())
y_true = water_flow.(x_test)
n_comp = 2
delta_x = 0.0001
extra_points = 2
initial_theta = [0.01 for i in 1:n_comp]
g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta)
y_pred = g.(x_test)
rmse = sqrt(sum(((y_pred - y_true) .^ 2) / n_test))
0.02339530805017101

Using GEKPLS With Surrogate Optimization

GEKPLS can also be used to find the minimum of a function with the optimization function. This next example demonstrates how this can be accomplished.

using Surrogates
using Zygote

function sphere_function(x)
    return sum(x .^ 2)
end

lb = [-5.0, -5.0, -5.0]
ub = [5.0, 5.0, 5.0]
n_comp = 2
delta_x = 0.0001
extra_points = 2
initial_theta = [0.01 for i in 1:n_comp]
n = 100
x = sample(n, lb, ub, SobolSample())
grads = gradient.(sphere_function, x)
y = sphere_function.(x)
g = GEKPLS(x, y, grads, n_comp, delta_x, lb, ub, extra_points, initial_theta)
x_point,
minima = surrogate_optimize!(sphere_function, SRBF(), lb, ub, g,
    RandomSample(); maxiters = 20,
    num_new_samples = 20, needs_gradient = true)
minima
0.12830177445407065