\section{Introduction}

For many applications in sciences and engineering, such as emulation-based studies, design of experiments, and automated machine learning, the goal is to optimize a complex black-box function $f(\boldsymbol{x})$ in a $D$-dimensional space, for which we have limited prior knowledge. The main challenge in such optimization problems is that we aim to efficiently find global optima rather than local optima, while the objective function $f$ is often gradient-free, multimodal, and computationally expensive to evaluate. Bayesian optimization (BO) offers a powerful statistical approach to these problems, particularly when the observation budgets are limited \cite{garnett2023bayesian,shahriari2015taking,gramacy2020surrogates}. A typical BO framework consists of two components to balance exploitation and exploration: the surrogate and the acquisition function (AF). The surrogate is a probabilistic model that  allows us to estimate $f(\boldsymbol{x})$ with uncertainty at a new location $\boldsymbol{x}$, and the AF is used to determine which location to query next.

Gaussian process (GP) regression is the most widely used surrogate for BO \cite{gramacy2020surrogates,williams2006gaussian}, thanks to its appealing properties in providing analytical derivations and uncertainty quantification (UQ). The choice of kernel/covariance function is a critical decision in GP models; for multidimensional BO problems, perhaps the most popular kernel is the ARD (automatic relevance determination)---Squared-Exponential (SE) or Mat\'{e}rn kernel \cite{williams2006gaussian}. Although this specification has certain numerical advantages and can help automatically learn the importance of input variables, a key limitation is that it implies/assumes that the underlying stochastic process is both stationary and separable, and the value of the covariance function between two random points quickly goes to zero with the increase of input dimensionality. These issues can be problematic for complex real-world processes that are nonstationary and nonseparable, as estimating the underlying function with a simple ARD kernel would require a large number of observations.
A potential solution to address this issue is to use more flexible kernel structures. The additive kernel, for example, is designed to characterize a more ``global'' and nonstationary structure by restricting variable interactions \cite{duvenaud2011additive}, and it has demonstrated great success in solving high-dimensional BO problems (see, e.g., \cite{kandasamy2015high, li2016high,rolland2018high}). However, in practice using additive kernels requires strong prior knowledge to determine the proper interactions and involves many kernel hyperparameters to learn \cite{binois2022survey}.

In this paper, we propose to use \emph{Bayesian Kernelized Tensor Factorization} (BKTF) \cite{lei2022bayesianBKMF,lei2022bayesianBCKL} as a flexible and adaptive surrogate model for BO in a $D$-dimensional Cartesian product space. BKTF is originally developed for modeling multidimensional spatiotemporal data with UQ, for tasks such as spatiotemporal kriging/cokriging. This paper adapts BKTF to the BO setting, and our key idea is to characterize the multivariate objective function $f\left(\boldsymbol{x}\right)=f\left(x_1,\ldots,x_D\right)$ for a specific BO problem using low-rank tensor CANDECOMP/PARAFAC (CP) factorization with certain basis functions. Unlike other decomposition methods that either define the basis functions or learn deterministic basis functions, BKTF uses a hierarchical Bayesian framework---GP priors are used to model the basis functions and hyperpriors are used to model kernel hyperparameters---to achieve  UQ.

\begin{figure*}[t]
\centering
\includegraphics[width=0.95\textwidth]{graphical/Intro.pdf}
\caption{BO for a 2D function: (a) true function surface, where the global maximum is marked; (b) comparison between BO models using GP surrogates (with two acquisition functions) and BKTF with 30 random initial observations, averaged over 20 replications; (c) specific results of one run, including the final mean surface for $f$, where green dots denote the locations of the search points during optimization, and the corresponding AF surface.}
\label{fig_Intro}
\end{figure*}

Figure~\ref{fig_Intro} shows the  comparison between BKTF and GP surrogates when optimizing a 2D function that is nonstationary, nonseparable, and multimodal. The details of this function and the BO experiments are provided in Appendix~\ref{appendix-introF}. For this case, GP becomes ineffective in finding the global solution, while BKTF offers great flexibility and adaptability to characterize the multidimensional process from limited data. It should be noted that, different from GP-based surrogate models, BKTF no longer has an analytical posterior; however, efficient inference and acquisition can be achieved through Markov chain Monte Carlo (MCMC), in which we can update basis functions using Gibbs sampling and kernel hyperparameters using slice sampling \cite{lei2022bayesianBCKL}. For the optimization, we can first use MCMC samples to approximate the posterior distribution of the whole tensor and then naturally use the upper confidence bound (UCB) as AF. This process is feasible for many real-world applications with a limited tensor product space, such as experiment design and automatic machine learning. We conduct extensive experiments on both standard optimization and hyperparameter tuning tasks. Our results show that BKTF achieves fast global search for optimizing complex objective functions under limited initial data and observation budgets. The remainder of this paper is organized as follows. We review related studies in Section~\ref{sec:RelatedWork}, introduce notations and basic BO procedures in Section~\ref{sec:preliminary}; in Section~\ref{sec:Methodology}, we present BKTF for BO and also discuss the kernel representations of different surrogates; in Section~\ref{sec:Experiments}, we conduct extensive experiments to evaluate BKTF as a new surrogate model; Section~\ref{sec:Conclusion} summarizes this work and discusses future research directions.

\section{Related Work}\label{sec:RelatedWork}

The key of BO is to effectively characterize the posterior distribution of the objective function from a limited number of observations. The most relevant work to our study is the \textit{Bayesian Kernelized Factorization} (BKF) framework, which has been mainly used for modeling large-scale and multidimensional spatiotemporal data with UQ. The key idea is to parameterize the multidimensional stochastic processes using a factorization model, in which specific priors are used to encode spatial
and temporal dependencies. Signature examples of BKF include spatial dynamic factor model (SDFM) \cite{lopes2008spatial}, variational Gaussian process factor analysis (VGFA) \cite{luttinen2009variational}, and Bayesian kernelized matrix/tensor factorization (BKMF/BKTF) \cite{lei2022bayesianBKMF,lei2022bayesianBCKL,lei2021scalable}. In SDFM, independent GP priors are used to model each column in the spatial factor matrix to encode spatial dependency. In VGFA, spatial and temporal factors are modeled with GP priors from a list of potential structures (e.g., exponential, squared-exponential, and quasi-periodic kernels). BKMF and BKTF also incorporate GP priors to model spatial and temporal factors; in addition, for multivariate data with more than one attribute, BKTF also introduces a Wishart prior to modeling the factors that encode the dependency among attributes. A key difference among these methods is how inference is performed. SDFM and BKMF/BKTF are fully Bayesian hierarchical models and they rely on MCMC for model inference, where the factors can be updated with Gibbs sampling with conjugate priors; for learning the posterior distributions of kernel hyperparameters, SDFM uses the Metropolis-Hastings sampling, while BKMF/BKTF uses the more efficient slice sampling. On the other hand, VGFA uses variational inference to learn factor matrices, while kernel hyperparameters are learned through maximum a posteriori (MAP) estimation without UQ. Overall, recent studies have demonstrated the superior performance of BKTF in modeling multidimensional spatiotemporal processes with high-quality UQ for 2D and 3D spaces \cite{lei2022bayesianBCKL} and conducting tensor regression \cite{lei2021scalable}; in this paper, we are interested in examining the performance of BKTF as a surrogate model for multidimensional BO problems.

