Lobachevsky surrogate tutorial

Lobachevsky splines function is a function that used for univariate and multivariate scattered interpolation. Introduced by Lobachevsky in 1842 to investigate errors in astronomical measurements.

We are going to use a Lobachevsky surrogate to optimize $f(x)=sin(x)+sin(10/3 * x)$.

First of all import Surrogates and Plots.

using Surrogates
using Plots
default()

Sampling

We choose to sample f in 4 points between 0 and 4 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) = sin(x) + sin(10/3 * x)
n_samples = 5
lower_bound = 1.0
upper_bound = 4.0
x = sample(n_samples, lower_bound, upper_bound, SobolSample())
y = f.(x)
scatter(x, y, label="Sampled points", xlims=(lower_bound, upper_bound))
plot!(f, label="True function", xlims=(lower_bound, upper_bound))

Building a surrogate

With our sampled points we can build the Lobachevsky surrogate using the LobachevskySurrogate function.

lobachevsky_surrogate behaves like an ordinary function which we can simply plot. Alpha is the shape parameters and n specify how close you want lobachevsky function to radial basis function.

alpha = 2.0
n = 6
lobachevsky_surrogate = LobacheskySurrogate(x, y, lower_bound, upper_bound, alpha = 2.0, n = 6)
plot(x, y, seriestype=:scatter, label="Sampled points", xlims=(lower_bound, upper_bound))
plot!(f, label="True function",  xlims=(lower_bound, upper_bound))
plot!(lobachevsky_surrogate, label="Surrogate function",  xlims=(lower_bound, upper_bound))

Optimizing

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

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

@show surrogate_optimize(f, SRBF(), lower_bound, upper_bound, lobachevsky_surrogate, SobolSample())
scatter(x, y, label="Sampled points")
plot!(f, label="True function",  xlims=(lower_bound, upper_bound))
plot!(lobachevsky_surrogate, label="Surrogate function",  xlims=(lower_bound, upper_bound))

In the example below, it shows how to use lobachevsky_surrogate for higher dimension problems.

Lobachevsky Surrogate Tutorial (ND):

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

function schaffer(x)
    x1=x[1]
    x2=x[2]
    fact1 = x1 ^2;
    fact2 = x2 ^2;
    y = fact1 + fact2;
end

schaffer (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, 8, and 0, 8 for the second dimension. We are taking 60 samples of the space using Sobol Sequences. We then evaluate our function on all of the sampling points.

n_samples = 60
lower_bound = [0.0, 0.0]
upper_bound = [8.0, 8.0]

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

60-element Array{Float64,1}:
  56.65625
  40.65625
  22.65625
  12.65625
  80.65625
  38.65625
  52.65625
   3.90625
  55.90625
  68.90625
   ⋮
  31.0078125
 126.0078125
  39.0078125
  39.0078125
   8.0078125
  71.0078125
  55.5078125
  31.5078125
   8.0078125

Building a surrogate

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

Lobachesky = LobacheskySurrogate(xys, zs,  lower_bound, upper_bound, alpha = [2.4,2.4], n=8)

(::LobacheskySurrogate{Array{Tuple{Float64,Float64},1},Array{Float64,1},Array{Float64,1},Int64,Array{Float64,1},Array{Float64,1},Array{Float64,1},Bool}) (generic function with 2 methods)

Optimizing

With our surrogate we can now search for the minimas of the 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(Lobachesky.x)

(60,)

surrogate_optimize(schaffer, SRBF(), lower_bound, upper_bound, Lobachesky, SobolSample(), maxiters=1, num_new_samples=10)

((0.9375, 0.9375), 1.7578125)

size(Lobachesky.x)

(60,)