Gaussian Process Optimization JuMP.jl

My goal is to design a trajectory (x,y coordinates) that minimizes the variance of a Gaussian Process at a specified set of query points. The objective function is the trace of the covariance matrix from the Gaussian Process conditioned on the measurements received at the x,y locations along the trajectory.

My question is: is this something that JuMP can handle? From this link it seems like the answer may be no. I’ve also looked into TrajectoryOptimization.jl, Optim.jl, NLopt.jl, Optimization.jl, CasADi.jl and a few other packages. JuMP seemed the most promising, but I would love suggestions on packages that are more well suited to this problem.

1 Like

Are there any other constraints? I guess there must be.

You probably shouldn’t use JuMP for this. You’ll also run into problems if your method of fitting the GP is stochastic.

As a next step, I’d take the time to formulate your problem in a bit more detail. It’s hard to offer advice at the moment. A function which takes x and y as input and returns the objective function would be useful, along with a description of any domain restrictions on x and y.

I’m using the Gaussian Process structure outlined in Chapter 15 of Algorithms for Optimization. My objective function is basically Algorithm 15.4 with slight modifications (note this is not a minimum working example, just trying to provide further detail on what I’m trying to do):

function objective_func(x, y, GP)    
    posterior_GP = deepcopy(GP)
    for i in 1:length(x)
        sample_loc = [[x[i],y[i]]]
        sample = GP.m([x[i],y[i]]) # just sampling the mean function
        posterior_GP = posterior(posterior_GP, sample_loc , sample, noise)
    tmp = posterior_GP.KXqX / (posterior_GP.KXX + diagm(posterior_GP.ν))
    S = posterior_GP.KXqXq - tmp*posterior_GP.KXqX'
    νₚ = diag(S) .+ eps() # eps prevents numerical issues
    return sum(νₚ)

In regard to other constraints that are present, it is mainly the dynamic constraints for the trajectory and the starting and ending position. In JuMP syntax, this would look like:

(I recognize this will not work based on how objective_func is defined above which is not compatible with JuMP)

# length of trajectory 
const n = 300

model = Model(Ipopt.Optimizer)

@variables(model, begin
               -1 ≤ xdotdot[1:n] ≤ 1              # acceleration
               -1 ≤ ydotdot[1:n] ≤ 1              # acceleration

# Fix initial conditions
fix(x[1], 0; force=true)
fix(y[1], 0; force=true)

fix(xdot[1], 0; force=true)
fix(ydot[1], 0; force=true)
fix(xdotdot[1], 0; force=true)
fix(ydotdot[1], 0; force=true)

# Fix final conditions
fix(x[n], 100; force=true)
fix(y[n], 100; force=true)

# System dynamics
for j in 2:n
    i = j - 1  # index of previous knot

    @NLconstraint(model, xdot[j] == xdot[i] + xdotdot[i])
    @NLconstraint(model, x[j] == x1[i] + xdot[i])
    @NLconstraint(model, ydot[j] == ydot[i] + ydotdot[i])
    @NLconstraint(model, y[j] == y[i] + ydot[i])

@NLobjective(model, Min, objective_func(x,y))

Tricky! How long does one evaluation of objective_func(x, y) take? Where is the noise coming in your GP? I guess don’t fully understand how that’s created or what it represents.

You essentially have three options:

  • Use a derivative-free method, but this is going to struggle with your dimensionality and the constraints
  • Derive gradients for objective_func, but since the GP is itself a local minimizer, you’ll probably need something like finite differences, which could be expensive in 600-dimensional space
  • Approximate the objective with something tractable and optimize that instead, don’t try to optimize the tricky objective

None of them are very nice, but someone else might have more ideas. This isn’t really my field.

1 Like

One evaluation of objective_func(x, y) takes roughly 0.15 sec. The noise (measurement variance) is the same for all samples just to indicate that a sample does not fully collapse the uncertainty of the GP at the evaluated sample location.

Thanks for the suggestions! I might give those a try and post back here if I find anything that works.