Gradient Enhanced Kriging

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 that increases proportionally with both the number of inputs and the number of sampling points.

Let's have a look to the following function to use Gradient Enhanced Surrogate: $f(x) = sin(x) + 2*x^2$

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

using Surrogates
using Plots
default()

Sampling

We choose to sample f in 8 points between 0 to 1 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 = 10
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)

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 = 1.4);

(::GEK{Vector{Float64}, Vector{Float64}, Int64, Int64, Float64, Float64, Float64, Matrix{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)

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)

Gradient Enhanced Kriging Surrogate Tutorial (ND)

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

Now, let's define the function:

function leon(x)
      x1 = x[1]
      x2 = x[2]
      term1 = 100*(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, 10, and 0, 10 for the second dimension. We are taking 80 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points.

n_samples = 80
lower_bound = [0, 0]
upper_bound = [10, 10]
xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
y1 = leon.(xys);

80-element Vector{Float64}:
   1577.8677217621007
      1.8071765303987595e6
      2.1269689354185887e7
  19491.21875369508
 397142.68190494325
      5.1836512480499245e7
      7.28543611393035e6
   1992.7400718308636
   5295.5550515377545
      1.2406020272666672e7
      ⋮
   4920.825319683173
   5947.404192859795
      1.1333865241816055e7
      7.268698331698339e7
 688180.3567782413
  90821.96402000086
      3.03575746835075e7
      3.53029568501928e6
   3643.980667055549

Building a surrogate

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

grad1 = x1 -> 2*(300*(x[1])^5 - 300*(x[1])^2*x[2] + x[1] -1)
grad2 = x2 -> 200*(x[2] - (x[1])^3)
d = 2
n = 10
function create_grads(n, d, grad1, grad2, y)
      c = 0
      y2 = zeros(eltype(y[1]),n*d)
      for i in 1:n
            y2[i + c] = grad1(x[i])
            y2[i + c + 1] = grad2(x[i])
            c = c + 1
      end
      return y2
end
y2 = create_grads(n, d, grad2, grad2, y)
y = vcat(y1,y2)

100-element Vector{Float64}:
   1577.8677217621007
      1.8071765303987595e6
      2.1269689354185887e7
  19491.21875369508
 397142.68190494325
      5.1836512480499245e7
      7.28543611393035e6
   1992.7400718308636
   5295.5550515377545
      1.2406020272666672e7
      ⋮
    200.0
    200.0
    200.0
    200.0
    200.0
    200.0
    200.0
    200.0
    200.0

my_GEK = GEK(xys, y, lower_bound, upper_bound, p=[1.9, 1.9])

(::GEK{Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Int64}, Vector{Int64}, Vector{Float64}, Vector{Float64}, Float64, Matrix{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)