Yes, it is possible to tune the hyperparameters of a Gaussian process using Bayesian optimization. In fact, there is at least one (very recent) library for fitting approximate GPs using a random features approximation that provides this as an option: https://xgpr.readthedocs.io/en/latest/ , and it is certainly possible to implement something like this using GPyTorch (https://gpytorch.ai/) as well.
The algorithm would at a high level work like this:
Randomly choose several sets of hyperparameter values (e.g. a specific lengthscale, amplitude etc.) and calculate the marginal likelihood for each set.
Fit a Gaussian process model with an RBF kernel (alternatively 5/2-Matern but I would argue RBF is a simple and perfectly acceptable choice for this task), using the hyperparameter values that have already been evaluated as x and the resulting marginal likelihood values as y.
Use the acquisition function to choose the next set of hyperparameters to evaluate, and calculate the marginal likelihood at this next set. (There are several possible acquisition functions -- more than I can cover in a StackOverflow answer; see this paper for details: https://arxiv.org/abs/1807.02811 . Personally, I like Thompson sampling and upper confidence bound (UCB), but not all practitioners would concur I think).
Repeat (2) followed by (3) until convergence or until a maximum number of iterations is reached.
There are in this scheme two GPs: the "main" GP whose hyperparameters we want to tune, and the "surrogate" GP, which predicts marginal likelihood as a function of the selected hyperparameters. The "surrogate" GP is cheap because the number of hyperparameter combinations we evaluate in the course of this optimization is generally small, while the "main" GP is of course expensive because it must fit the full dataset.
As far as choosing the kernel: For the surrogate GP, I would use RBF, although Matern-5/2 is also a possibility; RBF makes a strong assumption that the marginal likelihood of the main GP is a very smooth (infinitely differentiable) function of the hyperparameters of the main GP, while Matern somewhat relaxes this assumption. It would be challenging to directly select a kernel for the main GP using Bayesian optimization, since it is hard to adapt Bayesian optimization to operate on discrete variables (choice of kernel). You could however indirectly use this approach to select a kernel by tuning hyperparameters of the main GP for several different kernels and then prefer the kernel which achieved the best marginal likelihood during hyperparameter tuning.
This strategy is probably best for kernels that have relatively few hyperparameters. For an RBF kernel, for example, you will have to tune 3 hyperparameters, and you could use Bayesian optimization for this. For a spectral mixture kernel which might have many more than 3 hyperparameters, however, Bayesian optimization would traditionally not be preferred for this (it tends to be inefficient in high-dimensional spaces). There has been however some recent work on improving the efficiency of Bayesian optimization for high dimensional spaces, see this paper: https://arxiv.org/abs/1910.01739
Finally, one other thing to think about is the efficiency of this procedure compared with alternatives. For a GP, you can calculate both the marginal likelihood and the gradient of the marginal likelihood in closed form. It is therefore possible to use say L-BFGS with multiple restarts (to avoid becoming trapped in local minima) to tune hyperparameters rather than Bayesian optimization, and this might be more efficient (depending on the kernel and the problem) than Bayesian optimization, especially for kernels like the spectral mixture kernel which have a large number of hyperparameters. So it would also be a good idea to consider the available alternatives for a given problem before adopting this approach.
stan
documentation (mc-stan.org) and Bayesian Data Analysis 3rd edition. $\endgroup$