Gaussian Process Surrogate Tutorial

Gaussian Process regression in Surrogates.jl is implemented as a simple wrapper around the AbstractGPs.jl package. AbstractGPs comes with a variety of covariance functions (kernels). See KernelFunctions.jl for examples.

Note

The examples below demonstrate the use of AbstractGPs with out-of-the-box settings without hyperparameter optimization (i.e. without changing parameters like lengthscale, signal variance and noise variance.) Beyond hyperparameter optimization, careful initialization of hyperparameters and priors on the parameters is required for this surrogate to work properly. For more details on how to fit GPs in practice, check out A Practical Guide to Gaussian Processes.

Also see this example to understand hyperparameter optimization with AbstractGPs.

1D Example

In the example below, the 'gp_surrogate' assignment code can be commented / uncommented to see how the different kernels influence the predictions.

using Surrogates
using Plots
default()
using AbstractGPs #required to access different types of kernels

f(x) = (6 * x - 2)^2 * sin(12 * x - 4)
n_samples = 4
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)
#gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(SqExponentialKernel()), Σy=0.05) #example of Squared Exponential Kernel
#gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(MaternKernel()), Σy=0.05) #example of MaternKernel
gp_surrogate = AbstractGPSurrogate(x,y, gp=GP(PolynomialKernel(; c=2.0, degree=5)), Σy=0.25)
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!(0:0.001:1, gp_surrogate.gp_posterior; label="Posterior", ribbon_scale=2)

Optimization Example

This example shows the use of AbstractGP Surrogates to find the minima of a function:

using Surrogates
using Plots
f(x) = (x-2)^2
n_samples = 4
lower_bound = 0.0
upper_bound = 4.0
xs = lower_bound:0.1:upper_bound
x = sample(n_samples, lower_bound, upper_bound, SobolSample())
y = f.(x)
gp_surrogate = AbstractGPSurrogate(x,y)
@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, gp_surrogate, SobolSample())

(2.0012499999999998, 1.5624999999993784e-6)

Plotting the function and the sampled points:

scatter(gp_surrogate.x, gp_surrogate.y, label="Sampled points", ylims=(-1.0, 5.0), legend=:top)
plot!(xs, gp_surrogate.(xs), label="Surrogate function", ribbon=p->std_error_at_point(gp_surrogate, p), legend=:top)
plot!(xs, f.(xs), label="True function", legend=:top)

ND Example

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

hypot_func = z -> 3*hypot(z...)+1
n_samples = 50
lower_bound = [-1.0, -1.0]
upper_bound = [1.0, 1.0]

xys = sample(n_samples, lower_bound, upper_bound, SobolSample())
zs = hypot_func.(xys);

x, y = -2:2, -2:2
p1 = surface(x, y, (x1,x2) -> hypot_func((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) -> hypot_func((x1,x2)))
scatter!(xs, ys)
plot(p1, p2, title="True function")

Now let's see how our surrogate performs:

gp_surrogate = AbstractGPSurrogate(xys, zs)
p1 = surface(x, y, (x, y) -> gp_surrogate([x y]))
scatter!(xs, ys, zs, marker_z=zs)
p2 = contour(x, y, (x, y) -> gp_surrogate([x y]))
scatter!(xs, ys, marker_z=zs)
plot(p1, p2, title="Surrogate")

@show gp_surrogate((0.2,0.2))

1.9682380785540037

@show hypot_func((0.2,0.2))

1.8485281374238571

And this is our log marginal posterior predictive probability:

@show logpdf_surrogate(gp_surrogate)

64.91417186860279