I am have been trying to optimize the Gaussian process likelihood function (multivariate Gaussian likelihood) in R (optimx and nloptr) while having some box constraints for my hyper parameter estimates.

I am facing the problem of the solution rushing to the edges of the parameter space. Why this happening is? And what is a good optimizer in R that is good for constrained R optimization?

  • 2
    $\begingroup$ It depends on your objective function and the values of the bound (box) constraints. How did you determine what values to use for those bounds? The explanation could be that the variables on the bounds wind up being optimal in your problem. Stupid question: Are you maximizing or minimizing the likelihood? Many optimizers have a default or exclusive "setting" (behavior) to only do minimization. If the optimizer only does minimization, you need to multiply the objective value provided to the optimizer by -1, so that it will maximize what you want. So, are you actually maximizing likelihood? $\endgroup$ – Mark L. Stone Aug 29 '16 at 19:15
  • $\begingroup$ Assuming your optimization is implemented correctly and box constraints are reasonable -- the best approach to solve this problem is to use HMC sampling of model hyperparameters, and averaging predictions over alternative hyperparameter configurations. $\endgroup$ – Sycorax says Reinstate Monica Aug 29 '16 at 20:09
  • $\begingroup$ There are also purpose-build GP parameter optimizers for R. One is in the DICE family (Deep Inside Computer Experiments). I think there are at least 2 others but I can't recall their names atm. $\endgroup$ – Sycorax says Reinstate Monica Aug 29 '16 at 20:11
  • $\begingroup$ @GeneralAbrial Thank you for your comment, The HMC sampling is extremely time consuming so sometimes it is not an option in my case. I tried to use teh Dice family however they do not offer to insert your own covariance function. I want to know how to just to use the optimization techniques implemented in Dice in R ? Do you have any experience with that $\endgroup$ – Wis Aug 30 '16 at 17:08
  • $\begingroup$ @MarkL.Stone I am doing maximization thus I did multiply by -1. thank you for yoru comment $\endgroup$ – Wis Aug 30 '16 at 17:10

A strong and relatively straightforward routine to use for box-constraints is BOBYQA; it is available in R through the minqa package. It is the default optimiser for the lme4 package when it comes to box-constraints for the evaluation of (generalised) linear mixed models deviances ($-2$ log-likelihoods) so I suspect it will work well with the GP likelihood evaluations too. (I have used for this GP optimisation tasks successfully at some point but I actually found that re-writing my problem in an unconstrained form was more beneficial.)

As it has already been commented the fact the optimisation algorithm rushes to the boundaries it might well be due to the boundary values offering optimal values. I would suggest generating some data that you know coming from a known GP where the optimal parameters are not near the boundaries and then check your optimisation routine's behaviour. Please note that if the optimal parameters are indeed near the parameter space boundaries some asymptotic results might not hold.

  • $\begingroup$ Thank you for your comment, I will try the BOBYQA package, if I am right this algorithm is gradient free and uses quadratic approximations $\endgroup$ – Wis Aug 30 '16 at 17:11
  • 1
    $\begingroup$ BOBYQA is may be a poor choice because the likelihood surface of a GP as a function of hyper-parameters has many local optima. $\endgroup$ – Sycorax says Reinstate Monica Aug 30 '16 at 17:25
  • $\begingroup$ @raw5: Yes, that is one BOBYQA main advantages actually. $\endgroup$ – usεr11852 says Reinstate Monic Aug 30 '16 at 21:59
  • 1
    $\begingroup$ @GeneralAbrial: Local minima can be due to functionally identical models being derived by permuting/exchanging terms or input weighting. If the local optimum is not too different from the global one, probably the difference in performance might even go unnoticed. This is quite noticeable for NN for example. For what is worth, from the packages you comment about: gptk uses gradient information within its own CG solver, GPFDA uses nlminb and GPfit and tgp use R's optim (usually with L-BFGS-B); so clearly derivative-based algorithms for GPs are not a priori poor choices... $\endgroup$ – usεr11852 says Reinstate Monic Aug 30 '16 at 22:39
  • $\begingroup$ @usεr11852 To be clear, I'm not convinced that there's a great way to avoid the local minima problem. My intention was to note that there is, to the best of my knowledge, no panacea by highlighting the diversity of R packages and a potential problem with your suggestion. $\endgroup$ – Sycorax says Reinstate Monica Aug 30 '16 at 22:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.