Problem Bayesian Optimisation

[#1]

In an earlier R&D study I sucessfully used a 4 factor mixture DOE (3 mixture factors and one 3-level categorical "Isomer" effect) to find optimal settings for a chemical formulation with three quality attributes Y1, Y2 and Y3.

I used the case in a DOE course and used simulation formulas to compute the responses Y123. Good models are obtained and optimization was straightforward and OK, see "ANALYSIS MIXTURE DOE CHEMICAL FORMULATION" in attachment.

Starting from a 12 run mixture Kowalski model DOE I am struggling with  Bayesian optimization (see annex "BAYESIAN OPTIMIZATION") This keeps failing because of a very unreliable Y2 model? Why does Bayesian recommends "Max Desirability" with this very unreliable Y2 model? Starting from other screening & space filling DOE's  I keep getting the same unreliable Y2 model.

What is the problem? How to solve?

[#2]

Hi @frankderuyck, 

Not related to your specific use case, but if you have trouble fitting a Gaussian process model based on your data using the Bayesian Optimization platform, I would recommend these two options :

1. Force the exploration of your design space by creating candidate points using MaxPro criterion. This model-agnostic criterion ensures the maximum distance between previous existing points and the new candidate points, which allows you to cover the experimental space efficiently.
2. You can also change default settings regarding the Model Based Augmentation RSquare Threshold (if you increase this threshold from 0,25 (default) to 0,7 for example, you'll increase the exploration part), and possibly adjust the correlation type (Gaussian, Mattern 3/2 and 5/2, or exponential) at the beginning of your BO loop in case one of this correlation type is more adapted to your use case. Mattern 5/2 is a good default choice, but the Gaussian one might provide slightly better results for Y2.

In every case, the practitioner should stay in control of the platform and its outcomes, and choose the option most relevant to the situation (exploration/exploitation). I have presented a use case at a french Discovery Seminar about BO where the "auto" mode (blindly following the recommendations of the platform with the default settings) leads to suboptimal results. As soon as exploration is enforced (because of model fit inadequate and unreliable) through the addition of some exploration points, the model gets the relationships between inputs and outputs, and the exploitation can become much more interesting.
Looking at your response formula, what can be seen is that they are greatly different :

• Y1 has a model with main effects and interactions effects, approximately the same size (main effects slightly larger),
• Y2 has a model with main effects and some large interaction effects (X1*X3 and X2*X3)
• Y3 has a model with large main effects and some small interaction effects.

The formula formats of Y1 and Y3 are "normal" and easy to fit with a Gaussian Process.
For Y2, the very large interaction effects make the GP fit a lot more complex. Kernel equations involve the distance between points (see example with Gaussian Kernel) and the bandwidth (or also called lenghtscale) parameter. This parameter is learned by the GP for each factor/dimension, and informs about how a relative change in each factor affects the response. But as the dimensions are decoupled in the distance calculations (euclidian distance) and in the learning of this parameter, the GP model has hard times figuring out a reasonable bandwidth/lenghtscale parameter for each factor involved in a strong interaction, as it cannot evaluate well how each factor independantly affect the response.
In practice, when faced with a similar situation of very large interactions effects possible, an additive or polynomial kernel would be more suited, as the additive kernel decouples 1st-order, 2nd-order, 3rd order ... effects, and a polynomial kernel is explicitly set up to a fixed degree (in this case 2 or 3). 

Given the options of the JMP Bayesian Optimization platform, I would recommend forcing the platform to generate more exploration points with MaxPro criterion. Please find attached my test, forcing 4 new exploration runs with MaxPro criterion before letting the platform optimize with Maximize expected Improvement option. The optimum found seems to be closer to your objectives than the one you obtain.

Hope this answer will help you,

[#3]

Thanks for indepth explanation Victor, however there still are issues:

(1) the Y2 profiler shows that the Y2 model still is very unreliable. (2)  Important(!): the optimal setting here for Isomer = Isomer 2 which does not correspond with the DOE - and correct - solution that must be Isomer 1. 

Hard to understand why Bayesion is struggling here..

[#4]

I was seeing R2 values close or above 0,7 for Y2 so I wouldn't say the model is very unreliable ?
Bayesian Optimization will not always provide you the global optimum (if existing), but an optimal solution.
Since mixture designs are optimization designs and your categorical isomer factor has not large effect relatively to the mixture factors, there might be other solutions available with different isomer (and different mixture ratios). There could be several explanations to understand the differences : several optima available with similar good performances, and/or the noise added to Y2 could "hide" the optimum to the GP model. You could try removing the noise from the Y2 formula and see if the optimum is the same, as well as checking with the model from the DoE how close the BO solutions are from the optimum you have in the DoE context.

Hope this answer will help you,

[#5]

Profiler Y2 does not look good; DOE model was excellent

[#6]

After some extra runs, Bayesian consistently yields the optimal Isomer 1 solution, the Y2 profiler remains very grey.

Solution is quite equivalent to what I got, my Y2 R² was lower = 0,52. 

[#7]

The Y2 profiler tells me that one can expect -0,48 < Y2 < 1 correct? Should one  ingnore this result and only look at R²?

[#8]

Looking at your BO file, when opening the JMP BO platform, I have a R² = 0,7999 and a Profiler not so bad for Y2:

Same results on my BO file test : R² = 0,7321 and Profiler seems ok:

So which Profiler or at which step are you showing the "bad" Prediction Profiler ?

You could add some runs using the Maximize Multivariate Model Std Dev (or Maximize Bayesian Desirability Std Dev if the batch size is larger than one) to reduce the prediction uncertainty. This can help, but due to this extreme scenario with very large interaction effects, it will be hard for the Gaussian Process to compete with the multivariate regression from the DoE with interaction effects included. You can take a look at how each model handle interaction through the option "Overlaid Interactions" ; you'll see that the interaction traces for (DoE) regression model are far more "diverse" than the ones from Gaussian Process, due to the limitation I explained in my first post.

• Regression model: 
  
• Gaussian Process:
  

The Gaussian Process interpolates the behavior of the response based on the proximity of the points, whereas the regression uses specific effects (main effects, interaction effects, etc...) to describe the phenomenon. The model fitting and "philosophy" is different: assume a model to be verified for regression, whereas with GP you "build the model as you experiment".

Note that even if you're using Bayesian Optimization platform, if you encounter difficulties fitting the GP model even with the addition of exploration points, you can explore more model options using the platform Model Screening with K-folds crossvalidation, and you'll see that regression methods appear at the top of the best model fitting options (here are the results using only the first construction points of your data table) :

Hope this answer may still help you,