GitHub

Kriging Surrogate Tutorial (1D)

Kriging or Gaussian process regression, is a method of interpolation in which the interpolated values are modeled by a Gaussian process.

We are going to use a Kriging surrogate to optimize $f(x)=(6x-2)^2sin(12x-4)$.

First of all, import Surrogates and Plots.

using Surrogates
using Plots

Sampling

We choose to sample f in 100 points between 0 and 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.

f(x) = (6 * x - 2)^2 * sin(12 * x - 4)

n_samples = 100
lower_bound = 0.0
upper_bound = 1.0

xs = lower_bound:0.001:upper_bound

x = sample(n_samples, lower_bound, upper_bound, SobolSample())
y = f.(x)

scatter(
    x, y, label = "Sampled points", xlims = (lower_bound, upper_bound), ylims = (-7, 17))
plot!(xs, f.(xs), label = "True function", legend = :top)
Example block output

Building a surrogate

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

kriging_surrogate behaves like an ordinary function, which we can simply plot. A nice statistical property of this surrogate is being able to calculate the error of the function at each point. We plot this as a confidence interval using the ribbon argument.

kriging_surrogate = Kriging(x, y, lower_bound, upper_bound)

plot(x, y, seriestype = :scatter, label = "Sampled points",
    xlims = (lower_bound, upper_bound), ylims = (-7, 17), legend = :top)
plot!(xs, f.(xs), label = "True function", legend = :top)
plot!(xs, kriging_surrogate.(xs), label = "Surrogate function",
    ribbon = p -> std_error_at_point(kriging_surrogate, p), legend = :top)
Example block output

Optimizing

Having built a surrogate, we can now use it to search for minima in our original function f.

To optimize using our surrogate, we call surrogate_optimize! method. We choose to use Stochastic RBF as the optimization technique and again Sobol sampling as the sampling technique.

surrogate_optimize!(
    f, SRBF(), lower_bound, upper_bound, kriging_surrogate, SobolSample())

scatter(x, y, label = "Sampled points", ylims = (-7, 7), legend = :top)
plot!(xs, f.(xs), label = "True function", legend = :top)
plot!(xs, kriging_surrogate.(xs), label = "Surrogate function",
    ribbon = p -> std_error_at_point(kriging_surrogate, p), legend = :top)
Example block output

Kriging Surrogate Tutorial (ND)

First of all, let's define the function we are going to build a surrogate for. Notice how its argument is a vector of numbers, one for each coordinate, and its output is a scalar.

using Plots
default(c = :matter, legend = false, xlabel = "x", ylabel = "y")
using Surrogates

function branin(x)
    x1 = x[1]
    x2 = x[2]
    a = 1
    b = 5.1 / (4 * π^2)
    c = 5 / π
    r = 6
    s = 10
    t = 1 / (8π)
    a * (x2 - b * x1 + c * x1 - r)^2 + s * (1 - t) * cos(x1) + s
end
branin (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 -5, 10, and 0, 15 for the second dimension. We are taking 50 samples of the space using Sobol sequences. We then evaluate our function on all the sampling points.

n_samples = 100
lower_bound = [-5.0, 0.0]
upper_bound = [10.0, 15.0]

xys = sample(n_samples, lower_bound, upper_bound, GoldenSample())
zs = branin.(xys)
100-element Vector{Float64}:
 158.7018168634276
   1.5528571097928943
  24.772675028600002
  85.5272702554848
 288.4548510839307
  21.60277793374007
  77.23290576102468
  23.78934303608679
  54.16023068493325
  85.1116736050258
   ⋮
  16.94176534139285
  43.474352401085284
 262.4648490860312
  37.018413538543754
  62.492158802600215
  27.698672660634816
 433.59782852353635
  99.7930248958614
 151.05907672958017
x, y = -5:10, 0:15
p1 = surface(x, y, (x1, x2) -> branin((x1, x2)))
xs = [xy[1] for xy in xys]
ys = [xy[2] for xy in xys]
scatter!(xs, ys, zs)
p2 = contour(x, y, (x1, x2) -> branin((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.

kriging_surrogate = Kriging(
    xys, zs, lower_bound, upper_bound, p = [2.0, 2.0], theta = [0.03, 0.003])
(::Kriging{Vector{Tuple{Float64, Float64}}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Vector{Float64}, Float64, Vector{Float64}, Float64, Matrix{Float64}}) (generic function with 2 methods)
p1 = surface(x, y, (x, y) -> kriging_surrogate([x y]))
scatter!(xs, ys, zs, marker_z = zs)
p2 = contour(x, y, (x, y) -> kriging_surrogate([x y]))
scatter!(xs, ys, marker_z = zs)
plot(p1, p2, title = "Surrogate")
Example block output

Optimizing

With our surrogate, we can now search for the minima of the branin function.

Notice how the new sampled points, which were created during the optimization process, are appended to the xys array. This is why its size changes.

size(xys)
(100,)
surrogate_optimize!(branin, SRBF(), lower_bound, upper_bound, kriging_surrogate,
    SobolSample(); maxiters = 100, num_new_samples = 10)
((-3.1474976188433437, 10.610450147441451), 0.39811343961242684)
size(xys)
(132,)
p1 = surface(x, y, (x, y) -> kriging_surrogate([x y]))
xs = [xy[1] for xy in xys]
ys = [xy[2] for xy in xys]
zs = branin.(xys)
scatter!(xs, ys, zs, marker_z = zs)
p2 = contour(x, y, (x, y) -> kriging_surrogate([x y]))
scatter!(xs, ys, marker_z = zs)
plot(p1, p2)
Example block output