Let's start with something easy to get our hands dirty. I want to build a surrogate for $f(x) = \log(x) \cdot x^2+x^3$. Let's choose the radial basis surrogate.
using Surrogates f = x -> log(x)*x^2+x^3 lb = 1.0 ub = 10.0 x = sample(50,lb,ub,SobolSample()) y = f.(x) my_radial_basis = RadialBasis(x,y,lb,ub) #I want an approximation at 5.4 approx = my_radial_basis(5.4)
Let's now see an example in 2D.
using Surrogates using LinearAlgebra f = x -> x*x lb = [1.0,2.0] ub = [10.0,8.5] x = sample(50,lb,ub,SobolSample()) y = f.(x) my_radial_basis = RadialBasis(x,y,lb,ub) #I want an approximation at (1.0,1.4) approx = my_radial_basis((1.0,1.4))
Let's now use the Kriging surrogate, which is a single-output Gaussian process. This surrogate has a nice feature: not only does it approximate the solution at a point, it also calculates the standard error at such point. Let's see an example:
using Surrogates f = x -> exp(x)*x^2+x^3 lb = 0.0 ub = 10.0 x = sample(100,lb,ub,UniformSample()) y = f.(x) p = 1.9 my_krig = Kriging(x,y,lb,ub,p=p) #I want an approximation at 5.4 approx = my_krig(5.4) #I want to find the standard error at 5.4 std_err = std_error_at_point(my_krig,5.4)
Let's now optimize the Kriging surrogate using Lower confidence bound method, this is just a one-liner:
Surrogate optimization methods have two purposes: they both sample the space in unknown regions and look for the minima at the same time.
The Lobachevsky surrogate has the nice feature of having a closed formula for its integral, which is something that other surrogates are missing. Let's compare it with QuadGK:
using Surrogates using QuadGK obj = x -> 3*x + log(x) a = 1.0 b = 4.0 x = sample(2000,a,b,SobolSample()) y = obj.(x) alpha = 2.0 n = 6 my_loba = LobachevskySurrogate(x,y,a,b,alpha=alpha,n=n) #1D integral int_1D = lobachevsky_integral(my_loba,a,b) int = quadgk(obj,a,b) int_val_true = int-int @assert int_1D ≈ int_val_true
Basic example of fitting a neural network on a simple function of two variables.
using Surrogates using Flux using Statistics f = x -> x^2 + x^2 bounds = Float32[-1.0, -1.0], Float32[1.0, 1.0] # Flux models are in single precision by default. # Thus, single precision will also be used here for our training samples. x_train = sample(100, bounds..., SobolSample()) y_train = f.(x_train) # Perceptron with one hidden layer of 20 neurons. model = Chain(Dense(2, 20, relu), Dense(20, 1)) loss(x, y) = Flux.mse(model(x), y) # Training of the neural network learning_rate = 0.1 optimizer = Descent(learning_rate) # Simple gradient descent. See Flux documentation for other options. n_epochs = 50 sgt = NeuralSurrogate(x_train, y_train, bounds..., model=model, loss=loss, opt=optimizer, n_echos=n_epochs) # Testing the new model x_test = sample(30, bounds..., SobolSample()) test_error = mean(abs2, sgt(x) - f(x) for x in x_test)