The proposed BKTF surrogate models the objective function using low-rank tensor factorization. This specification is closely related to multidimensional Karhunen-Lo\`{e}ve (KL) expansion \cite{wang2008karhunen}, which can be considered principle component analysis (PCA) for stochastic (spatial, temporal, and spatiotemporal) processes or functional principle component analysis (FPCA). The empirical analysis of KL expansion is also known as proper orthogonal decomposition (POD). With a known kernel/covariance function, truncated KL expansion allows us to approximate the underlying random process using a set of eigenvalues and eigenfunctions derived from the kernel function. Numerical KL expansion is often referred to as the Garlekin method, and the key challenge is to obtain or select a set of orthogonal basis functions. Practically, the basis functions are often chosen as prespecified and deterministic functions, such as using the Fourier basis, wavelet basis, orthogonal polynomials, B-splines, empirical orthogonal functions, and Laplacian eigenfunctions (see, e.g., \cite{beylkin2009multivariate}, \cite{wikle1999dimension}, \cite{chevreuil2015least}, \cite{solin2020hilbert}). BKTF uses a hierarchical Bayesian framework to better capture uncertainty, in which GP priors are used to model the latent factors (i.e., basis functions are also random processes) on different dimensions, and hyperpriors are introduced on the kernel hyperparameters. Therefore, we can consider BKTF as a fully Bayesian version of multidimensional KL expansion for stochastic processes with unknown kernel/covariance from partially observed data, however, without imposing any orthogonality constraint on the basis functions.

\section{Preliminaries}\label{sec:preliminary}

Throughout this paper, we use lowercase letters to denote scalars, e.g., $x$,  boldface lowercase letters to denote vectors, e.g., $\boldsymbol{x}=\left(x_1,\ldots,x_D\right)^{\top}\in\mathbb{R}^{D}$, and boldface uppercase letters to denote matrices, e.g., $\boldsymbol{X}\in\mathbb{R}^{M\times N}$. For a matrix $\boldsymbol{X}$, we denote its determinant by $\det{\left(\boldsymbol{X}\right)}$. We use $\bd{I}_{N}$ to represent an identity matrix of size $N$. Given two matrices $\boldsymbol{A}\in\mathbb{R}^{M\times N}$ and $\boldsymbol{B}\in\mathbb{R}^{P\times Q}$, the Kronecker product is defined as $\boldsymbol{A}\otimes\boldsymbol{B}=
\footnotesize{\begin{bmatrix}
  a_{1,1}\boldsymbol{B} & \cdots & a_{1,N}\boldsymbol{B} \\
  \vdots & \ddots & \vdots \\
  a_{M,1}\boldsymbol{B} & \cdots & a_{M,N}\boldsymbol{B}
\end{bmatrix}}\in\mathbb{R}^{MP\times NQ}$. The outer product of two vectors $\boldsymbol{a}$ and $\boldsymbol{b}$ is denoted by $\boldsymbol{a}\circ\boldsymbol{b}$. The vectorization operation $\operatorname{vec}(\bd{X})$ stacks all column vectors in $\bd{X}$ as a single vector. Following the tensor notation in \cite{kolda2009tensor}, we denote a third-order tensor by $\boldsymbol{\mathcal{X}}\in\mathbb{R}^{M\times N\times P}$ and its mode-$k$ $(k=\text{1},\text{2},\text{3})$ unfolding by $\boldsymbol{X}_{(k)}$, which maps a tensor into a matrix. Higher-order tensors can be defined in a similar way.

Let $f:\mathcal{X}=\mathcal{X}_1\times\ldots\times\mathcal{X}_D\rightarrow\mathbb{R}$ be a black-box function that could be nonconvex, derivative-free, and expensive to evaluate. BO aims to address the global optimization problem:
\begin{equation}
\begin{aligned}
\boldsymbol{x}^{\star}=\argmax_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x}),\ \  f^{\star}=\max_{\boldsymbol{x}\in\mathcal{X}}f(\boldsymbol{x})=f(\boldsymbol{x}^{\star}).
\end{aligned}
\end{equation}
BO solves this problem by first building a probabilistic model for $f(\boldsymbol{x})$ (i.e., surrogate model) based on initial observations and then using the model to decide where in $\mathcal{X}$ to evaluate/query next. The overall goal of BO is to find the global optimum of the objective function using as few evaluations as possible. Most BO models rely on a  GP prior for $f(\boldsymbol{x})$ to achieve prediction and UQ:
\begin{equation}
\begin{aligned}
f(\boldsymbol{x})=f(x_1,x_2,\ldots,x_D)\sim\mathcal{GP}\left(m\left(\boldsymbol{x}\right),k\left(\boldsymbol{x},\boldsymbol{x}'\right)\right), \  x_d\in \mathcal{X}_d, \ d=1,\ldots,D,
\end{aligned}
\end{equation}
where $k$ is a valid kernel/covariance function, and $m()$ is a mean function that can be generally assumed to be 0. Given a finite set of initial observations $\left\{\boldsymbol{x}_{i}\right\}_{i=1}^{n}$ with $\boldsymbol{x}_{i}=\left(x_1^i,\ldots,x_D^i\right)^{\top}$, the vector of function values $\boldsymbol{f}=\left(f(\boldsymbol{x}_{1}),\ldots,f(\boldsymbol{x}_{n})\right)^{\top}$ has a multivariate Gaussian distribution $\boldsymbol{f}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{K}\right)$, where $\boldsymbol{K}$ denotes the $n\times n$ covariance matrix. For a set of observations $\mathcal{D}_{n}=\left\{\boldsymbol{x}_i,y_i\right\}_{i=1}^{n}$ with i.i.d. Gaussian noise, i.e., $y_{i}=f(\boldsymbol{x}_{i})+\epsilon_i$ where $\epsilon_i\sim\mathcal{N}(0,\tau^{-1})$, GP gives an analytical posterior distribution of $f(\boldsymbol{x})$ at an unobserved point $\boldsymbol{x}^{\ast}$:
\begin{equation}
\left.f(\boldsymbol{x}^{\ast})\;\middle|\;\right.\mathcal{D}_{n}\sim\mathcal{N}\left(\boldsymbol{k}_{\boldsymbol{x}^{\ast}\boldsymbol{X}}\left(\boldsymbol{K}+\tau^{-1}\boldsymbol{I}_{n}\right)^{-1}\boldsymbol{y},~k_{\boldsymbol{x}^{\ast}\boldsymbol{x}^{\ast}}-\boldsymbol{k}_{\boldsymbol{x}^{\ast}\boldsymbol{X}}\left(\boldsymbol{K}+\tau^{-1}\boldsymbol{I}_n\right)^{-1}\boldsymbol{k}_{\boldsymbol{x}^{\ast}\boldsymbol{X}}^{\top}\right),
\end{equation}
where $k_{\boldsymbol{x}^{\ast}\boldsymbol{x}^{\ast}}$, $\boldsymbol{k}_{\boldsymbol{x}^{\ast}\boldsymbol{X}}\in\mathbb{R}^{1\times n}$ are variance of $\boldsymbol{x}^{\ast}$, covariances between $\boldsymbol{x}^{\ast}$ and $\left\{\boldsymbol{x}_i\right\}_{i=1}^{n}$, respectively, and $\boldsymbol{y}=\left[y_1,\ldots,y_n\right]^{\top}$.
\begin{wrapfigure}[10]{R}{0.44\textwidth}
\vspace{-0.8cm}
  \begin{algorithm}[H]
\SetCustomAlgoRuledWidth{0.44\textwidth}
    \caption{Basic BO process}
    \label{alg_BO}
    \KwIn{Initial dataset $\mathcal{D}_{0}$ and a trained surrogate model; total budget $N$.}
    \For{$n=1,\ldots,N$}
    {Compute the posterior distribution of $f$ using all available data; \\
    Find next evaluation point $\boldsymbol{x}_{n}\in\mathbb{R}^{D}$ by optimizing the AF; \\
    Augment data $\mathcal{D}_{n}=\mathcal{D}_{n-1}\cup\{\boldsymbol{x}_{n},f(\boldsymbol{x}_{n})\}$. \\
    }
  \end{algorithm}
\end{wrapfigure}

Based on the posterior distributions of $f$, one can compute an AF, denoted by $\alpha:\mathcal{X}\rightarrow\mathbb{R}$, for a new candidate $\boldsymbol{x}^{\ast}$ and evaluate how promising $\boldsymbol{x}^{\ast}$ is. In BO, the next query point is often determined by maximizing a selected/predefined AF, i.e., $\boldsymbol{x}_{n+1}=\argmax_{\boldsymbol{x}\in\mathcal{X}}\alpha\left(\boldsymbol{x}\;\middle|\;\mathcal{D}_{n}\right)$. Most AFs are built on the predictive mean and variance; for example, a commonly used AF is the \textbf{expected improvement} (EI) \cite{garnett2023bayesian}:
\begin{equation}   \alpha_{\text{EI}}\left(\boldsymbol{x}\;\middle|\;\mathcal{D}_{n}\right)=\sigma(\boldsymbol{x})\varphi\left(\frac{\Delta(\boldsymbol{x})}{\sigma(\boldsymbol{x})}\right)+\left|\Delta(\boldsymbol{x})\right|\Phi\left(\frac{\Delta(\boldsymbol{x})}{\sigma(\boldsymbol{x})}\right),
\end{equation}
where $\Delta(\boldsymbol{x})=\mu(\boldsymbol{x})-f_n^{\star}$ is the expected difference between the proposed point $\boldsymbol{x}$ and the current best solution, $f_{n}^{\star}=\max_{\boldsymbol{x}\in\{\boldsymbol{x}_i\}_{i=1}^{n}}f(\boldsymbol{x})$ denotes the best function value obtained so far; $\mu(\boldsymbol{x})$ and $\sigma(\boldsymbol{x})$ are the predictive mean and predictive standard deviation at $\boldsymbol{x}$, respectively; and $\varphi(\cdot)$ and $\Phi(\cdot)$ denote the probability density function (PDF) and the cumulative distribution function (CDF) of standard normal, respectively. Another widely applied AF for maximization problems is the \textbf{upper confidence bound} (UCB) \cite{auer2002using}:
\begin{equation}\label{eq:GPucb}
\alpha_{\text{UCB}}\left(\boldsymbol{x}\;\middle|\;\mathcal{D}_{n},\beta\right)=\mu(\boldsymbol{x})+\beta\sigma(\boldsymbol{x}),
\end{equation}
where $\beta$ is a tunable parameter that balances exploration and exploitation. The general BO procedure can be summarized as Algorithm~\ref{alg_BO}.

\section{Methodology}\label{sec:Methodology}

\subsection{Model Specification}

Before introducing BKTF, we first construct a $D$-dimensional Cartesian product space corresponding to the search space $\mathcal{X}$. We define it over $D$ sets $\{S_{1},\ldots,S_{D}\}$ and denote as $\prod_{d=1}^{D}S_{d}$: ${S}_1\times\cdots\times{S}_{D}=\left\{(s_1,\ldots,s_D)\;\middle|\;\forall{d}\in\{1,\ldots,D\},s_d\in S_d\right\}$. For $d=1,\ldots,D$, the coordinates set $S_d$ is formed by $m_d$ interpolation points that are distributed over of $\mathcal{X}_d$, represented by $\boldsymbol{c}_{d}=\left\{c_{1}^{d},\ldots,c_{m_d}^{d}\right\}$, i.e., $S_d=\left\{c_{i}^{d}\right\}_{i=1}^{m_d}$. The size of $S_d$ becomes $|S_d|=m_d$, and the entire space owns $\prod_{d=1}^{D}\left|S_d\right|$ samples. Note that $S_d$ could be either uniformly or irregularly distributed.

We randomly sample an initial dataset including $n_0$ input-output data pairs from the pre-defined space, $\mathcal{D}_{0}=\left\{\boldsymbol{x}_{i},y_i\right\}_{i=1}^{n_0}$ where $\left\{\boldsymbol{x}_{i}\right\}_{i=1}^{n_0}$ are located in $\prod_{d=1}^{D}S_d$, and this yields an incomplete $D$-dimensional tensor $\boldsymbol{\mathcal{Y}}\in\mathbb{R}^{|S_1|\times\cdots\times|S_D|}$ with $n_0$ observed points. BKTF approximates the entire data tensor $\boldsymbol{\mathcal{Y}}$ by a kernelized CANDECOMP/PARAFAC (CP) tensor decomposition:
\begin{equation} \label{eq:bktf}
\boldsymbol{\mathcal{Y}}=\sum_{r=1}^{R}\lambda_r\cdot\boldsymbol{g}_{1}^{r}\circ\boldsymbol{g}_{2}^{r}\circ\cdots\circ\boldsymbol{g}_{D}^{r}+\boldsymbol{\mathcal{E}},
\end{equation}
where $R$ is a pre-specified tensor CP rank, $\boldsymbol{\lambda}=\left(\lambda_1,\ldots,\lambda_R\right)^{\top}$ denote weights that capture the magnitude/importance of each rank in the factorization, $\boldsymbol{g}_{d}^{r}=\left[g_d^r(s_d):s_d\in S_d\right]\in\mathbb{R}^{|S_d|}$ denotes the $r$th latent factor for the $d$th dimension, entries in $\boldsymbol{\mathcal{E}}$ are i.i.d. white noises from  $\mathcal{N}(0,\tau^{-1})$. The function approximation for $\boldsymbol{x}=\left(x_{1},\ldots,x_{D}\right)^{\top}$ can be written as:
\begin{equation}\label{eq:functional_decomposition}
f(\boldsymbol{x})=\sum_{r=1}^{R} \lambda_r g_{1}^{r}\left(x_{1}\right)g_2^r\left(x_{2}\right)\cdots g_{D}^{r}\left(x_{D}\right)=\sum_{r=1}^{R}\lambda_r\prod_{d=1}^{D}g_{d}^{r}\left(x_{d}\right).
\end{equation}

For priors, we assume $\lambda_r\sim\mathcal{N}\left(0,1\right)$ for $r=1,\ldots,R$ and use a GP prior on the latent factors:
\begin{equation}
\left.g_{d}^{r}\left(x_d\right)\;\middle|\; l_d^r\right.\sim\mathcal{GP}\left(1,k_{d}^{r}\left(x_d,x_d';l_d^r\right)\right),~r=1,\ldots,R,~d=1,\ldots,D,
\end{equation}
where $k_d^r$ is a valid kernel function. The reason we use ``1'' as the mean is that we would expect $g_d^r\approx 1$ if variable $x_d$ is not important/relevant. We fix the variances of $k_{d}^{r}$ as $\sigma^2=1$, and only learn the length-scale hyperparameters $l_{d}^{r}$, since the variances of the model can be captured by $\boldsymbol{\lambda}$. In the learning process, for simplicity, we can also assume the lengthscale parameters to be identical, i.e., $l_{d}^{1}=l_{d}^{2}=\ldots=l_{d}^{R}=l_d$, for each dimension $d$. The prior for the corresponding latent factor $\boldsymbol{g}_{d}^{r}$ is then a Gaussian  distribution: $\boldsymbol{g}_{d}^{r}\sim\mathcal{N}\left(\boldsymbol{0},\boldsymbol{K}_{d}^{r}\right)$, where $\boldsymbol{K}_{d}^{r}$ is the $|S_d|\times|S_d|$ correlation matrix computed from $k_{d}^{r}$. We place Gaussian hyper-priors on the log-transformed kernel hyperparameters to ensure positive values, i.e., $\log\left(l_{d}^{r}\right)\sim\mathcal{N}\left(\mu_{l},\tau_{l}^{-1}\right)$. For noise precision $\tau$, we assume a conjugate Gamma prior $\tau\sim\text{Gamma}\left(a_0,b_0\right)$.

For observations, based on Eq.~\eqref{eq:functional_decomposition} we assume each $y_i$ in the initial dataset $\mathcal{D}_{0}$ to be:
\begin{equation}
\left.y_i\;\middle|\; \{g_d^r\left(x_d^i\right)\},\{\lambda_r\},\tau\right.\sim\mathcal{N}\left(f\left(\boldsymbol{x}_i\right),\tau^{-1}\right).
\end{equation}

\subsection{Kernel Representation}

If we consider $\{\lambda_r,g_1^r,\ldots,g_D^r\}_{r=1}^R$ as parameters and rewrite the functional decomposition in Eq.~\eqref{eq:functional_decomposition} as a linear function $f\left(\boldsymbol{x};\{\xi_r\}\right)=\sum_{r=1}^{R}\xi_r|\lambda_r|\prod_{d=1}^{D}g_{d}^{r}\left(x_{d}\right)$ with $\xi_r\overset{\mathrm{iid}}\sim \mathcal{N}\left(0,1\right)$, we can marginalize $\{\xi_r\}$ and obtain a fully symmetric multilinear kernel/covariance function for any two data points  $\boldsymbol{x}=\left(x_{1},\ldots,x_{D}\right)^{\top}$ and $\boldsymbol{x}'=\left(x_{1}',\ldots,x_{D}'\right)^{\top}$:
\begin{equation}\label{eq:k2}
k\left(\boldsymbol{x},\boldsymbol{x}';\{\lambda_r,g_1^r,\ldots,g_D^r\}_{r=1}^R\right) = \sum_{r=1}^R \lambda_r^2 \left[\prod_{d=1}^D g_d^r \left(x_d\right) g_d^r \left(x_d'\right) \right].
\end{equation}

As can be seen, the kernel function is actually parameterized by $\{\lambda_r,g_1^r,\ldots,g_D^r\}_{r=1}^R$. There are some properties to highlight: (i) the kernel is \textbf{nonstationary} since the value of $g_d^r(\cdot)$ is location-specific, and (ii) the kernel is \textbf{nonseparable} as long as $R>1$. Therefore, this specification is very different from traditional GP surrogates:
\begin{equation}\notag
\begin{cases}
    \text{GP with SE-ARD:}& k\left(\boldsymbol{x},\boldsymbol{x}'\right)=\sigma^2\prod_{d=1}^{D}k_{d}\left(x_{d},x_{d}'\right), \\
    & \text{kernel is stationary and separable}\\
    \text{additive GP (1st/2nd order): }& k\left(\boldsymbol{x},\boldsymbol{x}' \right) = \sum_{d=1}^D k_d^{\text{1st}}\left(x_{d},x_{d}'\right)+ \sum_{d=1}^{D-1}\sum_{e=d+1}^{D} k_d^{\text{2nd}} \left(x_{d},x_{d}'\right)k_e^{\text{2nd}} \left(x_{e},x_{e}'\right),\\
    & \text{kerenl is stationary and nonseparable}
\end{cases}
\end{equation}
where $\sigma^2$ represents the kernel variance, and each kernel function of $\left\{k_d(\cdot),k_{d}^{\text{1st}}(\cdot),k_{d}^{\text{2nd}}(\cdot),k_{e}^{\text{2nd}}(\cdot)\right\}$ has its own hyperparameters (e.g., length scale and variance). Compared with GP-based kernel specification, the multilinear kernel in Eq.~\eqref{eq:k2} is more flexible and adaptive to the data. From a GP perspective, Eq.~\eqref{eq:k2} has many hyperparameters to learn; however, we can achieve efficient inference of  $\{\lambda_r,g_1^r,\ldots,g_D^r\}_{r=1}^R$ under a tensor factorization framework. Based on the derivation in Eq.~\eqref{eq:k2}, we can consider BKTF as a ``Bayesian'' version of the multidimensional Karhunen-Lo\`{e}ve (KL) expansion \cite{wang2008karhunen}, in which the basis functions $\{g_d^r\}$ are random processes (i.e., GPs) and $\{\lambda_r\}$ are random variables.

\subsection{Model Inference}

Unlike GP, BKTF no longer enjoys an analytical posterior distribution. Based on the aforementioned prior and hyper-prior settings, we adapt the MCMC updating procedure in \cite{lei2022bayesianBKMF, lei2022bayesianBCKL} to an efficient element-wise Gibbs sampling algorithm for model inference. This allows us to accommodate observations that are not located on the grid space $\prod_{d=1}^D S_d$. The detailed derivation of the sampling algorithm is given in Appendix~\ref{appendix-inference}.

\subsection{Prediction and AF Computation}

In each step of function evaluation, we run the MCMC sampling process $K$ iterations for model inference, where the first $K_0$ samples are taken as burn-in and the last $K-K_0$ samples are used for posterior approximation. The predictive distribution for any entry $f^{\ast}$ in the defined grid space conditioned on the observed dataset $\mathcal{D}_0$ can be obtained by the Monte Carlo approximation:
\begin{equation}
p\left(f^{\ast}\;\middle|\;\mathcal{D}_0,\boldsymbol{\theta}_0\right)\approx\frac{1}{K-K_0}\sum_{k=K_0+1}^{K}p\left(f^{\ast}\;\middle|\;\left(\boldsymbol{g}_{d}^{r}\right)^{(k)},\boldsymbol{\lambda}^{(k)},\tau^{(k)}\right),
\end{equation}
where $\boldsymbol{\theta}_0=\{\mu_{l},\tau_{l},a_0,b_0\}$ is the set of all parameters used in hyper-priors. Although direct analytical predictive distribution does not exist in BKTF, the posterior mean and variance estimated from MCMC samples at each location naturally offer us a Bayesian approach to define the acquisition functions.

BKTF provides a fully Bayesian surrogate model. We define a Bayesian variant of UCB as the AF by adapting the predictive mean and variance (or uncertainty) in ordinary GP-based UCB with the values calculated from MCMC sampling. For every MCMC sample after burn-in, i.e., $k>K_0$, we can estimate a output tensor $\tilde{\boldsymbol{\mathcal{F}}}^{(k)}$ over the entire grid space using the latent factors $\left(\boldsymbol{g}_{d}^{r}\right)^{(k)}$ and the weight vector $\boldsymbol{\lambda}^{(k)}$: $\tilde{\boldsymbol{\mathcal{F}}}^{(k)}=\sum_{r=1}^{R}\lambda_{r}^{(k)}\left(\boldsymbol{g}_{1}^{r}\right)^{(k)}\circ\left(\boldsymbol{g}_{2}^{r}\right)^{(k)}\circ\cdots\circ\left(\boldsymbol{g}_{D}^{r}\right)^{(k)}$. We can then compute the corresponding mean and variance tensors of the $K-K_0$ samples $\{\tilde{\boldsymbol{\mathcal{F}}}^{(k)}\}_{k=K_0+1}^{K}$; we denote the mean tensor and variance tensor by $\boldsymbol{\mathcal{U}}$ and $\boldsymbol{\mathcal{V}}$, respectively. The approximated predictive distribution at each point $\boldsymbol{x}$ in the space becomes $\tilde{f}(\boldsymbol{x})\sim\mathcal{N}\left(u({\boldsymbol{x}}),v(\boldsymbol{x})\right)$. Following the definition of UCB in Eq.~\eqref{eq:GPucb}, we define Bayesian UCB (B-UCB) at location $\boldsymbol{x}$ as $\alpha_{\text{B-UCB}}\left(\boldsymbol{x}\;\middle|\;\mathcal{D},\beta,\boldsymbol{g}_{d}^{r},\boldsymbol{\lambda}\right)=u(\boldsymbol{x})+\beta\sqrt{v(\boldsymbol{x})}$. The next search/query point can be determined via  $\boldsymbol{x}_{\text{next}}=\argmax_{\boldsymbol{x}\in{\{\prod_{d=1}^{D}S_d}-\mathcal{D}_{n-1}\}}{\alpha_{\text{B-UCB}}}\left(\boldsymbol{x}\right)$.

\begin{wrapfigure}[18]{R}{0.5\textwidth}
\centering
  \begin{algorithm}[H]
    \SetCustomAlgoRuledWidth{0.45\textwidth}
    \caption{BKTF for BO}
    \label{alg_BKFBO}
    \SetKwBlock{DoParallel}{for $k=K_0+1:K$ do in parallel}{end}
    \KwIn{Initial dataset $\mathcal{D}_{0}$.}
    \For{$n=1:N$}{
    \For{$k=1:K$}{
    \For{$r=1:R$}{

    \For{$d=1:D$}{
    Draw kernel length-scale hyperparameters $l_{d}^{r}$; \\
    Draw latent factors $\left(\boldsymbol{g}_{d}^{r}\right)^{(k)}$;
    }}
    Draw model noise precision $\tau^{(k)}$; \\
    Draw weight vector $\boldsymbol{\lambda}^{(k)}$; \\
    \uIf{$k>K_0$}{
    Compute and collect $\tilde{\boldsymbol{\mathcal{F}}}^{(k)}$. \\

    }}

    Compute mean $\boldsymbol{\mathcal{U}}$ and variance $\boldsymbol{\mathcal{V}}$ of $\{\tilde{\boldsymbol{\mathcal{F}}}^{(k)}\}$; \\
    Compute $\alpha_{\text{B-UCB}}\left(\boldsymbol{x}\;\middle|\;\mathcal{D}_{n-1},\beta\right)$ as a tensor; \\
    Find next $\boldsymbol{x}_n=\argmax_{\boldsymbol{x}}\alpha_{\text{B-UCB}}\left(\boldsymbol{x}\;\middle|\;\mathcal{D}_{n-1},\beta\right)$; \\
    Augment the data $\mathcal{D}_{n}=\mathcal{D}_{n-1}\cup\{\boldsymbol{x}_{n},f(\boldsymbol{x}_{n})\}$.
    }
  \end{algorithm}
\end{wrapfigure}

We summarize the implementation procedure of BKTF for BO in Algorithm~\ref{alg_BKFBO}. Given the sequential nature of BO, when a new data point arrives at step $n$, we can start the MCMC with the last iteration of the chains at step $n-1$. The main computational and storage cost of BKTF is to interpolate and save the tensors $\tilde{\boldsymbol{\mathcal{F}}}\in\mathbb{R}^{|S_1|\times\cdots\times|S_{D}|}$ over  $(K-K_0)$ iterations for Bayesian AF estimation. This could be prohibitive when the MCMC sample size is large or the dimensionality of input space is large. To avoid saving the tensors, in practice, we can simply use the maximum values of each entry over the $(K-K_0)$ iterations through iterative pairwise comparison. We adopt this simple AF in our numerical experiments.

\section{Experiments}\label{sec:Experiments}

In this section, we analyze the performance of BKTF for BO through extensive experiments on two problems, optimization for benchmark test functions and hyperparameters tuning of ML algorithms.

\subsection{Optimization for Benchmark Test Functions}\label{subsec:ExpF}

\subsubsection{Experimental settings}

\paragraph{Datasets.}

We test the proposed BKTF model for BO on six benchmark functions that are used for global optimization problems \cite{jamil2013literature}, which are summarized in Table \ref{TestFunctions}. Figure~\ref{fig_testF} shows those functions with 2-dimensional inputs together with the 2D Griewank function. All the selected standard functions are multimodal, more detailed descriptions can be found in Appendix~\ref{appendix-testF}. In fact, we can visually see that the standard Damavandi/Schaffer/Griewank functions in Figure~\ref{fig_testF} indeed have a low-rank structure. For each function, we assume the initial dataset $\mathcal{D}_0$ contains $n_0=D$ observed data pairs, and we set the total number of query points to $N=80$ for 4D Griewank and Hartmann function and $N=50$ for others. We rescale the input search range to $[0,1]$ for all dimensions and normalize the output data using z-score normalization.

\paragraph{Model configuration.}

When applying BKTF on the continuous test functions, we introduce $m_d$ interpolation points $\boldsymbol{c}_d$ in the $d$th dimension of the input space. The values of $m_d$ used for each benchmark function are predefined and given in Table~\ref{TestFunctions}. In practice, setting the resolution grid will require certain prior knowledge (e.g., smoothness of the function); and it also depends on the available computational resources and the number of entries in the tensor grows exponentially with $m_d$. In practice, we find that setting $m_d=10\sim 100$ is sufficient for most problems. We set the CP rank $R=2$, and for each BO function evaluation run 400 MCMC iterations for model inference where the first 200 iterations are taken as burn-in, i.e., $K=400$ and $K_0=200$. We use Mat\'{e}rn 3/2 kernel as the covariance function for all the test functions. Since we build a fully Bayesian model, the hyperparameters of the covariance functions can be updated automatically from the data likelihood and hyper-prior.

\paragraph{Baselines.}
We compare BKTF with the following BO methods that use GP as the surrogate model:
\begin{itemize}
    \item GP $\alpha_{\text{EI}}$: GP as the surrogate model and EI as the AF in continuous space $\prod_{d=1}^D{\mathcal{X}_d}$;
    \item GP $\alpha_{\text{UCB}}$: GP as the surrogate model with UCB as the AF with $\beta=0.999$, in $\prod_{d=1}^D{\mathcal{X}_d}$;
    \item GPgrid $\alpha_{\text{EI}}$: GP as the surrogate model with EI as the AF, in Cartesian grid space $\prod_{d=1}^{D} S_d$;
    \item GPgrid $\alpha_{\text{UCB}}$: GP as the surrogate model with UCB as the AF, where $\beta=0.999$, in $\prod_{d=1}^{D} S_d$;
\end{itemize}
We use the same  Mat\'{e}rn 3/2 kernel function for all the GP surrogates. For AF optimization in GP $\alpha_{\text{EI}}$ and GP $\alpha_{\text{UCB}}$, we firstly use the DIRECT algorithm and then apply the Nelder-Mead algorithm to further search if there exist better solutions.

\begin{table}[!t]
\footnotesize
  \centering
  \caption{Summary of the studied benchmark functions.}
    \begin{tabular}{l|rrrrr}
    \toprule
    Function & $D$ & Search space & $f(\boldsymbol{x}^{\star})$ & $\boldsymbol{x}^{\star}$ & $m_d$  \\
    \midrule
    \multirow[c]{3}{*}{Branin} & \multirow[c]{3}{*}{2} & \multirow[c]{3}{*}{$[-5,10]\times[0,15]$}  & \multirow[c]{3}{*}{0.3978873} & $(-\pi,12.275)$ & \multirow[c]{3}{*}{14} \\
     & & & & $(\pi,2.275)$ \\
     & & & & $(3\pi,2.425)$ \\

    Damavandi & 2 & $[0,14]^2$ & 0 & $(2,2)$ & 71 \\
    Schaffer & 2 & $[-10,10]^2$ & 0 & $(0,0)$ & 11 \\

    \multirow[c]{2}{*}{Griewank} & 3 & $[-10,10]^3$ & 0 & $(0,0,0)$ & 11 \\
     & 4 & $[-10,10]^4$ & 0 & $(0,0,0,0)$ & 11 \\

    \multirow[c]{2}{*}{Hartmann} & \multirow[c]{2}{*}{6} & \multirow[c]{2}{*}{$[0,1]^6$} & \multirow[c]{2}{*}{-3.32237} & $(0.20169, 0.150011, 0.476874,$ & \multirow[c]{2}{*}{12} \\
     & & & & $0.275332, 0.311652, 0.6573)$ \\
    \midrule
     & \multicolumn{5}{c}{Characteristics} \\
    \midrule
    Branin & \multicolumn{5}{c}{3 global minima, flat} \\
    Damavandi & \multicolumn{5}{c}{multimodal, global minimum located in small area} \\
    Schaffer & \multicolumn{5}{c}{multimodal, global optimum located close to local minima} \\
    Griewank & \multicolumn{5}{c}{multimodal, many widespread and regularly distributed local optima} \\
    Hartmann & \multicolumn{5}{c}{multimodal, multi-input} \\
    \bottomrule
    \end{tabular}
  \label{TestFunctions}
\end{table}

\begin{figure*}[!t]
\centering
\includegraphics[width=0.99\textwidth]{graphical/testF.pdf}
\caption{Illustrations of the tested benchmark functions, where the global optima are marked on the surfaces. Note that for the Schaffer function, we enlarge the search space to $[-50,50]^2$ to clearly show the feature of this function from a more global view.}
\label{fig_testF}
\end{figure*}

\subsubsection{Results}

\paragraph{Effects of model parameters.}

Since the observation data is scarce, the model performance of BKTF highly depends on the hyper-prior settings of the kernel length-scales of the latent factors and the model noise precision when proceeding estimation for the unknown points, i.e., $\boldsymbol{\theta}_0=\left\{\mu_{l},\tau_{l},a_{0},b_{0}\right\}$. A proper hyper-prior becomes rather important in such optimization scenarios. We discuss the effects of hyper-priors on kernel hyperparameters for optimization performance on Branin function in Appendix~\ref{appendix-prior}. We see that for the re-scaled input space, a reasonable setting is to suppose the mean prior of the kernel length-scales is around half of the input domain, i.e., $\mu_l=\log{(0.5)}$. The effects of hyper-priors on other functions are similar, and we choose an appropriate setting relevant to the input range. The hyper-prior on $\tau$ impacts the uncertainty of the latent factors, for example a large model noise assumption results in more variances in the factors. The uncertainty of BKTF is explained by tow components, uncertainty in latent factors and variances in white noise. The role of $\{a_0, b_0\}$ becomes more important when the objective function is complex that BKTF cannot well describe the function with limited observations. Generally we select the priors that make the noise variances not quite large. An example of the uncertainty provided by BKTF is explained in Appendix~\ref{appendix-introF}.

\paragraph{Comparison with baseline models.}

To compare optimization performances of different models on the benchmark functions, we consider the absolute error between the global optimum $f^{\star}$ and the current estimated global optimum $\hat{f}^{\star}$, i.e., $\left|f^{\star}-\hat{f}^{\star}\right|$, w.r.t. the number of function evaluations. We run the optimization 10 times for every test function with a different set of initial observations. The corresponding results are summarized in Figure~\ref{fig_testFComp}, where the medians and associated 25\% and 75\% quartiles of the 10 runs are illustrated. We see that for the 2D functions Branin and Schaffer, BKTF clearly finds the global optima much faster than those GP-based baselines. For the Damavandi function, where the global minimum ($f(\boldsymbol{x}^{\star})=0$) is located at a small sharp area while the local optimum ($f(\boldsymbol{x})=2$) is located at a large smooth area (see Figure~\ref{fig_testF}), we can see that those GP-based models are trapped around the local optima in most cases, i.e., $\left|f^{\star}-\hat{f}^{\star}\right|=2$, and cannot jump out. On the contrary, BKTF explores the global characteristics of the objective function over the entire search space and reaches the global optimum within 10 iterations of function evaluations. For higher dimensional Griewank and Hartmann functions, BKTF successfully arrives at the global optima under the given observation budgets, while GP-based models are prone to be trapped around local optima.

We compare the absolute error between $f^{\star}$ and the final estimated $\hat{f}^{\star}$ in Table~\ref{Synthetic-comTable}. The enumeration-based GP surrogates, i.e., GPgrid $\alpha_{\text{EI}}$ and GPgrid $\alpha_{\text{UCB}}$, perform a little better than direct GP-based search, i.e., GP $\alpha_{\text{EI}}$ and GP $\alpha_{\text{UCB}}$ on the Damavandi function, but worse on other functions. This means that the discretization, to some extent, offers possibilities for searching all the alternative points in the space, since in each function evaluation, every sample in the space is equally compared solely based on the predictive distribution. Overall, BKTF reaches the global optimum for every test function and shows superior performance for complex objective functions with a faster convergence rate.

To intuitively compare the overall performances of different models across multiple experiments/functions, we further estimate performance profiles (Dolan-Mor\'e curves) \cite{dolan2002benchmarking} and compute the area under the curve (AUC) for quantitative analyses. We consider the number of function evaluations to find the global optimum as the performance measure. Let $t_{p,a}$ denote the number of evaluations used by method/solver $a$ to reach the global solution in experiment $p$ (a lower value is better). The value equals to $N_p+100$ if the method cannot find the global optimum with $N_p$ being the given observation budget for experiment $p$. The performance ratio $\gamma_{p,a}=\frac{t_{p,a}}{\min\{t_{p,a}:a\in\mathcal{A}\}}$, where $\mathcal{A}$ represents the set including all comparing models, and the performance profile for each method is the distribution of $p\left(\gamma_{p,a}\leq\rho\right)$ in terms of a factor $\rho$. We set $\rho=1:N_{\text{max}}+1$, where $N_{\text{max}}=\max{N_{p}}$ is the largest observation budget assumed for the compared experiments. In Figure~\ref{fig_testF_PPComp}, we define the problem set $\left\{\mathcal{P}\;\middle|\;\forall{p}\in\mathcal{P}\right\}$ as the 10 runs for every function and draw the performance profiles of each model, also set $\mathcal{P}$ as all 60 experiments across the 6 tested functions and estimate the overall performance profiles. The AUC values of these curves are given in Table~\ref{Synthetic-comTable}, where the AUC values of the overall performance profiles are taken as the metric to compare the overall performances. As can be seen, BKTF obtains the best performance across all the considered functions. In addition, for most of the test functions, the AUC of grid-based baseline models is comparable with those of continuous GP-based models, suggesting that the discretization of the continuous space is feasible to simplify the optimization problem.

\begin{figure*}[!t]
\centering
\includegraphics[width=0.94\textwidth]{graphical/testF_Comp.pdf}
\caption{Comparison of the optimization results on the six test functions. The results show medians with 25\% and 75\% quartiles of 10 runs with different initialization points.}
\label{fig_testFComp}
\end{figure*}

\begin{table}[!t]
\footnotesize
  \centering
  \caption{Comparison of $\left|f^{\star}-\hat{f}^{\star}\right|$ when $n=N$ (mean $\pm$ std.), and AUC of performance profiles on standard benchmark functions.}
    \begin{tabular}{l|rrrrr}
    \toprule
    Function & GP $\alpha_{\text{EI}}$ & GP $\alpha_{\text{UCB}}$ & GPgrid $\alpha_{\text{EI}}$ & GPgrid $\alpha_{\text{UCB}}$ & BKTF $\alpha_{\text{B-UCB}}$  \\
    \midrule
    Branin ($D=2$) & 0.01$\pm$0.02/36.08 & 0.001$\pm$0.002/36.08 & 0.31$\pm$0.62/46.51 & 4.09$\pm$8.19/43.39 & \textbf{0.00}$\pm$0.00/\textbf{49.91} \\
    Damavandi ($D=2$) & 2.00$\pm$0.00/16.36 & 2.00$\pm$0.00/16.36 & 1.60$\pm$0.80/24.14 & 1.80$\pm$0.60/21.29 & \textbf{0.00}$\pm$0.00/\textbf{50.47} \\
    Schaffer ($D=2$) & 0.02$\pm$0.02/44.89 & 0.01$\pm$0.02/44.30 & 0.10$\pm$0.15/38.73 & 0.11$\pm$0.05/36.73 & \textbf{0.00}$\pm$0.00/\textbf{49.66} \\
    Griewank ($D=3$) & 0.09$\pm$0.09/48.90 & 0.16$\pm$0.12/47.69 & 0.23$\pm$0.13/47.69 & 0.31$\pm$0.19/47.69 & \textbf{0.00}$\pm$0.00/\textbf{50.78} \\
    Griewank ($D=4$) & 0.10$\pm$0.07/79.49 & 0.19$\pm$0.12/77.78 & 0.39$\pm$0.19/77.78 & 0.44$\pm$0.19/77.78 & \textbf{0.00}$\pm$0.00/\textbf{80.46}  \\
    Hartmann ($D=6$) & 0.12$\pm$0.07/78.11 & 0.09$\pm$0.09/78.11 & 0.70$\pm$0.70/79.18 & 1.42$\pm$0.66/78.11 & \textbf{0.00}$\pm$0.00/\textbf{80.78} \\
    \midrule
    Overall performance & -/69.86 & -/69.27 & -/71.16 & -/69.66 & -/\textbf{80.34} \\
    \bottomrule
    \multicolumn{5}{l}{{Best results are highlighted in bold fonts.}}
    \end{tabular}
  \label{Synthetic-comTable}
\end{table}

\begin{figure*}[!t]
\centering
\includegraphics[width=0.99\textwidth]{graphical/testF_PPComp.pdf}
\caption{Performance profiles on the standard test functions.}
\label{fig_testF_PPComp}
\end{figure*}

\subsection{Hyperparameter Tuning for Machine Learning}

In this section, we evaluate the performance of BKTF for automatic machine-learning tasks. Specifically, we compare different models to optimize the hyperparameters of two machine learning (ML) algorithms---random forest (RF) and neural network (NN)---on classification for the MNIST database of handwritten digits\footnote{\url{http://yann.lecun.com/exdb/mnist/}} and housing price regression for the Boston housing dataset\footnote{\url{https://www.cs.toronto.edu/~delve/data/boston/bostonDetail.html}}. The details of the hyperparameters that need to learn are given in Appendix~\ref{appendix-tuning}. We assume the number of data points in the initial dataset $\mathcal{D}_{0}$ equals the dimension of hyperparameters need to tune, i.e., $n_0=4$ and $n_0=3$ for RF and NN, respectively. The total budget is $N=50$. We implement the RF algorithms using scikit-learn package and construct NN models through Keras with 2 hidden layers. All other model hyperparameters are set as the default values.

\paragraph{Model configuration.} We treat all the discrete hyperparameters  as samples from a continuous space and then generate the corresponding Cartesian product space $\prod_{d=1}^{D}S_d$. One can interpret the candidate values for each hyperparameter as the interpolation points in the corresponding input dimension. According to Appendix~\ref{appendix-tuning}, the size of the spanned space $\prod S_d$ is $91\times 46\times 64\times 10$ and $91\times 46\times 13\times 10$ for RF classifier and RF regressor, respectively; for the two NN algorithms, the size of parameter space is $91\times 49\times 31$. Similar to the settings on standard test functions, we set the tensor rank $R=2$, assume $K=400$ and $K_0=200$ for MCMC inference, and use Mat\'{e}rn 3/2 kernel as the kernel function.
\paragraph{Baselines.} As for baseline models, other than the GP surrogate based GP $\alpha_{\text{EI}}$ and GP $\alpha_{\text{UCB}}$, we also compare with Tree-structured Parzen Estimator (BO-TPE) \cite{bergstra2011algorithms}, which is a widely applied BO approach for hyperparameter tuning in ML algorithms. For baselines, we exclude the grid-based GP models as sampling the entire grid using GP becomes infeasible.
\paragraph{Results.} We compare the accuracy for MNIST classification and MSE (mean squared error) for Boston housing regression both in terms of the number of function evaluations and still run the optimization processes 10 times with different initial datasets $\mathcal{D}_0$. The medians along with the 25\% and 75\% quantiles of the 10 runs obtained by different BO models are given in Figure~\ref{fig_MLComp}, and the final classification accuracy and regression MSE are compared in Table~\ref{ML-comTable}. For BKTF, we can see from Figure~\ref{fig_MLComp} that the width between the two quartiles of the accuracy and error decreases as more iterations are evaluated, and the median curves present superior convergence rates compared to baselines. For example, BKTF finds the hyperparameters of NN that achieve 100\% classification accuracy on MNIST using less than 4 function evaluations in all 10 runs. Table~\ref{ML-comTable} also shows the proposed BKTF surrogate achieves the best final mean accuracy and regression error with small standard deviations. All these demonstrate the advantage of BKTF as a surrogate for black-box function optimization.

\begin{figure}[!t]
\centering
\subfigure[MNIST classification.]{
\includegraphics[width=0.48\textwidth]{graphical/MNIST_Comp.pdf}
}
\subfigure[Boston housing regression.]{
\includegraphics[width=0.48\textwidth]{graphical/BH_Comp.pdf}
}
\caption{Comparison of the hyperparameter tuning results for automated ML on MNIST classification and Boston housing regression. The figure shows medians with 25\% and 75\% quartiles of 10 runs.}
\label{fig_MLComp}
\end{figure}

\begin{table}[!t]
\footnotesize
  \centering
  \caption{Comparison of final accuracy (for MNIST classification) and MSE (for Boston housing regression).}
    \begin{tabular}{ll|rrrr}
    \toprule
    Dataset & Algorithm & GP $\alpha_{\text{EI}}$ & GP $\alpha_{\text{UCB}}$ &  BO-TPE & BKTF $\alpha_{\text{B-UCB}}$  \\
    \midrule
    \multirow{2}{*}{MNIST} & RF & 94.12$\pm$0.16 & 94.04$\pm$0.19 & 94.14$\pm$0.20 & \textbf{94.44}$\pm$0.15 \\

     & NN & 99.96$\pm$0.03 & 99.99$\pm$0.02 & 99.96$\pm$0.02 & \textbf{100.00}$\pm$0.00 \\

    \midrule
    \multirow{2}{*}{Boston housing} & RF & 26.19$\pm$0.45 & 26.35$\pm$0.35 & 26.27$\pm$0.31 & \textbf{25.03}$\pm$0.18 \\

     & NN & 38.46$\pm$3.31 & 36.78$\pm$1.91 & 36.40$\pm$4.72 & \textbf{30.84}$\pm$1.13 \\

    \bottomrule
    \multicolumn{5}{l}{The values are presented as mean$\pm$std. {Best results are highlighted in bold fonts.}}
    \end{tabular}
  \label{ML-comTable}
\end{table}

\section{Conclusion and Discussions}\label{sec:Conclusion}

In this paper, we propose to use Bayesian Kernelized Tensor Factorization (BKTF) as a new surrogate model for Bayesian optimization. Compared with traditional GP surrogates, the BKTF surrogate is more flexible and adaptive to data given its Bayesian hierarchical nature. The tensor factorization model behind BKTF offers an effective solution to capture global/long-range correlations and cross-dimension correlations. Therefore, it shows superior performance in characterizing complex multidimensional stochastic processes that are nonstationary, nonseparable, and multimodal. The inference of BKTF is achieved through MCMC, which provides a natural solution for acquisition. Experiments on both standard test function optimization and ML hyperparameter tuning confirm the superiority of BKTF as a surrogate for BO. A limitation of BKTF is the fact that we restrict BO to Castersian grid space to leverage tensor factorization; however, for many real-world optimizations, we believe it is not difficult to design a compatible grid space (or its resolution) based on our prior knowledge.

There are several directions to be explored for future research. A key computational issue of BKTF is that we need to reconstruct the posterior distribution for the whole tensor to obtain the AF. This is problematic for high-dimensional problems due to the curse of dimensionality. It would be interesting to see whether we can achieve efficient acquisition directly using the basis functions and corresponding weights without constructing the tensors explicitly. In terms of rank determination, we can introduce the multiplicative gamma process (MGP) prior to automatically learn the rank; this will create a Bayesian nonparametric model that can automatically adapt to the data. In terms of surrogate modeling, we can further integrate a local (short-scale) GP component to construct a more precise surrogate model, as presented in \cite{lei2022bayesianBCKL}. The combined framework would be more expensive in computation, but we expect the combination to provide better UQ performance. In terms of parameterization, we also expect that introducing orthogonality prior to the latent factors (basis functions) will improve the inference. This can be potentially achieved through more advanced prior specifications such as the Matrix angular central Gaussian \cite{jauch2021monte}. In addition, for the tensor factorization framework, it is straightforward to adapt the model to handle categorical variables as input and multivariate output by placing a Wishart prior on the latent factors for the categorical/output dimension.

\subsubsection*{Acknowledgements}

This work was supported in part by the Natural Sciences and Engineering Research Council (NSERC) of Canada Discovery Grant RGPIN-2019-05950, in part by the Fonds de recherche du Qu\'{e}bec--Nature et technologies (FRQNT) Research Support for New Academics \#283653, and in part by the Canada Foundation for Innovation (CFI) John R. Evans Leaders Fund (JELF). M. Lei would like to thank the Institute for Data Valorization (IVADO) for providing the Excellence Ph.D. Scholarship.

\bibliographystyle{unsrt}