GitHub

Gradient Enhanced Kriging Surrogate Tutorial

Gradient-enhanced Kriging is an extension of kriging which supports gradient information. GEK is usually more accurate than kriging. However, it is not computationally efficient when the number of inputs, the number of sampling points, or both, are high. This is mainly due to the size of the corresponding correlation matrix, which increases proportionally with both the number of inputs and the number of sampling points.

Let's have a look at the following function to use Gradient Enhanced Surrogate: $f(x) = x^3 - 6x^2 + 4x + 12$

First of all, we will import Surrogates and Plots packages:

using Surrogates
using Plots

Sampling

We choose to sample f in 100 points between 2 and 10 using the sample function. The sampling points are chosen using a Sobol sequence, this can be done by passing SobolSample() to the sample function.

n_samples = 100
lower_bound = 2
upper_bound = 10
xs = lower_bound:0.001:upper_bound
x = sample(n_samples, lower_bound, upper_bound, SobolSample())
f(x) = x^3 - 6x^2 + 4x + 12
y1 = f.(x)
der = x -> 3 * x^2 - 12 * x + 4
y2 = der.(x)
y = vcat(y1, y2)
scatter(x, y1, label = "Sampled points", xlims = (lower_bound, upper_bound), legend = :top)
plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top)
Example block output

Building a surrogate

With our sampled points, we can build the Gradient Enhanced Kriging surrogate using the GEK function.

my_gek = GEK(x, y, lower_bound, upper_bound, p = 0.03, theta = 0.3)

scatter(x, y1, label = "Sampled points", xlims = (lower_bound, upper_bound), legend = :top)
plot!(f, label = "True function", xlims = (lower_bound, upper_bound), legend = :top)
plot!(my_gek, label = "Surrogate function", ribbon = p -> std_error_at_point(my_gek, p),
    xlims = (lower_bound, upper_bound), legend = :top)
Example block output

Gradient Enhanced Kriging Surrogate Tutorial (ND)

First of all, let's define the function we are going to build a surrogate for.

using Plots
using Surrogates

Now, let's define the function:

function leon(x)
    x1 = x[1]
    x2 = x[2]
    term1 = (x2 - x1^3)^2
    term2 = (1 - x1)^2
    y = term1 + term2
end
leon (generic function with 1 method)

Sampling

Let's define our bounds, this time we are working in two dimensions. In particular, we want our first dimension x to have bounds 0, 1, and 0, 1 for the second dimension. We are taking 100 samples of the space using Sobol Sequences. We then evaluate our function on all the sampling points.

n_samples = 100
lower_bound = [0, 0]
upper_bound = [1, 1]
xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
y1 = leon.(xys)
100-element Vector{Float64}:
 1.1124164985287734
 0.7971708308525649
 0.15007553258897133
 0.9222684086473691
 0.36346270954368265
 0.05102696733797529
 0.12359681268412714
 1.3183157051082617
 0.6632428030300161
 0.20718450295885305
 ⋮
 0.7813352313541237
 0.5076610619732129
 0.05600952949467697
 0.25017019456678113
 1.6610263616386511
 0.890628840729395
 0.3485896751477675
 0.08169559776627366
 1.097441250041019
x, y = 0:1, 0:1
p1 = surface(x, y, (x1, x2) -> leon((x1, x2)))
xs = [xy[1] for xy in xys]
ys = [xy[2] for xy in xys]
scatter!(xs, ys, y1)
p2 = contour(x, y, (x1, x2) -> leon((x1, x2)))
scatter!(xs, ys)
plot(p1, p2, title = "True function")
Example block output

Building a surrogate

Using the sampled points, we build the surrogate, the steps are analogous to the 1-dimensional case.

grad1 = x -> 2 * (x[2] - x[1]^3) * (-3x[1]^2) - 2 * (1 - x[1])
grad2 = x -> 2 * (x[2] - x[1]^3)
d = 2
n = 100
function create_grads(n, d, grad1, grad2, y1)
    c = 0
    y2 = zeros(eltype(y1[1]), n * d)
    for i in 1:n
        y2[i + c] = grad1(xys[i])
        y2[i + c + 1] = grad2(xys[i])
        c = c + 1
    end
    return y2
end
y2 = create_grads(n, d, grad1, grad2, y1)
y = vcat(y1, y2)
300-element Vector{Float64}:
  1.1124164985287734
  0.7971708308525649
  0.15007553258897133
  0.9222684086473691
  0.36346270954368265
  0.05102696733797529
  0.12359681268412714
  1.3183157051082617
  0.6632428030300161
  0.20718450295885305
  ⋮
  1.726403832435608
 -1.8841762827123603
  0.1324101686477661
 -1.6169057588904252
  0.7842198610305786
  0.44991703246523684
 -0.42454326152801514
 -1.832499806454507
  1.5740376710891724
my_GEK = GEK(xys, y, lower_bound, upper_bound)
(::GEK{Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Float64, Vector{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)
p1 = surface(x, y, (x, y) -> my_GEK([x y]))
scatter!(xs, ys, y1, marker_z = y1)
p2 = contour(x, y, (x, y) -> my_GEK([x y]))
scatter!(xs, ys, marker_z = y1)
plot(p1, p2, title = "Surrogate")
Example block output