Standard Gaussian Process is All You Need for High-Dimensional Bayesian Optimization

To re-state it as a question: Are you claiming that practically relevant optimization problems do not have an exponentially growing (in $d$) number of local optima? Or perhaps less explicitly, that practically relevant problems have some other property that means they can be globally optimized with a sample complexity that grows sub-exponentially with $d$? Which part of your work supports this claim? For example, would you say that the problems listed in Figure 4 are a good surrogate for "practically relevant" problems? How so?

R1: Thanks for the insightful questions. We do agree that even there is no training optimization issues, to realistically find the global optimum, the objective function itself must be nice enough and additional assumptions must be made, such as the one you mentioned: "practically relevant optimization problems do not have an exponentially growing (in $d$) number of local optima". However, similar to other high-dimensional BO studies, we take a more pragmatic view. That is, we do not require BO to always find the global optimum in practice. In other words, we are open to applying BO to a broad range of problems, even when the global optimum may be challenging to reach. This perspective is common in gradient-based optimization, where, for many applications (e.g., neural network training), a good local optimum suffices for practical use. We believe this expectation broadens BO's practical applicability.

In our experiments, we actually have tested both the nice and tough problems regarding the number of local minima. Among the synthetic benchmarks, Hartmann6 and Rosenbrock are problems that do not have exponentially growing local optima. Hartmann6 has a fixed dimensionality of 6, while Rosenbrock is unimodal. One can see from Figure 3 and 5, BO successfully converges to the global maximum (3.3237 for Hartmann6 and 0 for Rosenbrock).

In contrast, Ackley and StybTang exhibit an exponential growth in local optima with the input dimension $d$. For Ackely, the proliferation of local optima stems from the term $\exp(\frac{1}{d}\sum_{i=1}^d \cos(2\pi x_i))$, and for StybTang, from each summand $(x_i-c_i)^4-16(x_i-c_i)^2 + 5(x_i-c_i)$; see Section B.1 in page 14. Since both objective functions are factorized across dimensions, and each dimension $i$ independently contributes multiple local optima for $x_i$. The combination leads to exponential growth in the total number of local optima with increasing $d$. These functions are more challenging to optimize. The global maximum is 0 for Ackely and 7833.198 for StybTang. As shown in Figure 3 and 4, after 500 BO iterations, the best function values found are still far away from the global maxima.

mini issues, easy to fix: Section 2, line 091... Proposition 4.3 ...Eq. 6: the parentheses...

R2: We greatly appreciate you pointing out these issues. We will address and correct them in our manuscript.

I'm not convinced by the experiments design in Section 4 ... To my understanding neither the $L_2$ difference nor the average gradient norm at the first training step reveals the presence of vanishing gradient throughout the training process ... Vanishing gradient could still occur in the following steps.... The values of length-scale initialization ($\ell_0$) is not robustly tested...

R2: There appears to be a misunderstanding about the motivation and design of our experiments. To clarify, the theoretical analysis in Section 4 is aimed at demonstrating that the gradient vanishing can occur at the very first step of training due to an inappropriate length-scale initialization. When gradient vanishing occurs at this initial step, it persists throughout training, as each subsequent step continues to compute gradients at the same initial values, preventing any updates to the length scales and ultimately causing training failure. Thus, we specifically use (1) the gradient norm at the first step to examine whether gradient vanishing occurs initially, and (2) the relative $L_2$ difference between the length-scales before and after training to confirm that, once gradient vanishing happens at the first step, it persists across all training iterations, resulting in no change after training. Our paper jointly interprets the results of (1) and (2) to verify our theory, rather than treating these metrics independently.

In other words, our experiments are not designed to assess whether gradient vanishing might occur later in the training process. In such cases, the length-scales would already be effectively updated (before any gradient vanishing), meaning training would not necessarily fail.

Second, we employed both commonly used length-scale initialization values and our proposed values for verification. The results sufficiently demonstrate that smaller values are more prone to gradient vanishing and training failure, while larger initializations handle higher dimensions robustly (our results are averaged over 20 runs); see Lines 245-260.

Moreover, with $\ell_0 = 1.0$, there has already been no gradient vanishing observed for the Matern kernel across all the test dimensions. This suggests that even larger values, such as $\ell_0 = 2$ or $3$, would also prevent vanishing gradients across all dimensions, making it unnecessary to duplicate results with larger initializations. While further testing with additional values is feasible, we believe our current test range is sufficient to verify our theoretical analysis.

