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,)