How to find the "right" hyperparameters for the Gaussian process used Bayesian optimization Suppose we are using Gaussian process as a surrogate model in Bayesian optimization.To compute the acquisition function using Gaussian process we need to know the right hyperparameter of the Gaussian process, for example, the length scale used in the kernel, noise in the kernel etc. One way to compute the hyperparameters would be to use MLE or MAP estimate using observed data. Is there any other known way to do this? Also, is there any prefered way?
 A: An alternative to MAP or MLE estimates is to sample hyperparameters (i'll call them $\theta$) from a distribution proportional to their probability, given your data and some prior assumptions. This is often more robust than point estimates such as MAP or MLE as it avoids assuming the distribution of $\theta$ can be represented by a single point.
Usually this is accomplished using Markov-Chain Monte-Carlo (MCMC). One then "marginalises" the acquisition function over these hyperparameter samples.
In terms of MCMC sampling for Gaussian Processes, there is an excellent three part series by Michael Betancourt (1,2,3) discussing how to do this with Stan, which explains this procedure (and it's nuances) far more effectively than I am capable here (skip to the third instalment if you are already comfortable with Gaussian Process Regression). Basically the idea is to draw lots of samples of your hyperparameters proportional to how likely they are, given your data and some initial beliefs of what the hyperparameters might be. 
What this means for Bayesian optimization is that rather than computing the acquisition function for a single, fixed instance of the hyperparameters, $\theta^*$, one has access to $m$ samples of $\theta$ drawn approximately from $p(\theta | \mathbf{x}, \mathbf{f}, k, \pi_{\theta})$ (with $\mathbf{x}$ and $\mathbf{f}$ being the observation locations and responses of the target function, respectively, $k$ being the kernel and $\pi_{\theta}$ being the prior distribution for $\theta$).
In this case, one optimizes the expectation of the acquisition function over this distribution to obtain the next point, using the $m$ samples of $\theta$ obtained by MCMC. For an arbitrary acquisition function $\alpha(x|\theta)$, this looks like this:
$
\mathbb{E}[\alpha(x)] = \frac{1}{m} \displaystyle\sum_{i=1}^{m} \alpha(x|\theta_i)
$
The disadvantage of this procedure relative to MAP or MLE estimates is a speed penalty: Not only is MCMC a more computationally intensive procedure than a gradient based optimisation, but the acquisition function must now be calculated $m$ times rather than just once. 
(However, this additional effort may be negligible in comparison to the time required to query your target function if you are willing to consider Bayesian optimization)
In terms of empirical results, I'm aware of one comparison here (See "PES-NB" vs PES for point estimates and marginalisation, respectively) to assess the benefits of marginalisation), but I would wager there are more elsewhere.
