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.
Paper
Similar papers
Peer review
Summary
This paper proposes using Bayesian Kernelized Tensor Factorization (BKTF) as a new surrogate model for Bayesian optimization (BO). BKTF approximates the objective function using a low-rank tensor factorization, with Gaussian process priors placed on the latent factors to capture dependencies and enable uncertainty quantification. The key advantages are the ability to handle complex functions that are non-stationary and non-separable, and the sharing of information across dimensions to enable a more global search compared to standard GP models with local kernels. Inference is performed via MCMC sampling. Experiments on benchmark optimization functions and hyperparameter tuning tasks demonstrate improved performance over GP-based BO, especially when the initial sample size and evaluation budget are limited.
Strengths
The BKTF surrogate is a novel and creative approach to extend BO to handle more complex objective functions. Modeling the objective as a tensor factorization with GP priors on the factors is an elegant way to introduce non-stationarity and multi-scale correlations in a principled Bayesian framework. The method is grounded in a clear mathematical framework, with full details of the model specification and MCMC-based inference procedure provided. Positioning BKTF as a type of deep GP offers useful insight into its expressive power. The paper includes extensive experiments on a range of synthetic functions and real-world hyperparameter tuning tasks. The results convincingly demonstrate the advantages of BKTF over standard GP-BO in terms of optimization performance and sample efficiency, especially in the realistic setting of a very limited evaluation budget. The paper is clearly written, with the methodology explained in detail and the experimental setup and results presented thoroughly. The authors also discuss the limitations of their work, including the scalability challenges and the restriction to a grid-based search space.
Weaknesses
The main weakness is that the proposed BKTF method is only compared against basic GP-based BO with standard kernels. To fully demonstrate the advantage of the BKTF surrogate, comparisons should be made to more advanced GP models such as deep kernel learning, deep GPs, and other scalable GP variants. Without these comparisons, it's difficult to assess how much of the performance gain is due to the specific BKTF approach vs. simply being a more flexible GP. The BKTF model bears significant similarities to existing works on scalable GPs, such as "Gaussian Processes for Big Data" and "Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)", which also use inducing points on grids to obtain tractable approximations. The relationship of BKTF to these methods is not discussed, and it's unclear whether BKTF provides any substantial advantages over these existing scalable GP approaches. The experiments on the synthetic test functions are somewhat limited in their dimensionality (only up to 10d). Given that BO is most useful for optimizing high-dimensional black-box functions, testing on some higher-dimensional benchmarks would be valuable. The scalability of BKTF as the dimensionality and grid size increase is not fully investigated. The MCMC inference procedure may become prohibitively slow for high-dimensional spaces or large evaluation budgets. The paper does not report the wall-clock time of the experiments, which makes it hard to assess the computational feasibility of BKTF in practice, especially compared to alternative approaches. For the hyperparameter tuning experiments in Section 5.2, the strong performance of BKTF with very few iterations seems counterintuitive and is not fully explained, since the BKTF model has many parameters and would be expected to require a substantial amount of data to train effectively. This seems to contradict the results on the synthetic test functions, where GP-EI performs equally well in the first few iterations.
Questions
How does BKTF differ from existing scalable GP methods like "Gaussian Processes for Big Data" and "Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)", which also exploit grid structure and inducing points? Does BKTF offer any advantages over these approaches in terms of performance or flexibility? Beyond standard GPs, have you considered comparing BKTF to other more advanced probabilistic regression models such as neural processes, Bayesian neural networks, or tree-based models like BART or XGBoost? These would seem to be the most relevant competitors for black-box optimization. Could you report the running time of the BKTF method, especially the cost of the MCMC inference as the dimensionality and evaluation budget grow? This is critical for assessing the practical feasibility of the approach. Have you investigated the robustness of BKTF to the choice of grid size and spacing? Does performance degrade gracefully if the grid is misspecified or can this lead to severe failure modes? The strong performance of BKTF on the hyperparameter tuning tasks, especially with very few iterations, is quite surprising and seems to conflict with the synthetic results. Could you provide some intuition or explanation for this? Is there something fundamentally different about these tasks that enables BKTF to perform well with little data, or is it more of an artifact of the experimental setup?
Rating
4
Confidence
4
Soundness
3
Presentation
3
Contribution
2
Limitations
The relationship of BKTF to existing scalable GP methods is not thoroughly discussed. The paper does not make clear how BKTF differs from or improves upon approaches like "Gaussian Processes for Big Data" and "Kernel Interpolation for Scalable Structured Gaussian Processes (KISS-GP)", which also exploit grid structure. A more thorough comparison is needed to justify BKTF as a novel contribution. The cubic scaling of the covariance matrix computations in the number of grid points, which could make BKTF infeasible for high-dimensional or very fine-grained grids. Some discussion of potential ways to scale up BKTF, e.g., by exploiting grid structure or using sparse approximations, would be valuable. The fact that BKTF relies heavily on a sensible grid specification, and may fail badly if the grid is poorly chosen. Some experiments showing the sensitivity of the results to the grid choice would help assess this risk. The lack of comparison to a wider range of flexible surrogate models beyond standard GPs, including more sophisticated GP models as well as other probabilistic regression approaches. The current experiments are not sufficient to establish BKTF as the best choice for BO in practice.
Rebuttal 2nd Round (Reply to Reviewer qG8Q)
Dear reviewer qG8Q, Thank you for your response, and we sincerely appreciate your previous review. We explained the reasoning for **not using scalable GP models (e.g., KISS-GP)** and clarified the **experiments, time cost (which also demonstrates the contribution), and overall contribution** below. - 1. Regarding **not using scalable GP models, e.g., `KISS-GP`** [1]. **(1)** First of all, we re-checked related papers, sorry a **correction** about previous reply: `KISS-GP` does provide predictive mean and variance; however, as stated in the second and third paragraphs of the paper [2]: ``` Although predictive uncertainties are a primary advantage of GP models, they have recently become their computational bottleneck. Historically, the use of GPs has been limited to problems with small datasets, since learning and inference computations naively scale cubically with the number of data points (n). Recent advances in inducing point methods have managed to scale up GP training and computing predictive means to larger datasets. Kernel Interpolation for Scalable Structured GPs (KISS-GP) is one such method that scales to millions of data points. However, these computational savings do not extend naturally to predictive uncertainties. With KISS-GP, computing the predictive covariance between two points requires O(n + m log m) computations, where m is the number of inducing points (see Table 1). While this asymptotic complexity is lower than standard GP inference and alternative scalable approaches, it becomes prohibitive when n is large, or when making many repeated computations. Additionally, drawing samples from the predictive distributions – a necessary operation in many applications – is similarly expensive. ``` **KISS-GP primarily focuses on scalable posterior mean estimation, and the computational advantages diminish when it comes to uncertainty quantification (UQ)**. In their original paper [1], the authors also do not provide posterior variance (UQ) results, focusing instead only on posterior mean. While scalable uncertainty estimation can be achieved through approximation methods such as those in [2], we believe that **scalable GP inference is not the main focus of this work (important)**. We provide a more detailed explanation below. [1] Andrew Gordon Wilson, and Hannes Nickisch. "Kernel interpolation for scalable structured Gaussian processes (KISS-GP)." International conference on machine learning. PMLR, 2015. [2] Pleiss, Geoff, Jacob Gardner, Kilian Weinberger, and Andrew Gordon Wilson. "Constant-time predictive distributions for Gaussian processes." International Conference on Machine Learning, pp. 4114-4123. PMLR, 2018. **(2)** `KISS-GP` was proposed to address the cubic computational cost challenge in GP inference for large-scale data. Specifically, Structured Kernel Interpolation (`SKI/KISS-GP`) [1] which interpolates inducing points on a regularly spaced grid, is **designed for low-dimensional data (dimension < 4, as stated in the paper) and stationary kernels**. A strong/necessary assumption in `KISS-GP` is the use of a stationary and separable kernel function (i.e., utilizing Kronecker structure). For the BO problem considered in this paper, we **do not face scalability issues**, as the number of observations is generally small (a limited budget of fewer than 100 observations, along with severely limited initial data points). The more crucial aspect here is **the specification (representational ability) of the covariance function**. The stationary and separate kernel function in `KISS-GP` limits its capacity to characterize global correlations and, consequently, to find the global optimum of complex nonstationary objective functions. **(3)** Additionally, to reduce the computational cost, scalable GP models often apply certain approximation methods, such as the Kronecker property used in `KISS-GP`. These models should be considered **approximations to the standard GP model in terms of kernel representation** property. With the same kernel function setting and the same observation/training data, **the posterior distribution (UQ) provided by these scalable methods should not surpass that of the original GP models**. Therefore, there is no justification for using these approximated methods, especially when the computational cost of using a standard GP is affordable. - 2. Regarding **Comparison to SOTA methods**. For the optimization tasks on benchmark test functions (Section 5.1), we consider **Bayesian deepGP** [3] as a SOTA model. We will also include grid-based GP combined with a random sampling strategy for the acquisition function as baselines for high-dimensional problems. For the ML hyperparameter tuning tasks (Section 5.2), after carefully reviewed your comments and those of other reviewers, we will include another baseline for mixed inputs, the work in [4], in the revised paper.
Rebuttal 2nd Round (Reply to Reviewer qG8Q) [Continued, Second Part]
We are also reviewing recent literature for more SOTA models and would welcome any suggestions from the reviewer to include in the revised paper. Given the time constraints, we may not be able to provide the updated results in this rebuttal, but we believe the advantage of our proposed BKTF surrogate will remain evident in the considered tasks, as these tasks are specifically chosen to highlight its contributions. [3] Sauer, Annie, Robert B. Gramacy, and David Higdon. "Active learning for deep Gaussian process surrogates." Technometrics 65.1 (2023): 4-18. [4] Ru, Binxin, et al. "Bayesian optimisation over multiple continuous and categorical inputs." International Conference on Machine Learning. PMLR, 2020. - 3. About **Computational Cost**. Please note that we have updated the **wall-clock time comparison** table and the **theoretical model complexity** table **in an additional response provided to all reviewers**. Please refer to the updated tables for more detailed information about the time cost of the proposed model. These time comparison results further demonstrate the advantage of our proposed BKTF surrogate in **optimizing nonstationary complex functions with a severely limited budget and initial observations**. We will also clearly clarify the motivation/focus of our work in the revised paper, as the reviewers suggested. - 4. About **Contribution** of this work. The contribution of this work was highlighted in **Point 1** of the previous General Response, not Point 3. We apologize for the typo. We would like to further emphasize the contribution of this work (as we mentioned in our reply to Reviewer CTJW): to the best of our knowledge, **this is the first time fully Bayesian kernelized low-rank factorization has been introduced to the BO community**. BKTF offers an elegant solution for BO, with a contribution comparable to the work referenced in [5]. BKTF is capable of addressing several challenging problems in BO, particularly in optimizing nonstationary processes and **discrete/categorical variables** with a limited budget. [5] Fusi, Nicolo, Rishit Sheth, and Melih Elibol. "Probabilistic matrix factorization for automated machine learning." Advances in neural information processing systems 31 (2018). We appreciate your understanding and hope these clarifications provide a more precise explanation of our work. If you have any further questions or comments, we would be happy to address them.
Thank you very much for your detailed response. I am still concerned about the novelty and practicality of this work. Particularly, how is low-rank factorization fundamentally different from other methods like KISS-GP, deep GP, and deep kernel learning to improve BO in general? The author claims to be the first in "fully Bayesian kernelized low-rank factorization" in BO but I do find it fundamentally different from implementing a {KISS-GP, deep GP, and deep kernel learning} with some adaption for BO. Thus I will keep my rating.
Reply to Reviewer qG8Q: "Bayesian" and "Kernelized" are the key distinctions of our model BKTF
Thank you for your response. ### Regarding the differences between the proposed surrogate model `BKTF` (Bayesian Kernelized Tensor Factorization) and the methods you mentioned: We would like to emphasize two key terms in the model name that the reviewer may miss: "**Bayesian**" and "**Kernelized**", which corresponds to the **model inference** and **model structure** aspects of the proposed surrogate model `BKTF`. These are the key components that enables `BKTF`'s performance in BO (the key contribution of this model) and distinguish `BKTF` from the methods you referenced. - "**Kernelized**" tensor factorization: In terms of model structure, (1) firstly `BKTF` (Tensor Factorization) utilizes low-rank factorization to model the objective function, capturing global correlations/structures in the data. (2) Secondly, compared to other low-rank models, `BKTF` introduces a GP prior on each column of the latent factors to characterize local correlations in each dimension. In this way, both global and local correlations and dependencies across dimensions can be effectively captured, enabling the model to find the global optimum in the objective function. As stated in Abstract of the paper (lines 8-13): ```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, the information from each sample can be shared not only with neighbors but also across dimensions, thus fostering a more global search strategy.``` We believe the introduced model structure has been clearly formulated in Section 3 (Methodology part) of the paper, see Equations (6)-(8). For example, Equation (8) explains the GP prior placed on the latent factors. - "**Bayesian**": `BKTF` is a fully-Bayesian model: We place priors and hyperpriors on model parameters (e.g., latent factors) and hyperparameters (e.g., kernel hyperparameters), respectively, and develop an efficient MCMC (Markov chain Monte Carlo) algorithm for model inference. This fully Bayesian model allows `BKTF` to avoid overfitting, even with very sparse observations. The fully-Bayesian sampling strategy is also the key reason why `BKTF` can deliver high-quality uncertainty quantification with severely limited observations, which contributes the primary advantage of `BKTF` as a surrogate model. We provide a clear illustration of the model inference process in Appendix A of the paper. The reason for using words like "first" in the rebuttal is to highlight the contribution of introducing `BKTF` **as a surrogate model for BO**. We have avoided using such language in the paper. - Lastly, BO requires high-quality uncertainty quantification from limited data, which is the goal of `BKTF`. This study demonstrates the superiority of `BKTF` as a surrogate model, particularly under conditions of a limited budget and sparse initial observations. `BKTF` offers high-quality uncertainty quantification, effectively captures the correlation structure across dimensions, and provides a straightforward yet effective solution for accommodating both continuous and categorical decision variables. Thank you.
Summary
This paper proposes the Bayesian Kernelized Tensor Factorization (BKTF) as a surrogate for Bayesian optimization (BO). This model uses a CP decomposition to define a set of random basis functions drawn from a GP prior. These latent functions are then weighted by another set of random variables. This defines a probabilistic model that can perform uncertainty quantification, which can then be trained by performing MCMC sampling over the hyperparameters. The acquisition function is calculated by computing the first and second moments of the samples, and calculating the upper confidence bound (UCB). This procedure results in a non-stationary, non-separable model that can capture complex functions. This model is tested on a range of synthetic and ML hyperparameter BO benchmarks, each of which are non-separable functions that exhibit high degrees of interaction between variables.
Strengths
**Originality:** This paper is the first to use kernelized tensors in the Bayesian optimization setting. Gaussian processes are a common surrogate in this setting, and this proposes an alternate surrogate with good arguments for its adoption. **Quality:** The paper demonstrates the performance of the BKTF surrogate on a wide range of benchmarks. The performance is strong, justifying the claims in the paper. Further, many supplementary results are provided to further investigate the modelling decisions made. **Clarity:** The explanation of the BKTF model and fitting process is explained well, providing a clear description of the model. Specifically, the 2D example in Figure 1 provides a clear motivation for modelling functions that have a high degree of interaction. **Significance:** The proposed method is a strong surrogate for Bayesian optimization, one that can model functions with mixed input spaces and high degrees of interaction between variables. This presents a good addition to the range of available surrogates in the field.
Weaknesses
**Comparison to other methods:** The authors claim that a strength of their method is that their method is non-stationary and non-separable. However, I do not feel that the paper sufficiently justifies this explanation for the model's success, for the following reasons: I do not believe that the GP ARD is separable. The authors present that the ARD kernel is the product of $D$ independent kernels in 3.2. However, this is not how these kernels are implemented. Instead, (specifically for stationary kernels) they are expressed as functions of weighted distance [5], where $d=\sqrt{\sum_{d=1}^D (x_d - x'_d)^2/l_d}$. These kernels cannot be written as products of 1-dimensional kernels, and are non-separable. The authors also suggest that these experience the curse of dimensionality as the dimension increases, but this effect is not severe in the low-dimension regime studied in the paper. The authors claim that the flaw of using additive GP kernels is that they: > [require] strong prior knowledge to determine the proper interactions and [involve] a large number of kernel hyperparameters to learn I do not find this argument convincing. For low dimensional problems, additive kernels can include all interactions up to order D, and learn the weighting of each order of interaction [1]. In fact, I believe that the number of kernel hyperparameters is of a similar order to the BKTF method. Additive kernels also work well with MAP estimation of the hyperparameters (especially for the low-dimensional problems investigated here), although I do not see why MCMC could not also be used for additive kernels if the number of kernel hyperparameters is considered too large. The additive GP baseline should therefore be order D, not order 1, to provide a fair comparison against existing non-separable modelling methods. It is unclear why the authors compare to SaaSBO, a technique designed for high dimensional (D>100) spaces that places a strong prior on the lengthscales of the inputs (with the prior belief that few of the inputs are important, which is not the case for the test functions used). Moreover, the authors do not provide comparisons against methods that are designed for non-stationary settings e.g. [1, 2]. (Minor comment) I would like to see MCMC over the GP kernel hyperparameters to obtain a fully Bayesian acquisition function, as in [6]. **Hyperparameter choices with BKTF:** I disagree with the authors that BKTF is robust to rank misspecification. Figure 13 shows that the CRPS doesn't converge until 30 observations for the rank 2 model. Including the initial 30 datapoints, this model is fit on 60 datapoints, for a 2D problem - the GP model provides a much better fit to the data for <60 datapoints. Since these models are used in a BO setting where few datapoints is common, this behaviour suggests that rank *is* an important parameter, and further that the model does not fit well with few datapoints. I would also want to see this experiment repeated on higher dimension test problems, to see if the problem is exacerbated in high dimensions. Moreover, the CRPS of the rank 2 model seems to converge only to that of the GP models, suggesting the performance over GP may not be solely due to modeling quality: this should be further investigated. **Experiments:** It is unclear how the authors handle categorical inputs for the baselines. The authors should be using methods designed for mixed input spaces, such as [4]. Following from the discussion on CRPS, this paper would benefit from some experiments on the quality of the fit of the model. (Minor comment) It would be interesting to see the (arguably more popular) Matern 5/2 kernel compared to the 3/2 kernel, especially for the GP baselines. [1] Snoek et al. "Input Warping for Bayesian Optimization of Non-stationary Functions" [2] Eriksson et al. "Scalable Global Optimization via Local Bayesian Optimization" [3] Duvenaud et al. "Additive Gaussian Processes" [4] Ru et al. "Bayesian optimisation over multiple continuous and categorical inputs" [5] Williams and Rasmussen. "Gaussian Processes for Machine Learning" [6] Snoek et al. "Practical Bayesian Optimization of Machine Learning Algorithms"
Questions
Appendix F claims that BKTF has the lowest computational complexity. Is this for the entire MCMC chain, or just for each sample? How does the wall-clock time compare against these other methods? Why do you use $R=2$ for synthetic experiments, and $R=4$ for the hyperparameter tuning problems? Do you observe similar results when you use $R\in\{6,8\}$? What are the disadvantages of using a higher rank? In E.3, you claim that BKTF learns periodicity of the function. However, a Matern kernel that fits to data sampled from a sine function will appear to have learned periodicity in the input domain. Has the function truly learned periodicity - if you test the function well outside the domain of the function, does it still exhibit periodic behaviour? The authors note that the BKTF model is a special case of a deep GP. Could you provide insight as to why this model still outperforms deepgp? (Minor question) How many MCMC iterations are generally needed to converge? 400 iterations are used in 5.1, 600 in 5.2, and 1000 samples are used in C.2. Do you have any insight for the convergence of the sampling procedure? Do you run multiple chains in parallel? (Minor question) How do you select values for $\{a_0, b_0\}$, and what values do you use in your experiments? (Minor question) BKTFRandom selects a subset of 20,000 points from the domain. However, the search space for the synthetic functions is $(m_d)^D$, which is less than 20,000. Why does BKTFRandom then perform worse than BKTF on these functions?
Rating
6
Confidence
3
Soundness
2
Presentation
3
Contribution
3
Limitations
The authors provide good discussion on the limitations of their approach.
Reply to Reviewer jRPp (Continued, Second Part)
We have updated the relevant table (Table 3 in the paper) accordingly, as provided below, where $n$ is the number of observed data, $m_d$ is the number of interpolation points, $J$ and $J'$ represent the optimization iterations taken for fitting a standard GP and an additive GP model, respectively, $K$ denotes the number of MCMC iterations, and $R$ is the rank. | Model | Complexity | | -------- | -------- | | GP $\alpha_{\text{EI}}$ | $\mathcal{O}\left(Jn^3\right)$ | | GP $\alpha_{\text{UCB}}$ | $\mathcal{O}\left(Jn^3\right)$ | |GPgrid $\alpha_{\text{EI}}$|$\mathcal{O}\left(Jn^3\right)$| |GPgrid $\alpha_{\text{UCB}}$|$\mathcal{O}\left(Jn^3\right)$| | additive GP | $\mathcal{O}\left(J'n^3\right)$ | | BKTF/BKTFrandom | $\min\left\(\mathcal{O}\left(Kn^3\right),\mathcal{O}\left(KR\sum_{d=1}^D m_d^3\right)\right\)$ | | deepGP | $\mathcal{O}\left(K\left(\prod_{d=1}^D m_d\right)^3\right)$ | Please note that this table illustrates the complexity of surrogate models in different BO methods. We will update Table 3 in the revised paper and have also provided this updated table in the comment added for the response to all reviewers. We will also include a table comparing the costs of acquisition function computation. Thank you very much for the detailed comment. (2) $n$ represents the number of observations, which pertains to the cost of GP regression. Sorry the complexity for `BKTFrandom` should be the same as that for `BKTF` and is not dependent on the number of candidate points $n_{\text{random}}$. This was a typo/mistake. (The time provided in wall-clock comparison table also reflects this point.) We have corrected this in the table above and will correct Table 3 of the paper. Thanks for the detailed comment. (3) We agree with the reviewer and will remove the sentence in the revised paper. Thanks a lot for the comment. - **Grid** (1) We set the grid to be uniformly distributed across the search space in each dimension. With this setup, for most test function optimization cases considered in the paper, including the Damavandi function, the global optimum is located within the defined grid. Generally, one can define a grid resolution as finely as needed (as long as the cost remains acceptable), ensuring that the global optimum on the grid is very close to the true optimum, which allows us to find a nearly optimal solution. Additionally, based on the results from BKTF, one can further combine it with a local GP or adjust the grid accordingly to find the global solution. (2) For the GP baselines, whether using a grid or not is unlikely the reason they do not find the global optimum. We provide the results for both continuous non-grid GP and corresponding grid-based GP models with the same grid specification as BKTF, and both non-grid and grid-based GP baselines fail to find the global solution. Additionally, as shown in Table 2 and the updated costs in the above wall-clock time comparison table (where the stopping criterion has been relaxed to a small error as discussed), the continuous GP baselines perform worse compared to BKTF, with a much larger final _regret_ under the limited budget. (3) We did not use adaptive grids in the paper, as the existing results were already better than the baselines. We will consider exploring adaptive grids in future work. - **Dimensionality** (1) Thanks for the comment. Yes, we agree with the reviewer. We will clearly specify the motivation and focus of this work, avoiding mention of high-dimensional problems in the revised paper. (2) (minor) Yes, 10D cannot be stated as "relatively high-dimensional". We will not mention "high/relatively high dimension.." in the paper. For SAASBO, we will consider removing this baseline and instead consider using grid-based GP with a random sample acquisition function (see the below point) as baselines. (3) (minor) Thanks for the suggestion. This is a good idea for comparing grid-based GP in ML hyperparameter tuning tasks. We will include a random sampling strategy for grid-based GP as a baseline for both the ML tasks and the 10D Griewank function optimization. Thanks again for the valuable feedback.
I appreciate the detailed responses to the issues raised. I will therefore increase my score to 6, and my confidence to 3. Ultimately, the proposed method for low-dimensional Bayesian optimization is a strong contribution. The performance of this method clearly outperforms competing baselines on both synthetic functions and hyperparameter tuning for ML problems. While these models are slower to train than simple GPs, the trade-off can be worth it in the significantly improved performance. Some comments that I leave with the authors for consideration: - I do not think the provided stopping criterion is correct. Specifically, it requires knowledge of the true optimum, which a BO algorithm does not have access to. - I am surprised by linear term in $R$ in the complexity; the authors suggest that increasing $R$ is expensive, however this does not seem to be particularly costly. Some more discussion on rank would be useful in the paper. - The training procedures for GPgrid and GP should be the same, and should therefore have the same training cost (they have different training cost in Table 3). - Additional baselines on the synthetic problems would make the paper more convincing. - The peak in the Damavandi function looked sharper to me than it truly is - I am now happy that a grid of 'good-enough' resolution would be able to easily include the optimum (or points close enough). I thank the authors for the time taken to provide responses to my concerns. I do hope that my comments have been taken and used to develop the paper.
Reply to Reviewer jRPp: Thank you for recognizing the contribution of this work
Thank you for recognizing the contribution of this work. We have carefully reviewed your comments and greatly appreciate your valuable feedback and suggestions. A brief response to the further comments: - 1. **Stopping criteria**: Yes, we agree with the reviewer. For realistic problems, a more suitable way could be to evaluate the acquisition function increment over two evaluation iterations to define the stopping criteria. We will cite relevant references in the paper to clarify this point. Thanks a lot for the comment. - 2. **Complexity relative to rank**: The model complexity is linearly proportional to rank because we sample the latent factors column by column. If one were to sample the entire vectorized factor, the complexity would be related to $R^3$. This complexity is discussed in paper [1], as shown in Table 1. We will clarify this point in the revised paper. Thanks for the comment. [1] Luttinen, Jaakko, and Alexander Ilin. "Variational Gaussian-process factor analysis for modeling spatio-temporal data." Advances in neural information processing systems 22 (2009). - 3. Yes, we agree that for GP fitting, the costs for GP and GPgrid methods are the same, scaling cubically with the number of observation points. We have updated the table above accordingly, and clarified that the table represents the complexity of the surrogate models in different BO methods. We will also include a table illustrating the costs for acquisition function computation. Thanks again for the detailed comment. - 4. We will consider more recent and relevant baselines. Thank you again for raising your score. Your comments help us improve this work a lot. Thank you very much.
Summary
The paper introduces a new surrogate model for Bayesian optimization, based on a functional tensor factorization. The approach discretizes the model to a pre-specified grid and uses MCMC sampling for inference. Bayesian optimization is carried out by selecting promising points from the pre-specified grid, as quantified by an acquisition function. The paper includes experimental results on synthetic functions as well as ML hyper-parameter tuning tasks.
Strengths
- Nice experimental results on the ML hyper-parameter tuning task. - Generally well written paper.
Weaknesses
- Optimization method appears to be restricted to an a-priori defined grid of possible candidate points. - I am concerned about reproducibility of the benchmarks, as the code submission not complete, e.g. appears to be missing implementations of basic functions like baseline GP fitting (`fitrgp`) and predicting (`predict`), the additive GP model mentioned in line 274 of the paper, the continuous optimization of the acquisition function. - Grid-based GP-UCB, GP-EI baselines are missing for ML tuning tasks (Figure 3). This is notable, because these experiments contain discrete variables, which requires a rounding operation if they are relaxed to a continuous space, as is done by the paper. This rounding operation is likely to degrade the performance of "continuous" GP-UCB and GP-EI, as it will can be prone to sampling "between" integers it has already seen, reducing its sample efficiency.
Questions
1) > BKTFrandom. BKTFrandom can be applied to functions with higher dimensions (e.g., D > 10). Did you try this out, and do you have any high-dimensional benchmark results to share? 2) > The predictive distribution for any entry f in the defined grid space conditioned on the observed dataset Dn can be obtained by the Monte Carlo approximation Can you clarify if the algorithm only works efficiently by amortizing MC inference to a fixed grid, and then optimizing the acquisition function on that grid? Is the implementation capable of numerical ("continuous) optimization of the acquisition function, and did you test this in any experiments? 3) > This allows us to accommodate observations that are not located in the grid space QD d=1 Sd. The detailed derivation of the sampling algorithm is given in Appendix A. It is not clear to me how Appendix A addresses the problem of off-grid inference, since e.g. it says: > "we define a binary tensor O with the same size of Y indicating the locations of the observation data". which only seems possible when restricted to a grid. 4) How is the acquisition function optimized for the methods that operate on the continuous space, like the non-grid GP EI, GP UCB, and additive GP methods (i.e. using any restarts, optimization algorithm, etc)? 5) Do you have any intuition why the results in Figure 2 are more mixed compared to the baseline, whereas BKTF seems to exhibit a clearer advantage on the the ML tuning tasks in Figure 3? In addition, please provide the results for the grid-based GP-EI and GP-UCB baselines, like for Figure 2, as this will determine if (part of) the performance gap can be explained by the continuous relaxation. Also, note that the the code for this experiment experiment is not included in the submission, in addition to the other missing code mentioned above. **Miscellaneous** 1) > perhaps the most popular kernel is the ARD (automatic relevance determination)— Squared-Exponential (SE) or Matérn kernel [4]. Although this specification has certain numerical advantages and can help automatically learn the importance of input variables, a key limitation is hat it implies/assumes that the underlying stochastic process is stationary and separable The SE kernel is equivalent to a separability (product kernel) assumption, but the Matérn kernel is not. ARD pertains to having separate lengthscales, which is not the same as an additive or multiplicative separability assumption. Please make the language more precise here. 2) > Another emerging solution is to use deep GP Typo: Either "deep GPs" or "a deep GP". 3) > we propose using Bayesian Kernelized Tensor Factorization (BKTF) as a flexible and adaptive surrogate model for BO in a D-dimensional Cartesian product space (i.e., grid) when D is relatively small (say D ≤ 10) The introduction also mentions that "the value of the covariance function between two random points quickly goes to zero with the increase of input dimensionality. These assumptions can be problematic". Since the introduction just stressed the limitation of canonical BO methods in higher dimensions, it appears important to stress that the method is not designed to mitigate this particular limitation, if it is really limited to D < 10. 4) > GP ARD: kernel function is stationary and separable Can you provide a reference for a body of work that equates GP ARD with a multiplicative separability assumption? In the GP literature, ARD is associated with having separate lengthscales for each input dimension (see e.g. Section 5.1 of https://gaussianprocess.org/gpml/chapters/RW5.pdf). This *happens* to imply separability for the SE / RBF kernel, but it does not imply separability in general, e.g. for the commonly used Matern kernel.
Rating
5
Confidence
4
Soundness
3
Presentation
3
Contribution
2
Limitations
> A limitation of BKTF is that we restrict BO to a grid search space in order to leverage tensor factorization; however, we believe that designing a compatible grid space based on prior knowledge is not a challenging task. An important limitation to highlight here once more that it goes from "not challenging" to prohibitively expensive as the dimension increases.
Thank you for your response. I think this work is interesting and has promise, but the rebuttal does not allow me to champion this paper. In particular, > (2) We have only submitted the code for the experiments in the Introduction section (see Figure 1 of the paper), where we do not compare with additive GP and other baselines. We will provide the full code after this work has been accepted. Why would you do that? The illustration in Figure 1 is neat, but the more important parts of evaluating a contribution are the experiments. > (2) We agree that the performance of continuous GP baselines could be impacted when discrete variables are included. In the revised paper, we will add a GP baseline model specifically designed for mixed spaces. This seems like something you could have readily addressed in the rebuttal. I am raising my score one point out of goodwill.
Reply to Reviewer MJcW
Thank you for your response and for raising your score. We are glad that you find this work promising. We explain to the further concerns below. - 1. Yes, we agree that the experiments in Section 5.1 and 5.2 are crucial for demonstrating the performance and contribution of our model. We have provided detailed illustration of the experiment setups and results in the paper and Appendices to ensure as much clarity and transparency as possible. All code will be provided later as well. - 2. About baselines used in the ML hyperparameter tuning tasks. We carefully reviewed the baselines we used. BO-TPE, in particular, can be considered as a GP baseline that specifically designed for problems with mixed continuous and categorical/discrete inputs. For other baselines, we are looking into additional literature to ensure appropriate comparisons are made. - 3. In addition, please note that we have updated the wall-clock time comparison table and the theoretical model complexity table in an additional response provided to all reviewers. Please refer to the updated tables for more detailed information about the time cost of the proposed model.
Summary
This paper proposes a method for Bayesian optimization where the prior is a low-rank sum of tensor products of GPs. An MCMC scheme is developed for approximate updating and UCB sampling. Several experiments show strong performance relative to baselines on artificial function optimization and ML hyperparam tuning.
Strengths
This is a seemingly new approach to BO with a more flexible type of prior, which shows good empirical performance.
Weaknesses
It's not clear what is new relative to previous BKTF papers [10,11], other than the scheme for using UQ in UCB sampling. The clearest potential drawbacks to the approach are the memory and compute costs. These should be reported for the experiments.
Questions
Would it be possible to store the means and variances just of the tensor factors ($g$) and use a mean field approximation to more cheaply approximate posterior? (6) and elsewhere: I think $\otimes$ would be more standard notation than $\circ$. 136: It would help to mention that each input dimension is standardized to $[0,1]$, because otherwise assuming equal length scales seems problematic.
Rating
6
Confidence
3
Soundness
3
Presentation
3
Contribution
3
Limitations
The paper motivates the approach in part because standard methods assume the generating process is stationary, but the BKTF is also stationary. It’s nonstationary only after conditioning on values of $g$.
Thank you to the authors for addressing my comments. I am happy with the reply regarding add-GP and non-stationary kernels. The discussion on rank is useful, and I appreciate the authors clarifying robustness. I further appreciate the time taken to reply to my questions. I am willing to raise my score. However, there are some discussion points below that I would like to be addressed, specifically: 1. Provide some clarifications on the time taken to run, and the time complexity of the models 2. Provide clarification on how the true optimum seemingly always lies in the grid **Wall-clock time, and Table 3** Thank you for providing the wall-clock times. I do **not** believe that slow model fitting harms the contribution of this work; there are many BO settings where the most significant cost is evaluation of the objective, and this method seems to minimize that. Nevertheless, it is important to provide clarity in this discussion in the paper. I note two issues here: 1. The paper does not provide stopping criteria for the algorithm. Therefore, one would use the entire BO budget with BKTF. So "time to optimum" is not necessarily the best metric to quote here (since a practitioner would likely evaluate many more samples as they do not know that the optimum has been found). This is especially the case if an adaptive grid is used, that changes the available set of points at each iteration. I would rather quote the time taken per iteration, then allow the reader to determine if that cost per iteration is worth the reduced number of iteration for their specific use case. 2. The authors quote the times for 2-dimensional problems. It would be better to quote these for some of the higher dimensional problems too. > The time for the entire MCMC sampling should be linearly proportional to the time provided in Table 3. It is important to note the scales of the terms. For example, when $n=10$, and the number of MCMC samples $K=1000$, then the linear term of $O(K\cdot n^3)$ is no longer negligible. The authors also note that the drawback of higher rank is in the computational cost, but the rank is not included in the BKTF complexity term. Moreover, I also note that complexity given for BKTFrandom in Table 3, $O(n_\text{random}^3)$, is significantly bigger than $O(n^3)$, since $n_\text{random}=20000$. To be clear, in the complexity term for BKTF is $n$ the number of observations, or the number of candidates? If the latter, then this is significantly worse than the GP methods. I therefore strongly disagree with the comment made in Appendix F, > theoretically the proposed model BKTF/BKTFrandom has the lowest computational complexity **Grid** There has been discussion across the reviewers about the grid, and the specification of the grid. This is an issue when the optimum does not lie on the grid. In the synthetic functions, do the authors always ensure that the global optimum lies in the grid? For example, a random grid would be **highly unlikely** to include the optimum of the Damavandi synthetic function. Otherwise, we should expect to see that the non-grid based methods (eg $\alpha_{EI}$) should eventually reach a lower optimum than any of the grid based methods. The authors refer to the possibility that the grid space being adaptive. Is this done during the experiments in the paper? **Dimensionality** The model shows strong performance on low dimensional problems. This is still a strong contribution, and many BO problems lie in this dimensionality of problem. I would appreciate if the authors could make clearer (for example, in the abstract) that this is the specific regime in which this model performs, in the final paper. (Minor) SAASBO is intended for high dimensional problems (see the paper, where the problems are $D>100$). I would argue that 6 and 10 dimensional are still low-dimensional, not "relatively high-dimensional". (Minor) It is possibly to use the approach for BKTFrandom with the other baselines as well (For GPgrid $\alpha_{EI}$, only evaluate on a subset of the grid).
Reply to Reviewer jRPp
Thank you for your timely and thoughtful response. We are pleased to have addressed several of your concerns and appreciate your positive feedback on our work. Regarding the further concerns on **time complexity** and **grid** assumptions of the proposed model, we respond to each comment below. - **Wall-clock time and Table 3** 1. Regarding **Wall-clock time**: (1) About `stopping criteria`. For optimization on benchmark test functions, we define **_regret_** as the absolute error between the true global optimum and the current estimated global optimum, as detailed in lines 284-287 of Section 5.1 Results in the paper: "To compare the optimization performance of different models on the benchmark functions,we define 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|$, as the regret, and examine how regret varies with the number of function evaluations." One way to define the `stopping criterion` is based on the regret; specifically, it can be defined as the _regret_ being smaller than a threshold $\epsilon$, e.g., $\left|f^{\star}−\hat{f}^{\star}\right|<\epsilon=1e-6$. In the paper, for test function optimization, we compute the average costs (budgets) or time for evaluation as the number of evaluations or the time taken to achieve _regret_=0, respectively. (Note that we provide the final _regret_ and average costs of evaluations for experiments on test functions in **Table 2, Appendix E.2**.) Therefore, the `stopping criterion` can be considered as setting $\epsilon=0$. However, we agree that using a relatively more relaxed threshold, such as $1e-6$, would more fairly quantify the performance of continuous baseline models. We have made this adjustment in the updated wall-clock table. Additionally, there are also other ways to define stopping criteria, such as based on the incremental change between acquisition functions over two iterations. We will clarify the `stopping criterion` in the revised paper. On the other hand, a careful examination of the results in Table 2 indicates that the evaluation costs may not significantly change regardless of the stopping criteria, particularly for the baselines that fail to find the global solution. These baselines still do not achieve the global optimum within the budget, as evidenced by the final _regret_ being much larger than $1e-6$. For example, all GP-based baselines for Damavandi function show final regrets of 2, indicating that these methods are stuck at local optima. The updated wall-clock table below also reflects this point. In addition, for ML tuning tasks, one can use as many evaluations as their budget allows to find the final results. (2) We have **updated the wall-clock comparison table** (see below), where "-" still denotes that the model does not find the global solution within the given budget. Specifically, as the reviewer suggested, we provide the **average time taken per (evaluation) iteration** and **average costs of evaluation** (using a stopping criterion of _regret_<1e-6) to more clearly illustrate the wall-clock time costs of the proposed method. In addition, we include the results for **all the test functions** discussed in Section 5.1. We also added a comment providing the updated table as an additional response to all reviewers. Thanks a lot for the comment. |$f(\boldsymbol{x})~(D)$|GP$\alpha_{\text{EI}}$|GP$\alpha_{\text{UCB}}$|GPgrid$\alpha_{\text{EI}}$|GPgrid$\alpha_{\text{UCB}}$|additive GP|BKTFrandom|BKTF|deepGP| | ----------------------- | ----------------------- | ------------------------ | --- | --- | --- | --- | --- | --- | | Branin (2)|0.35s/37iter|0.37s/36iter|0.33s/23iter|0.28s/36iter|0.82s/100iter|0.19s/47iter|0.33s/4iter|1.29s/47iter| | Damavandi (2)|0.82s/-|0.29s/-|0.19s/-|0.14s/-|1.96s/-|2.68s/48iter|2.64s/5iter|3.78s/-| | Schaffer (2)|0.34s/36iter|0.34s/44iter|0.17s/>50iter|0.05s/>50iter|1.53s/43iter|0.20s/54iter|0.19s/22iter|-| | Griewank (3)|0.28s/>100iter|0.57s/>100iter|0.18s/>100iter|0.11s/>100iter|1.75s/>100iter|0.37s/47iter|0.41s/43iter|-| | Griewank (4)|0.25s/>100iter|0.59s/>100iter|0.40s/>100iter|0.26s/>100iter|2.03s/>100iter|0.74s/87iter|0.85s/68iter|-| | Hartmann (6)|0.52s/>100iter|0.56s/>100iter|5.62s/>100iter|5.55s/>100iter|3.33s>100iter|9.26s/154iter|39.05s/60iter|-| | Griewank (10)|0.44s/>200iter|0.57s/>200iter|-|-|3.97s/>150iter|6.63s/124iter|-|-| 2. Regarding the theoretical time cost in **Table 3**: (1) Yes, we agree that both the number of MCMC iterations $K$ and the rank $R$ should be considered in the computational complexity. Additionally, as reminded by the reviewer, the number of optimization iterations required to fit/solve the GP model parameters/hyperparameters should also be considered for GP baseline models.
Thank you for the detailed response. However, the rebuttal does not allow me to change my rating significantly. In particular, ``` We'd like to clarify that the two works you mentioned, e.g., KISS-GP, are not Bayesian models. They do not provide predictive/posterior distributions for unobserved data, cannot offer uncertainty quantification (UQ), and cannot be used for Bayesian optimization (BO). These deterministic methods are not within the scope of this work and should not be discussed. ``` No. There are bayesian methods that provide predictive posterior. Also, the comparison to SOTA methods is not resolved, which hopefully will be resolved in the future.
Statement of Computational Cost
Dear reviewers, We would like to provide two tables to clarify/restate the **computational cost** of the proposed BKTF surrogate, addressing both the **theoretical (model complexity)** and **practical (wall-clock time cost)** aspects, respectively. - **Theoretical model complexity**: Following the reviewer's comments, we corrected Table 3 in the paper as provided below, where $n$ is the number of observed data, $m_d$ is the number of interpolation points, $J$ and $J'$ represent the optimization iterations taken for fitting a standard GP and an additive GP model, respectively, $K$ denotes the number of MCMC iterations, and $R$ is the rank. | Model | Complexity | | -------- | -------- | | GP $\alpha_{\text{EI}}$ | $\mathcal{O}\left(Jn^3\right)$ | | GP $\alpha_{\text{UCB}}$ | $\mathcal{O}\left(Jn^3\right)$ | |GPgrid $\alpha_{\text{EI}}$|$\mathcal{O}\left(Jn^3\right)$| |GPgrid $\alpha_{\text{UCB}}$|$\mathcal{O}\left(Jn^3\right)$| | additive GP | $\mathcal{O}\left(J'n^3\right)$ | | BKTF/BKTFrandom | $\min\left( \mathcal{O}\left(Kn^3\right),\mathcal{O}\left(KR\sum_{d=1}^D m_d^3\right) \right)$ | | deepGP | $\mathcal{O}\left(K\left(\prod_{d=1}^D m_d\right)^3\right)$ | Please note that this table illustrates the complexity of surrogate models in different BO methods. We will update Table 3 in the revised paper and also include a table comparing the costs of acquisition function computation. - **Practical wall-clock time cost**: We provide a more comprehensive wall-clock comparison table as shown below, where "-" still denotes that the model does not find the global solution within the given budget. The table now includes: (1) **average time taken per evaluation iteration**; (2) **average costs of evaluation** (using a stopping criterion of regret<1e-6); and we provide the results for **all test functions** in Section 5.1. These additions aim to illustrate the wall-clock time costs of the proposed method more clearly. |$f(\boldsymbol{x})~(D)$|GP$\alpha_{\text{EI}}$|GP$\alpha_{\text{UCB}}$|GPgrid$\alpha_{\text{EI}}$|GPgrid$\alpha_{\text{UCB}}$|additive GP|BKTFrandom|BKTF|deepGP| | ----------------------- | ----------------------- | ------------------------ | --- | --- | --- | --- | --- | --- | | Branin (2)|0.35s/37iter|0.37s/36iter|0.33s/23iter|0.28s/36iter|0.82s/100iter|0.19s/47iter|0.33s/4iter|1.29s/47iter| | Damavandi (2)|0.82s/-|0.29s/-|0.19s/-|0.14s/-|1.96s/-|2.68s/48iter|2.64s/5iter|3.78s/-| | Schaffer (2)|0.34s/36iter|0.34s/44iter|0.17s/>50iter|0.05s/>50iter|1.53s/43iter|0.20s/54iter|0.19s/22iter|-| | Griewank (3)|0.28s/>100iter|0.57s/>100iter|0.18s/>100iter|0.11s/>100iter|1.75s/>100iter|0.37s/47iter|0.41s/43iter|-| | Griewank (4)|0.25s/>100iter|0.59s/>100iter|0.40s/>100iter|0.26s/>100iter|2.03s/>100iter|0.74s/87iter|0.85s/68iter|-| | Hartmann (6)|0.52s/>100iter|0.56s/>100iter|5.62s/>100iter|5.55s/>100iter|3.33s>100iter|9.26s/154iter|39.05s/60iter|-| | Griewank (10)|0.44s/>200iter|0.57s/>200iter|-|-|3.97s/>150iter|6.63s/124iter|-|-| As demonstrated, the proposed BKTF surrogate successfully finds the global solution within the given budget in all cases, utilizing acceptable optimization time. This indicates the advantage of the proposed surrogate in optimization on nonstationary complex functions with a limited budget and observations.
Thanks for the clear replies. I continue to support acceptance.
Decision
Reject