Lightspeed Black-box Bayesian Optimization via Local Score Matching

Venue: NeurIPS_cc_2024_Workshop_BDU
Published: 2024-09-06 13:47:39

Abstract: Bayesian Optimization (BO) is a powerful tool for tackling optimization problems involving limited black-box function evaluations. However, it suffers from high computational complexity and struggles to scale efficiently on high-dimensional problems when fitting a Gaussian process surrogate model. We address these issues by proposing a fast acquisition function maximization procedure. We leverage the fact that Probability Improvement (PI) acquisition function is a likelihood function whose score can be estimated through a simple linear regression problem called local score matching. This enables fast gradient-based optimization of the acquisition function, and a competitive BO procedure which performs similarly to that of computationally expensive neural networks.

Main paper: 
 introduction: 
 Black-box optimization [1, 2] seeks to identify the optimum (maximum in this paper) of a function g(x) whose closed-form expression and gradient information are unknown, with as little computational resources as possible. Bayesian optimization (BO, [6] ) is particularly effective for this task.It utilizes a probabilistic surrogate, typically a Gaussian Process (GP, [19] ), to sequentially select new evaluation points based on the mean function and quantify uncertainty through the covariance function. Although BO can achieve relatively good accuracy with only a few function evaluations, the cost of each new proposal evaluation scales as O(n 3 ) with the number of function evaluations n, which becomes the dominant factor for overall cost in the optimization process.To reduce computational expenses, a new paradigm called Bayesian optimization via density ratio estimation (BORE, [4, 16] ) was introduced. BORE reformulates the improvement-based acquisition function (e.g. probability improvement, PI [10] , expected improvement, EI [11] ) as a problem of estimating the density ratio [15, 21] between two distributions of x conditioned on whether the corresponding function value y exceeds a certain threshold τ . Song et al. [13] further generalize BORE to likelihood-free Bayesian optimization (LFBO), where the acquisition function can be constructed in terms of more complex utility functions through variational representation [12] . This approach significantly reduces the computational complexity from O(n 3 ) to O(n). However, the acquisition function itself in both BORE and LFBO still remains demanding to optimize due to:• The involvement of a neural network typically requiring cumbersome numerical optimization.• Potential overfitting in the density ratio estimate, particularly in high dimensions.Motivated by these challenges, we present a new approach that maximizes the acquisition function by gradient ascent. The key observation is that the PI acquisition function is a likelihood function whose gradient could be estimated using score matching. Under mild regularity conditions, we utilize a simple regression model to learn the score without training any neural network. Moreover, our method maintains a computational complexity at O(n). 
backgroud: 
 Consider the global optimization problem for a black-box function f over a compact search space X ⊂ R d :EQUATIONwhere the expression and gradient of the objective function f : X → R are unknown. We can only access f through a set of noisy observationsD N = {(x n , y n )} N n=1, where the evaluations y = g(f (x); ϵ) are corrupted by noise ϵ. In many cases, we also have a computational budget B.By specifying a probabilistic surrogate model, we define the acquisition function α(x; D N , τ ) as the expected value of the utility function u(x, y, τ ) under its posterior predictive p(y|x, D N ):EQUATIONwhere τ is a hyperparameter that balances exploration and exploitation. BO selects the candidate solutions by maximizing acquisition x t+1 = arg max x α(x; τ ). In GP-based BO [5] , most acquisition functions have closed-form solutions that are easy to optimize [20] .We focus on a particular utility function, the probability improvement(PI): u(x, y, τ ) := I(y -τ ≥ 0), where τ is a threshold of the function value. By definition, the acquisition function isα(x; τ ) = E p(y|x,D N ) [I(y -τ ≥ 0)] = p(y ≥ τ |x).This acquisition function measures the probability that a new candidate point x yields an observation y greater than the threshold τ . τ is typically selected as the current maximum observed function value: τ = max n y n . 