Questions regarding $a$ in training and analysis

R3: Thank you for these questions. In our experiments, the initialization value of $a$ was set to $\text{SoftPlus}(0) = 0.693$, which is the default choice in the BoTorch library. We will clarify this in the paper. There is a typo in the gradient equation for the Matern kernel, where $a$ was omitted; $a$ should be a factor. We will correct this. Lastly, yes, the amplitude $a$ is treated as a constant during the theoretical analysis.

In Figure 2, the relative length-scale differences for both GP-Matern($\ell_0=0.693$) and GP-SE($\ell_0=\sqrt{d}$) are flat lines. Since $\ell_0$ is initialized at the same value independent of the number of iteration, does this mean the value of length-scale is all the same after training even for different GP (since the GP at each iteration is modeled on different data sets)?

R4: No. The figure just shows that the overall change of the length-scale values, $L_2=\frac{|| \ell_1 - \ell_0||}{||\ell_0||}$, is close across the BO iterations, but it does not mean the actual learned values are the same. For example, starting from $\mathbf{\ell}_0 = [1,1,1]$, consider two sets of the learned length-scales $\mathbf{\ell}_1 = [0.7, 0.6, 0.3]$ and $\mathbf{\ell}_2 = [0.5, 0.505, 0.505]$. They are obviously different, but the relative $L_2$ difference from the initialization is almost the same: $||\mathbf{\ell}_1 - \mathbf{\ell}_0|| = 0.86023$ and $||\mathbf{\ell}_2 - \mathbf{\ell}_0|| = 0.86026$.

We thank the reviewer for their time and detailed feedback. We are concerned, however, that certain key aspects of our work --- including its primary contributions, motivation, results, and the interpretation of our gradient vanishing experiment --- may not have been fully understood. We provide the following clarifications to address these points in detail.

This paper mainly focuses on evaluating existing method. Vanishing gradient problem has been widely discussed in gradient-based training. I think the novelty of this submission only comes from the proposed initialization method ...

R1: We respectfully disagree.

First, our paper does not "mainly focus on evaluating existing methods." As clearly summarized in the Introduction (Lines 58-81), this work includes three main contributions: Empirical Results, Theory, and Simple Robust Initialization. While we conducted extensive empirical evaluations, these form just one part of our contributions. The theoretical analysis and its numerical verification are central to our study, occupying significant space: the entirety of Section 4 (pages 4-6), two tables (Tables 1-2), two large figures (Figures 1 and 2), and all of Section A in the Appendix. We are uncertain as to why this substantial portion of the work was overlooked by the reviewer.

Second, we were surprised to see our theoretical analysis is dismissed as lacking novelty or contribution on the grounds that "the vanishing gradient problem has been widely discussed in gradient-based training". This critique seems overly reductive and sloppy. By that logic, the vast body of contemporary neural network research would also lack novelty simply because neural networks have been widely discussed since the 1980s.

To clarify, our contribution is not the concept of "gradient vanishing" itself. Rather, we identified why standard BO with a SE kernel often fails in high-dimensional problems --- specifically, because training fails due to gradient vanishing right from the initial training step. This discovery is non-trivial, as BO is a complex, iterative process involving interleaved GP training and acquisition function optimization. We not only identify that the vanishing gradient arises from an inappropriate yet commonly used length-scale initialization, but also provide a probabilistic bound that characterizes the likelihood of gradient vanishing, which grows rapidly with the dimensionality, $d$. To our knowledge, the identification of this failure mode in BO and the accompanying quantitative analysis are the first in this area, and therefore represent novel contributions of our work.

Thank you for the clarifications. I was a bit disappointed previously because when the authors state they are paying attention to the issue of gradient vanishing, I was expecting a more general approach to resolve this problem for the entire training process. I was also misled by the typo made in the authors' comment for their experiment results. Now after the authors' further clarification on the conclusion of the experiments and the scope of the problem they are studying, I acknowledge the novelty and contribution of this work and increase my score.

I would like to make two more comments:

