  Bayesian optimization (BO) primarily uses Gaussian processes (GP) as the key
surrogate model, mostly with a simple stationary and separable kernel function
such as the squared-exponential kernel with automatic relevance determination
(SE-ARD). However, such simple kernel specifications are deficient in learning
functions with complex features, such as being nonstationary, nonseparable, and
multimodal. Approximating such functions using a local GP, even in a
low-dimensional space, requires a large number of samples, not to mention in a
high-dimensional setting. In this paper, we propose to use Bayesian Kernelized
Tensor Factorization (BKTF) -- as a new surrogate model -- for BO in a
$D$-dimensional Cartesian product space. Our key idea is to approximate the
underlying $D$-dimensional solid with a fully Bayesian low-rank tensor CP
decomposition, in which we place GP priors on the latent basis functions for
each dimension to encode local consistency and smoothness. With this
formulation, information from each sample can be shared not only with neighbors
but also across dimensions. Although BKTF no longer has an analytical
posterior, we can still efficiently approximate the posterior distribution
through Markov chain Monte Carlo (MCMC) and obtain prediction and full
uncertainty quantification (UQ). We conduct numerical experiments on both
standard BO test functions and machine learning hyperparameter tuning problems,
and our results show that BKTF offers a flexible and highly effective approach
for characterizing complex functions with UQ, especially in cases where the
initial sample size and budget are severely limited.