methodology: 
 First, let z denote a binary variable:EQUATIONThen, the PI acquisition function can be rewritten as α(x; τ ) = p(z = 1|x). We propose to directly maximize the PI acquisition function by performing gradient ascent :EQUATIONWe modelEQUATIONWe leverage score matching techniques [8, 18] to learn the score of the acquisition function∇ x log p(z = 1|x) through a least squares regression:EQUATIONwhere x k-1 is from the previous value of x in the gradient ascent algorithm. The joint probabilityis factorized as p(z, x|x k-1 ) = p(z|x)q(x|x k-1 ) 1 , where q(x|x k-1 ) is a proposal distribution. It samples candidate x around x k-1 .Since the target score ∇ x log p(z = 1|x) in the objective ( 6) is unknown, we cannot directly minimize J(β). The following theorem provides a closed-form solution of ( 6) and forms the basis of our algorithm:Algorithm 1 Maximizing Acquisition Function via Gradient Ascent Require: Generator y = g(f (x); ϵ), hyperparameters σ, η, threshold τ , input x k=0 steps of gradient ascent K, number of samples M 1: for k = 1, ..., K do ▷ Perform gradient ascent K steps 2:while m ≤ M do ▷ Survey local landscape3: sample x (m) ∼ q(x|x k-1 ) 4:evaluate y (m) ← g(f (x); ϵ)5: z (m) ← I[y (m) ≥ τ ] 6: m ← m + 1 7:end while 8:D ′ k ← {(x (m) , y (m) )} M m=1 ▷ Collect sample 9:β 1 ← Monte Carlo approximate according to (9) ▷ Local score estimation 10:x k ← x k-1 + η β 1 ▷ Update location 11: end for 12: Output: the final candidate x K , used sample D ′ 1:K Theorem 1. Assuming that the probability p(z|x) and the proposal distribution q(x|x k-1 ) are differentiable. The objective (6) has the optimal solution:EQUATIONwhere the second equality follows from integration by parts under mild assumptions.Note that the minimizer β * 1 in Eq. ( 7) is a conditional expectation of the desired score∇ x log p(z = 1|x)| x=x k-1, implying that it is a biased estimator. The following result ensures that it is asymptotically unbiased given an appropriate choice of the proposal distribution.Corollary 1. Let the proposal distribution q(x|x k-1 ) be a normal distributionN (x|x k-1 , σ 2 I d ).Thenβ * 1 = E p(x|z=1,x k-1 ) σ -2 (x -x k-1 )andEQUATIONCorollary 1 shows the solution to the objective ( 6) is tractable and, when σ goes to zero, provides an asymptotically unbiased estimate of the desired score. The proofs of Theorem 1 and Corollary 1 can be found in the appendix.We summarize our gradient ascent algorithm in Algorithm 1 and defer the full BO algorithm to Appendix B. We now detail the computation of β 1 using D ′ k , i.e., paired function evaluations and their input values. According to Corollary 1,β * 1 = p(x|z = 1, x k-1 ) σ -2 (x -x k-1 ) dx.Once the threshold τ is determined, we can approximate the above integral using a Monte Carlo estimator:EQUATIONwhere x (m),1 is short for x (m) |z (m) = 1. If z (m) = 0 for all 1 ≤ m ≤ M we simply set β 1 = 0.Intuitively, β 1 can be interpreted as an average of (x -x k-1 )/σ 2 using local information that is {(x (m) , z (m) )} M m=1 sampled within a neighborhood of x k-1 with radius determined by σ. 
experiments: 
 Since we estimate the gradient of the acquisition function, we can resort to first-order optimizers such as ADAM [9] and RMSProp [17] to maximize the acquisition function. Note that Corollary 1 suggests that σ should be small. However, since σ controls the local survey region as well, we propose a linear annealing schedule to balance this trade-off.σ 2 t = σ 2 0 • max 0, 1 - t -0.1 TWe quantitatively evaluate the performance of our method by benchmarking against some standard black-box optimization problems, fixing a function evaluation budget of 250. We report the immediate regret as a metric, which measures the distance of the current best function evaluation from the optimum function evaluation. We compare against the neural-network based LFBO approach, which despite using a different acquisition function, is a similar approach to our work that does not require inference of a surrogate model. We also provide the results using a random-search algorithm as a reference. From these results, we observe that our local score matching method (SM) is broadly competitive with the LFBO method. It outperforms LFBO in higher dimensions in the Rosenbrock experiment when the neural network estimation in LFBO is likely to struggle. 
