JuliaGaussianProcesses / AbstractGPs.jl

Abstract types and methods for Gaussian Processes.
https://juliagaussianprocesses.github.io/AbstractGPs.jl/dev
Other
218 stars 20 forks source link

AbstractGP shares array data - impossible to add new points to grid (bayesian optimisation) #385

Closed filipporemonato closed 11 months ago

filipporemonato commented 11 months ago

I'm using AbstractGP to setup a Bayesian Optimisation framework. However, AbstractGP shares the array used to construct the GP, such that it's impossible to add new points to it (the ones decided by the acquisition function).

Ideally, the steps would be:

  1. Sample random points X
  2. Evaluate objective function (expensive) on X, get pairs (X,Y)
  3. Fit GP to those points
  4. Evaluate surrogate function on new set of points, X_2
  5. Acquisition function picks which point in X_2 is worth evaluating with the objective
  6. Add the selected point to X

Step 6 fails with error "ERROR: cannot resize array with shared data". The sharing of data happens during step 3.

Steps to reproduce (simplified):

using AbstractGPs
using Random

rng = MersenneTwister(0)

n_samples = 10
x_test = rand(Float64, n_samples) # Initial evaluation grid X
y_test = rand(Float64, n_samples) # Suppose these are the objective function's evaluaton on X
dummy_new_point = 0.111 #dummy point (supposedly coming from the acquisition function)

gp = GP(Matern32Kernel()) #Initialise GP
fx = gp(x_test, 0.01) # Place the GP at the x-locations with some noise

# Fit GP to create the posterior
model = posterior(fx, y_test)  # <--------------------- this shares the data in x_test

# <some acquisition function gives us a new point to evaluate>

push!(x_test, dummy_new_point)  # fails with ERROR: cannot resize array with shared data

Is this intended? Looks like it's rather limiting the usability in a stragithforward bayesian opt setting. If this is intended, what would be the correct course of action here? Should I keep a hard copy of x_test that I update with the new points?

Thanks

theogf commented 11 months ago

It is a known issue and unfortunately something to do with Julia...

However for your usecase you can use the following example: https://juliagaussianprocesses.github.io/AbstractGPs.jl/stable/concrete_features/#Sequential-Conditioning which should be faster as well

filipporemonato commented 11 months ago

Hey, thanks a lot for the blazing fast reply. I had a suspicion it was coming from Julia (facepalm...). So I guess my only option (while keeping the seup above) is doing as I suggested, by maintaing a hard copy of X at all times and updating that one. Thank you so much for the suggestion about sequential coding, I'll look into it!

filipporemonato commented 11 months ago

Btw, sorry for going off-topic here, but would you have a suggestion in the 2D case? The actual objective I need to optimise has 2 input variables. I'm not sure how to set up that case.

devmotion commented 11 months ago

Here's a longer explanation of the problem: https://github.com/JuliaGaussianProcesses/KernelFunctions.jl/issues/319 And this is the underlying Julia issue: https://github.com/JuliaLang/julia/issues/33143

theogf commented 11 months ago

Btw, sorry for going off-topic here, but would you have a suggestion in the 2D case?

It should work in the exact same way as the example mentioned before, just keep recomputing the posterior on the new data point you have

filipporemonato commented 11 months ago

Thanks a lot both of you @theogf and @devmotion . I believe this issue can be closed.