1. The phrase $**Runtime performance**$ (line 448) is misleading. Readers will expect to see actual running time reported. The authors might consider using a different word.
2. The authors used a very general title for their submission which gives readers impression that the challenge in high-dimensional BO has been solved by standard GP (also one of the reasons why I was expecting a lot previously $**:(**$ )

We thank you for your support, and the many insightful and constructive comments.

Insufficient replicates in the numerical experiments....Really, there should probably be closer to 50 replicates or more.

R1: Thank you for the great suggestion; we agree. As recommended, we will add more replicates --- around 40 additional runs or more --- to further reduce variance and strengthen our findings.

cite documents about "Criticism of the squared exponential and advocacy of the Matern kernel" and "trying to locate other discussion of this topic."

R2: Thank you for the helpful suggestions! We will add the provided reference along with additional relevant literature on this topic. We will also expand our discussion of this point.

Why not use a Matern kernel with the proposed initialization in your numerical experiments?

R3: Great question. We empirically found when Matern kernel with commonly used initialization does not exhibit gradient vanishing and training failure, its performance is close to that achieved with our robust initialization. Since the Matern kernel with commonly-used initialization performs well across all benchmarks in the paper (without gradient vanishing), adding additional results with our new initialization would make the figures overly cluttered --- particularly in Figures 3, 4, and 5, where numerous baselines have already been included --- and would reduce readability.

Here, we have included a comparison of the Matern kernel using both initializations on the Rover, DNA, SVM, and Hartmann6 benchmarks. The results are shown in final and running , demonstrating that the performance is indeed very similar between the two initializations. We will clarify this point in the paper.

We have also conducted an additional test on a new benchmark, Humanoid, with a dimensionality of $d = 6,392$. Please refer to our response to @Reviewer M8nK for the results. In such high dimensions, we observed that commonly used initialization leads to gradient vanishing and training failure even for the Matern kernel. However, our robust initialization approach dramatically enhances performance, achieving state-of-the-art results.

Under the assumption of uniform sampling, does any of the procedures of DiceKriging, hetGP or GPVecchia satisfy your robust initialization?

R4: Thanks for providing these excellent GP learning packages. We appreciate the opportunity to compare our initialization approach with those used in DiceKriging, hetGP, and GPVecchia.

For DiceKriging, the initializations are independently sampled across input dimensions, with values ranging between $1 \times 10^{-10}$ and twice the maximum possible value in the domain (in our case, this is $[1 \times 10^{-10}, 2]$). Therefore, these initialization values are independent of the input dimension $d$ and do not align with our approach.

In hetGP and GPVecchia, the initialization strategies rely on the training data, making them data-dependent. However, in expectation, their initializations can match our approach. For hetGP, the initialization can be set to the 50th percentile of the observed pairwise distances. This provides a stochastic estimate of the median pairwise distance, which is on a similar scale to the mean pairwise distance. In our analysis, we assume a uniform data distribution within $[0, 1]^d$, where the mean pairwise distance between data points is approximately $\sqrt{d}/6$. Thus, the hetGP initialization is, on average, consistent with our robust initialization.

Similarly, in GPVecchia, the initialization is a proportion of the mean distance between training points, which provides a stochastic estimate of the mean distance scaled by a constant. This, too, aligns with our robust initialization strategy in expectation.

We will add a discussion of these methods and the similarities and differences in our paper.

Have you done experiments on problems with even higher dimensions, e.g. in the thousands?

R1: While our benchmarks are drawn from widely used and recognized sources in high-dimensional optimization literature, we agree that testing in even higher dimensions, such as in the thousands, is valuable. To address this, we conducted an additional test on the Humanoid benchmark, a reinforcement learning task with a dimensionality of 6,392, as used in recent work by Hvarfner et al., (2024). We ran standard BO with both Matern and RBF kernels using MLE estimation and ADAM optimization over 1,000 BO iterations.

The results are available here. In this high-dimensional setting, we observed that commonly used initializations lead to gradient vanishing and training failure for both kernels. With our robust initialization method, however, performance for both kernels improved dramatically, achieving results close to VBO (Hvarfner et al., 2024), which employs a log-normal prior for training.

We will supplement our paper with these results on the Humanoid benchmark.

Like pG21 I wish to thank the authors for the interesting work and the reviewers for their insights.

@pG21: Thank you for your comments. Let me expand on some points and clarify further.

1. That is indeed an excellent point --- thank you for highlighting it, I was not aware of this. With that in mind, there is certainly a case to be made for considering it concurrent work.

I think you're absolutely right that the notion of "simpler" is important here, however, it can be interpreted in different ways depending on context.

I agree that it can sometimes be simpler for a practitioner to implement it this way.

However, the approach of "encouraging" a solution via its initialization primarily applies when employing local optimization methods for hyperparameters, such as gradient descent or Adam. It does not generalize well to global optimization algorithms. Defining an appropriate lengthscale prior is much more general, and in combination with other things, can often make things simpler rather than more difficult. This approach explicitly embeds the desired properties into the optimization objective for the lengthscale, making it adaptable to a wide range of optimization methods (local or global).

When relying solely on initialization to encourage a solution, you're implicitly assuming that the algorithm's movement away from the initial point will be limited. This assumption can become problematic. If the optimization algorithm is highly effective and explores globally optimal solutions, it might counteract the initialization's influence, potentially returning to less desirable lengthscales (which can be better according to the objective/loss function). This is particularly true for global optimization algorithms (which are popular and useful in these settings), where the influence of initialisation tend to be much weaker. To mitigate this, you might need to constrain the algorithm artificially (e.g., limiting exploration or stopping optimization early), which can add complexity rather than reduce it.

By contrast, modifying the objective function for the lengthscale optimization provides a clear, unambiguous goal, and is easily done via the hyperparameter prior. This approach aligns the optimization process with the desired outcome, which in turn can lead to a more robust system and a simpler experience for the practitioner.

Moreover, in BO it is common and often advantageous to rely on Bayesian inference, to consider a (data agreement-weighted) range of hypotheses for the hyperparameters, instead of relying on point estimates [1]. In this context, a well-defined prior is far more practical than relying on implicit biases introduced through initialization. The prior integrates seamlessly with inference procedures, avoiding the need for algorithmic workarounds.

"What if I already had a different prior in mind that I wanted to use?" This is a great question! The prior is simply a density function. Importantly, it doesn't need to be normalized for optimization or Bayesian inference, as the normalization constant often cancels out during these processes. This flexibility makes it easy to combine or adjust priors to reflect different assumptions or preferences.

My point here really is that encouraging a solution via its initialization is changing the prior, it just means doing so in an implicit way, where you also need to consider how the optimization/inference algorithm works. Sometimes this is a practical, pragmatic approach which makes things simpler to implement, but it has its limitations, and can make things harder instead.

In my opinion, now considering this concurrent work, I think the paper can be made ready for publication by merely smaller edits with regards to the other paper. In particular, under "Alternative method" (line 353-360): (1) making it clear that it is concurrent work while highlighting differences in focus, and (2) remove/lessen the argument about it being "simpler", which argued above, is needlessly ambiguous, and very context specific.

[1] Snoek, Jasper, Hugo Larochelle, and Ryan P. Adams. "Practical bayesian optimization of machine learning algorithms." Advances in neural information processing systems 25 (2012).

Thank you for your detailed and constructive comments. We would especially like to thank @Review pG21 for clarifying that our work is concurrent with (Hvarfner et al., 2024) and highlighting the differences between our approach and theirs. We will incorporate these clarifications into our paper.

The difference from (Hvarfner et al 2024).

R1: We agree that both (Hvarfner et al., 2024) and our work demonstrate that standard (or vanilla) Bayesian Optimization (BO) can perform effectively in high-dimensional optimization. Below, we provide a detailed summary of our unique contributions:

1. We identify that the failure of standard BO, particularly with the SE kernel, stems from training failure caused by gradient vanishing. The gradient vanishing arises from commonly used but inappropriate initialization.

2. We provide a mathematical analysis of the gradient vanishing issue, presenting several bounds, including a probability bound, and characterize its relationship with the input dimension $d$. Numerical verification further supports our analysis.

3. We find that the Matern kernel is a more robust choice than the SE kernel for high-dimensional optimization.

4. We propose a new initialization approach that aligns conceptually with (Hvarfner et al., 2024) but differs in implementation; specifically, we select a new initialization rather than constructing a prior distribution.

We will emphasize these distinctions in our paper to clarify the differences between our work and (Hvarfner et al., 2024).

Regarding "Simpler"

R2: Thank you for the detailed discussion; we agree with your points. We would like to clarify that we do not claim or imply that our initialization strategy is superior to constructing a prior. By "simpler," we mean that (1) it is more straightforward to implement, and (2) it avoids the need to design a distribution, which is often a nontrivial intellectual process. We are happy to remove the term "simpler" to prevent any potential confusion or debate.