limitations and future works: 
 As the gradient information is myopic, our algorithm may not achieve the global optimum of the acquisition function. In addition, our method is restricted to PI acquisition function. An interesting avenue for future work is to extend our method to a broad class of acquisition functions. We also notice that our algorithm is very similar to Covariance Matrix Adaptation -Evolutionary Strategy (CMA-ES, [7] ). Further investigations into the relationship between these two methods could be a promising direction. Now we set the proposal distribution q(x|x k-1 ) = N (x|x k-1 , σ 2 I d ) as Gaussian, whose limit distribution is is the dirac delta function centered on x k-1 as σ → 0. Under the same assumption in Proof of Theorem 1, we now have:EQUATIONwhen σ → 0. D t ← D t-1 ∪ K k=1 D ′ k ∪ (x t , y t ) ▷ Update dataset 9:t ← t + 1 10: end while Given T, K, M , the total number of evaluations required in our algorithm is B = T × (K × M + 1).Theoretically, x t-1 can be initialized to any location x 1:N in D N . Experimentally, we set it as the x of which corresponding evaluation is max y, y ∈ D N by the greedy heuristic.For all algorithms, we initialized the dataset with 4 initial points, which we do not count in our function evaluations, as is necessary for the LFBO algorithm. This initial dataset is the same across all algorithms. For the random search algorithm, the optimal point and function evaluation is initialized from this dataset.For our local score matching algorithm, the Adam optimizer is used for the maximizing of the acquisition function. We used a fixed step-size of 5 • 10 -1 throughout all experiments. 3 We used a budget of T = 5, K = 5, M = 10, which results in an overall function evaluation budget of 255, although only the first 250 function evaluations are shown in Figures 1 and 2 .For the random search algorithm, we sampled uniformly with bounds of [-5, 5] across all dimensions, for both experiments. We sampled 10 points, over a total of 25 iterations, for a total of 250 function evaluations in total.For the LFBO algorithm, we broadly followed the specifications provided by code in the original paper [13] . We used a two layer Neural Network, with 32 hidden units, which is optimized with the Adam optimizer with a fixed learning rate of 10 -3 . The acquisition function was optimized with the BOTORCH package [3] in Python.Submitted to Workshop on Bayesian Decision-making and Uncertainty, 38th Conference on Neural Information Processing Systems (BDU at NeurIPS 2024). Do not distribute.By our construction, z and x k-1 are conditionally independent given x k-1 , so p(z|x, x k-1 ) = p(z|x).Recall we assume z and x k-1 are independent and this assumption of independency p(z = 1|x) = p(z = 1|x, x k-1 ) is also required in proof of Corollary 1.Inspired by the training of Noise Conditional Score Network (NCSN)[14], when the budget is considerably large we can set the learning rate η ∝ 1/σ 2 to cancel the unstability induced by variance annealing. 
 
 Appendix:  



Reviews: review: title: title: Peer review for Lightspeed Black-box Bayesian Optimization via Local Score Matching
main_review: main_review: ## Summary

This paper presents a method for improving Bayesian Optimization (BO) by addressing computational complexity and scalability issues associated with high-dimensional black-box function evaluations. The authors propose a fast acquisition function maximization procedure using local score matching to estimate the gradient of the Probability Improvement (PI) acquisition function. The approach allows for gradient-based optimization without the computational burden of fitting typical neural networks. The effectiveness of this approach is supported by theoretical developments and proven by empirical evaluations.

The originality and significance of the paper are present. However, it is recommended the authors to strengthen the result discussion section, which further sheds light on the strengths and limitations of the proposed approach and paves the road for future research.

## Big-picture comments
* **Originality**: the originality of this paper comes from the use of local score matching for gradient estimation in Bayesian Optimization. While the paper specifically focuses on PI acquisition function, the authors pointed out potential future direction to expand into Covariance Matrix Adaptation - Evolutionary Strategy (CMA-ES). 
* **Research set-up**: the paper provides both a theoretical foundation and empirical evidence for local score matching method. The paper should, however, include more detailed discussion on the conditions under which the method performs to ensure reproducibility.
* **Result discussion**: The paper could benefit from a more detailed discussion of the results. Particularly, in figure 1, score matching (SM) shows a more promising results than LFBO, compared to figure 2. How should this result be interpreted and what strengths and limitations does it reveal?

## Questions
1. Where do the authors see the local score matching approach the most effective? 
2. How should the results in figure 1 and figure 2 be interpreted?
3. What challenges do authors foresee in adapting this approach to other acquisition functions beyond the PI?
title: title: Review of "Lightspeed Black-box Bayesian Optimization via Local Score Matching"
main_review: main_review: **Summary:**
This paper proposes a new and efficient approach to optimize the Probability Improvement (PI) acquisition function for Bayesian Optimization. The key contribution is the use of local score matching to estimate the gradient of the PI acquisition function, allowing for gradient ascent without the expensive covariance matrix inversion that typically makes GP-based methods computationally intensive. The method offers computational efficiency with complexity reduced to $O(n)$, an improvement over the standard $O(n^3)$ required for GP-based methods. This paper demonstrates the efficiency of their local score matching method through experiments on synthetic objectives, comparing it against the state-of-the art $O(n)$ Likelihood-Free Bayesian Optimization (LFBO) and random search.

**Pros:**
- Reducing the complexity to $O(n)$ using local score matching is a valuable improvement over traditional GP-based methods, which typically have $O(n^3)$ complexity due to covariance matrix inversion. While there are other state-of-the-art methods such as BORE and LFBO, these approaches often rely on neural networks and suffer from challenges like overfitting in high dimensions or complex numerical optimization. The proposed method addresses these issues by using local score matching, offering a simpler and more scalable alternative.
- Applying local score matching to estimate the gradient of the PI acquisition function is a novel technique that provides a new way to optimize acquisition functions without expensive surrogate model inference.
- Providing some derivations to support the soundness of the gradient estimation via local score matching. Specifically, the authors present a closed-form solution for the score matching objective and prove that the estimator is asymptotically unbiased under certain regularity conditions. 
- Validating the efficiency of their proposed method on two synthetic objective functions with dimensions ranging from 5 to 20, demonstrating its scalability in medium-high dimensions.

**Cons:** \
**Methodology:**
- The statement, "We propose to directly maximize the PI acquisition function by performing gradient ascent," could be misleading. While the approach indeed uses gradient ascent, the gradient computation is not direct but is instead estimated via local score matching. It's better to explicitly mention that "direct" refers to the difference from the standard approach, which involves GP regression to obtain the gradient.
- The paper lacks a clear explanation of why the proposed method has computational complexity $O(n)$. A brief discussion explaining the difference between the $O(n^3)$ complexity of the standard method (due to covariance matrix inversion) and the $O(n)$ complexity of the proposed method would make this more clear.
- While the paper admits that the method is restricted to the PI acquisition function, it would be useful to provide a short explanation of why this limitation exists. Specifically, this could involve discussing why the local score matching technique works well for PI but is not easily extendable to other acquisition functions like EI.

**Experiment:**
- The paper uses immediate regret as a performance metric and briefly explains what it is. However, it would be better to add a few words about why immediate regret is considered over other metrics, such as simple regret or cumulative regret.
- The authors mention that a different acquisition function was used for LFBO in the experiments, but there is no explanation of why this was done and what acquisition function was used. 
- While the paper mentions their proposed method is supposed to be more efficient than the standard PI optimization methods, it would still be helpful to add the standard PI optimization methods as baselines in the experiments.
- The plots in Figure 1 and Figure 2 lack clarity regarding how the results are aggregated and the variability across different runs, e.g., whether they are showing the mean with standard deviation/error bars or the median with other quartiles.
- More thorough experiments can be added. For example, higher dimensions, additional benchmarks, real-world problems.