Black-Box Gradient Matching for Reliable Offline Black-Box Optimization

1. Assumption that the oracle has gradient.

We agree with the reviewer that not all oracle functions possess a gradient.

However, we want to point out that the implicit premise of the entire offline optimization field is that the developed techniques are meant for oracle functions which are continuous and differentiable. This premise is set the moment we set out to optimize the oracle using a differentiable surrogate such as the neural network.

This is the common theme in most (if not all) prior work.

This is also not an unreasonable assumption because attempting to optimize an oracle function with discontinuities is otherwise an ill-posed problem. Indeed, without the assumed differentiability, the oracle can be discontinued at any unseen data points at which its output behavior can be arbitrary, giving us no solid ground to relate between the training and (unseen) test output.

With this, we hope the reviewer would agree with us that we had in fact not used any assumptions prior work had not used, and consider re-evaluating the rating of our work.

2. Discussing and providing empirical comparisons with suggested baselines.

We would like to point out that we did compare with CMA-ES. We will include in our revision the following positioning with BO-qEI, BDI, NEMO, and IOM.

(1) BDI uses forward and backward mappings to distill knowledge from the offline dataset to the design. BDI employs infinite width (neural tangent kernels) neural nets as surrogate models which might be limited in representation power for some applications.

(2) IOM uses domain adaptation methods to enforce representation invariance between the training dataset and the distribution of optimized designs. Figuring out the right notion of distribution over optimized designs is an open challenge and the one used by IOM (points found during gradient ascent with a given surrogate) might not be the most effective one.

(3) BO-qEI is standard Bayesian optimization (BO) run with a surrogate model as the objective function.

(4) NEMO is an uncertainty quantification approach that uses normalized maximum-likelihood (NML) estimator which requires quantizing the output scores y. This quantization might be sub-optimal.

Furthermore, we have also compared our method with BO-qEI & BDI and reported the results below.

100-th percentile

50-th percentile

Note: OOM means out of memory.

The results show that MATCH-OPT performs the best in 6 out of 12 cases (across both the 100-th and 50-th percentile settings) while BO-qEI only performs best in 1 out of 12 cases. BDI performs best in 5 out of 12 cases, runs out of memory in 2 out of 12 cases. Overall, MATCH-OPT appears to perform more stable than BDI and is marginally better than BDI. It is also more memory-efficient than BDI as it does run successfully in all cases, while BDI runs out of memory in 2 cases. MATCH-OPT also outperforms BO-qEI significantly in 11 out of 12 cases. These new results will be incorporated into our main text.

In addition, we would also be happy to run comparative analysis with NEMO and IOM but the code for those baselines has not been released. We hope the reviewer would sympathize with us on this point since offline optimization is a fast-growing field and providing an exhaustive comparison with all prior work, including those with no release code is too difficult. Nonetheless, we have tried our best to compare with most prior work that provides released code.

We hope the reviewer would take into account such practical difficulties and re-consider the rating of our work. We are looking forward to hearing back from you & will be happy to continue the discussion if you still have other questions for us.

3. Restructuring the experiment to make the introduction more concise. We agree with the reviewer and will adjust the content accordingly in our revision.

Thank you again for the detailed review & suggestions.

Thank you for the quick follow-up. From your response, you said that

COMs does not use the "oracle differentiable assumption" in its method design and only uses it for analysis. This also holds for BONET

This is true for us as well. If you take a look at Eq. (11), it does not involve the oracle's gradient anywhere. It only involves the gradient of the surrogate and the output difference between two sampled data points.

Furthermore, you have also agreed that:

Even if the objective function does not have gradient (for example, has some discontinuities), we can still use NN to approximate it

This means it is fine to use the NN's gradient in place of the objective function's gradient (even if it does not exist). This is no different from our approach. So, once again, we do not use any assumptions that prior methods did not use.

We hope the reviewer could reconsider the rating given the above point

Thank you for increasing the score.

We understand that the reviewer’s suggestion is although (implementation-wise) we can use the NN’s gradient in place of the objective function’s gradient, it is better to write in the paper that

We have an ideal differentiable proxy which is the closest one to the ground-truth and this ideal proxy is what is called the oracle.

We have put this sentence in a new revision of our main text (see Section 2). This will solve the problem. The additional experiments in our original response will also be included. With this, we hope the reviewer would still reconsider the assessment that this paper needs a lot of refinement to publish.

--

In addition, we respectfully do not agree with the reviewer that COMS does not use the gradient assumption in its method design. It is written at the beginning of its Section 3 that

COMs learn estimates of the true function that do not overestimate the value of the ground truth objective on out-of-distribution inputs in the vicinity of the training dataset. As a result, COMs prevent erroneous overestimation that would drive the optimizer (Equation 2) to produce out-of-distribution inputs with low values under the ground-truth objective function

The last part is true only if its main theorem holds, which is under the assumption that the true function is L-Lipschitz, which implies the gradient assumption. In fact, the regularizer term in the objective function -- see Eq.(4) -- is equivalently transformed to the more expressive form in Eq. (7) which is exactly what is bounded in the main theorem. Thus, similar to our paper, the intriguing theoretical analysis of COMS informs its algorithm design in this case.

What COMS assumed is fundamentally not different from the approximation step we used in Eq. (10). More importantly, we can remove the first equality in Eq. (10) so that there is no explicit writing of $\nabla g(\mathbf{x}) \simeq \nabla g_\phi(\mathbf{x})$ in the algorithm design and the rest of the technical narrative is still correct. But, we think the clarification statement above is already enough to make it clear.

Thanks for the authors' detailed response. Part of my concerns have been addressed. I am happy to see the additional insight provided by the authors. However, some of my concerns remain unsolved.

$1. **Optimizing Eq. (11) might not guarantee gradient match**$

It seems that the authors assume $\mathbf{F}_{\phi}(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) = 0$ for all $(\mathbf{x}, \mathbf{x}’)$.

However, we may only guarantee $\mathbf{F}_{\phi}(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) = 0$ for finite $(\mathbf{x}, \mathbf{x}’)$ in the training set, even if we assume a strong enough optimizer.

As a result, we only know $g(\mathbf{x}) - g_\phi(\mathbf{x}) = c$ for finite samples $\mathbf{x}$ in the training set instead of $\forall \mathbf{x} \in \mathbb{R}^d$.

Thus, we may not conclude $\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x}) = 0$ for all $\mathbf{x}$, especially for the high-dimensional cases that the training data is too sparse to cover the entire domain.

1. Optimizing Eq. (11) might not guarantee gradient match. Thank you for the question. We will prove below that optimizing Eq. (11) does necessarily guarantee gradient match to address this concern.

First, we understand the reviewer's key concern here is that there might exist a parameter $\phi$ such that the integrated vector

$\mathbf{F}_{\phi}(\mathbf{x},\mathbf{x}') \ = \ \int_0^1 \nabla g(t\mathbf{x} + (1-t)\mathbf{x}')\mathrm{d}t$

$- \int_0^1 \nabla g_\phi(t\mathbf{x} + (1-t)\mathbf{x}')\mathrm{d}t$

is non-zero but instead orthogonal to $\Delta\mathbf{x} = \mathbf{x} - \mathbf{x}'$. We can show that such a parameter does not exist in what follows.

First, following your argument, for such a parameter phi, we have

$\mathbf{F}_{\phi}(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) = 0$ for all $(\mathbf{x}, \mathbf{x}’)$

Next, by the line integration theorem, we have:

(1) $g(\mathbf{x}) - g(\mathbf{x}’) = (\mathbf{x} - \mathbf{x}’)^\top \int_0^1 \nabla g(t\mathbf{x} + (1-t)\mathbf{x}’)\mathrm{d}t$

(2) $g_\phi(\mathbf{x}) - g_\phi(\mathbf{x}’) = (\mathbf{x} - \mathbf{x}’)^\top \int_0^1 \nabla g_\phi(t\mathbf{x} + (1-t)\mathbf{x}’)\mathrm{d}t$

for all $(\mathbf{x}, \mathbf{x}’)$. As such, combining (1) and (2) leads to

$\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) =$

$g(\mathbf{x}) - g(\mathbf{x}’) - g_\phi(\mathbf{x}) + g_\phi(\mathbf{x}’)$

Or equivalently, $g(\mathbf{x}) - g(\mathbf{x}’) - g_\phi(\mathbf{x}) + g_\phi(\mathbf{x}’) = 0$ since $\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) = 0$ for all $(\mathbf{x}, \mathbf{x}’)$.

We can thus rewrite the above as $g(\mathbf{x}) - g_\phi(\mathbf{x}) = g(\mathbf{x}’) - g_\phi(\mathbf{x}’)$ for all $(\mathbf{x}, \mathbf{x}’)$, which means there must exist a constant $c$ such that $g(\mathbf{x}) - g_\phi(\mathbf{x}) = c$ for all $\mathbf{x}$.

Taking gradient with respect to $\mathbf{x}$ on both sides leads to

$\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x}) = 0$ or equivalently, $\nabla g(\mathbf{x}) = \nabla g_\phi(\mathbf{x})$

This means for any parameter $\phi$ such that

$\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’) = 0$

for all $(\mathbf{x}, \mathbf{x}')$, the gradients of $g(\mathbf{x})$ and $g_\phi(\mathbf{x})$ are guaranteed to match.

Thus, the case that the reviewer is concerned about does not exist. Please let us know if this has addressed your concern.

--

Furthermore, please also refer to Appendix C where we show an ablation study between the training objective with and without the regression term. This addresses the concern of the reviewer that the performance might come solely from the regression term. This is not true according to our experiment in Appendix C: including the regression term improves the performance as expected but even without it, the performance is still very competitive, asserting the impact of the gradient match term, which matches with the implication of Theorem 1.

2. The bound in Theorem 1 seems to grow exponentially in $m$ and hence, might be loose. Thank you for another critical question. We understand the concern here is that the bound seems to grow exponentially in the no. of iterations $m$, and might be loose as a result. To address this concern, we want to assure the reviewer that the bound will NOT grow exponentially in the no. of iterations with a standard choice of the learning rate $\lambda$. Here is how.

Let us choose $\lambda$ such that it is no more than $1/m$, which is the case in all our experiments where $m = 200$ and $\lambda = 10^-4$, as stated in the appendix. With this choice, we have:

$(1 + \lambda \cdot \mu)^{(m - 1)} \leq (1 + \mu / m)^{(m-1)} < (1 + \mu/m)^m$

which will approach $e^mu$ in the limit of $m$. Here, we use the known fact that $\mathrm{lim}_{m\rightarrow\infty}(1 + \mu/m)^m = e^\mu$ with $\mu > 0$.

As such, when $m$ is sufficiently large the bound in Theorem 1 is upper-bounded with

$m \cdot\lambda \cdot \ell \cdot (1 + \lambda \cdot \mu)^{(m - 1)} \cdot \text{gradient-gap} \simeq \ell \cdot e^\mu \cdot\text{gradient-gap}$

which asserts that the worst-case performance gap of our offline optimizer is approaching (in the limit of $m$) $\ell \cdot e^\mu \cdot \text{gradien-gap} = \mathbf{O}(\text{gradient-gap})$ which is not dependent on the number of gradient steps.

Thus, the bound is not growing exponential in $m$ and hence, it is not loose and trivial. Please let us know if this has addressed the reviewer’s concern.

We will address your remaining concerns in the next comment

Thank you for the prompt response. Previously, we had proved that in the limit of data if we can find the true optimal solution that zeroes out the generalized loss, the gradient is matched.

Now, suppose that $\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’)$ is not zero everywhere (i.e., the generalized loss is not zero), we will show that the gradient gap is still bounded by a decreasing quantity.

To see this, note that $\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’)$

$= (g(\mathbf{x}) - g_\phi(\mathbf{x})) - (g(\mathbf{x}’) - g_\phi(\mathbf{x}’))$

which was established in the 5th equation of our previous message.

Let $h(\mathbf{x}) \triangleq g(\mathbf{x}) - g_\phi(\mathbf{x})$. It follows that $|\mathbf{F}_\phi(\mathbf{x}, \mathbf{x}’)^\top(\mathbf{x} - \mathbf{x}’)| = |h(\mathbf{x}) - h(\mathbf{x}')|$ which means our loss function is working towards minimizing $|h(\mathbf{x}) - h(\mathbf{x}')|^2$ over $(\mathbf{x},\mathbf{x}')$. Intuitively, this will make $h(\mathbf{x})$ smoother as the output distance between different inputs are being reduced.

As a result, this process will reduce the Lipschitz constant $\epsilon$ of $h(.)$. Here, we are alluding to the hypothesis that if a model behaves smoothly on the training data, it will also be smooth on the test data. We believe this is reasonable in the broader context as well as this context because most work in offline optimization present some forms of smoothing out the behavior of the model on test data and it is often achieved via conditioning the model on the training data.

Now, we will show that the gradient gap $||\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x})||$ is bounded by $\epsilon$:

First, by definition, $|h(\mathbf{x}) - h(\mathbf{x}')| \leq \epsilon \cdot ||\mathbf{x} - \mathbf{x}'||$ which implies

$|(g(\mathbf{x}) - g(\mathbf{x}')) - (g_\phi(\mathbf{x}) - g_\phi(\mathbf{x}'))| \leq \epsilon ||\mathbf{x} - \mathbf{x}'||$ -- due to the definiton of $h(.)$

Dividing both sides by $||\mathbf{x} - \mathbf{x}'||$ leads to:

$\left|\frac{(g(\mathbf{x}) - g(\mathbf{x}'))}{||\mathbf{x} - \mathbf{x}'||} - \frac{(g_\phi(\mathbf{x}) - g_\phi(\mathbf{x}'))}{||\mathbf{x} - \mathbf{x}'||}\right| \leq \epsilon$

Since this holds for any $(\mathbf{x}, \mathbf{x'})$, we can choose $\mathbf{x'} = \mathbf{x} + t\cdot \mathbf{e_j}$, where $t$ is some scalar and $\mathbf{e_j}$ denotes the $d$-dimensional one-hot vector with the hot component at the $j$-th coordinate. Here, $d$ denotes the input dimension. This gives:

$-\epsilon \leq \frac{(g(\mathbf{x} + t \cdot \mathbf{e_j} ) - g(\mathbf{x}))}{t||\mathbf{e_j}||} - \frac{(g_\phi(\mathbf{x} + t \cdot \mathbf{e_j}) - g_\phi(\mathbf{x}))}{t||\mathbf{e_j}||} \leq \epsilon$

Now, taking $\mathrm{lim}_{t\rightarrow 0}$ on all terms in the inequality and using the definition of directional gradient, we will arrive at:

$-\epsilon \leq (\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x}))^\top \mathbf{e_j} \leq \epsilon$ , or $\left| (\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x}))^\top \mathbf{e_j} \right| \leq \epsilon$

For each value of $j = 1, 2, \ldots, d$ , we will obtain one such inequality. Summing all these inequalities gives:

$||\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x})||_1 \leq \epsilon \cdot d = \mathbf{O}(\epsilon)$

Finally, since norm-2 is upper bounded by norm-1, we also have $||\nabla g(\mathbf{x}) - \nabla g_\phi(\mathbf{x})||_2 \leq \epsilon$.

Thus, the gradient gap will be bounded by the Lipschitz constant of the function gap $h(\mathbf{x}) = g(\mathbf{x}) - g_\phi(\mathbf{x})$ and this constant decreases as we optimize the loss function. In addition, if $\mathbf{x}$ is within the training data, the gradient gap zeros out as established previously.

Please let us know if this has addressed your concern. Thank you again for keeping this discussion alive.

Thanks for the authors' detailed response. However, my concern has not been addressed.

(1). Finite Sample Case

There exists $g_\phi$ that achieves zero loss on the finite training dataset $\mathfrak{D}$, but the gradient gap is large even on the training samples.

Consider a finite dataset $\mathfrak{D}:= \\{ +1,-1 \\}^d \subset \mathbb{R}^d$,

The loss in Eq.(11) on $\mathfrak{D}$ is given below: (Note that $\mathfrak{D}$ denotes the training dataset under Eq.(11) in the paper )

$

\mathfrak{L}(\phi) = \sum _{\boldsymbol{x}, \boldsymbol{x}' \in \{ +1,-1 \}^d } ( g(\boldsymbol{x})- g _\phi (\boldsymbol{x}) - (g(\boldsymbol{x}')- g _\phi (\boldsymbol{x}') ) )^2 = \sum _{\boldsymbol{x}, \boldsymbol{x}' \in \{ +1,-1 \}^d } ( h(\boldsymbol{x}) - h(\boldsymbol{x}') )^2

$

A solution that achieves zero loss $\mathfrak{L}(\phi) =0$ on the training set $\\{ +1,-1 \\}^d$ is $h(\boldsymbol{x})= g(\boldsymbol{x}) - g _\phi (\boldsymbol{x}) = c$ for $\forall \boldsymbol{x} \in \\{ +1,-1 \\}^d$, where $c$ denotes a constant.

However, there exist functions $h(\boldsymbol{x})$ that achieves zero loss $\mathfrak{L}(\phi) =0$ on the training set $\\{ +1,-1 \\}^d$, but the gradient gap is large even on the training samples.

For example, a simple function $h(\boldsymbol{x})= (\boldsymbol{x}^\top\boldsymbol{x} )^{5}$. For $\forall \boldsymbol{x} \in \\{ +1,-1 \\}^d$, we have $h(\boldsymbol{x})= (\boldsymbol{x}^\top\boldsymbol{x} )^{5} = d^{5}$. Thus, the loss is zero $\mathfrak{L}(\phi) =0$.

But the gradient of $h(\boldsymbol{x})$ on any training sample $\boldsymbol{x}$ is $\nabla h(\boldsymbol{x}) = 10(\boldsymbol{x}^\top\boldsymbol{x} )^4 \boldsymbol{x}$.

The gradient gap on the training sample $\forall \boldsymbol{x} \in \\{ +1,-1 \\}^d$ is

$

\| \nabla g(\boldsymbol{x}) - \nabla g _\phi (\boldsymbol{x}) \|_2^2 = \| \nabla h(\boldsymbol{x}) \|_2^2 = 100 (\boldsymbol{x}^\top\boldsymbol{x} )^{8} \| \boldsymbol{x} \|_2^2= 100 d^{9}

$

For a small $d=100$, we have the gradient gap is

$

\bf{ \| \nabla g(\boldsymbol{x}) - \nabla g _\phi (\boldsymbol{x}) \|_2^2 = 100 ^{10} = 10 ^{20} }

$

In fact, minimizing loss Eq.(11) on the finite training dataset $\mathfrak{D}$ can not guarantee the local variation, especially for high-dimensional cases.

For the above example with a small $d=100$, even with $2^{100}$ training samples, the gradient gap $\\| \nabla g(\boldsymbol{x}) - \nabla g _\phi (\boldsymbol{x}) \\|_2^2$ on the training sample can grows up to $\bf{10 ^{20}}$.

(2) Infinite sample cases

The key assumption of the infinite sample cases is that we can achieve a constant function $h(\boldsymbol{x})=c$ for $\forall \boldsymbol{x} \in \mathbb{R}^d$ to ensure that the gradient $\nabla h(\boldsymbol{x}) = \boldsymbol{0}$. However, this assumption may be strong and impractical.

Thank you for the strong support of our work. We are very happy to know that you find our work insightful!

To answer your question, organizing training data into monotonically increasing trajectories encourages the model to learn the behavior of a gradient-based optimization algorithm, which in fact produces (mostly) monotonically increasing trajectories. This allows the gradient matching algorithm to focus more on strategic input pairs that are more relevant for gradient estimation.

Furthermore, Eq. (13) combines both the gradient match loss and the regression loss. The regression loss helps as a regularizer that amplifies the importance of accurate gradient matching along the sampled trajectories which encodes our bias that the optimal trajectories would be sufficiently close to one of those sampled trajectories.

1. Compare running time with other baselines. Thank you again for the suggestion. We will quote the memory and time complexity of the baselines in the revised version. Some of these are not made explicitly in their corresponding paper but we can look into their formulation to deduce their complexities.

We do, however, want to remark that such complexity comparison is only tangential to our main contribution. Our main focus is on building optimizer with better and more stable performance overall, even at an affordable increase of running time. Furthermore, we want to point out that as some of the baselines (such as BONET) use an entirely different model which has a different number of parameters than ours, the complexity comparison might not be apple-to-apple.

But regardless, we want to report the averaged running time of the baselines below.

All reported running times are in seconds. Our algorithm incurs more time than other baselines but its total running time is still affordable in the offline setting: 4785s = 1.32hr.

Once again, we want to emphasize that (1) it is not our claim contribution to be faster than the baseline; and (2) as we are optimizing expensive-to-evaluate objective functions (e.g., evaluating candidate materials through physical experiments), providing more stable performance (i.e., MNR) in exchange for an affordable increase (around 1hr) of offline computational time is still preferred by the practitioners.

2. Limitations compared to baselines. One potential limitation of our approach in comparison to other baselines is that our gradient match algorithm learns from pairs of data points. Thus, the total number of training pairs it needs to consume grows quadratically in the number of offline data points. For example, an offline dataset with N examples will result in a set of O(N^2) training pairs for our algorithm, which increases the training time quadratically. However, an intuition here is that training pairs are not equally informative and, in our experiments, it suffices to get competitive performance by just focusing on pairs of data along the sampled trajectories with monotonically increasing objective function values. This allows us to keep training cost linearly with respect to N.

On another note, while it is true that none of the existing baselines (including our algorithm) outperform others on all tasks, we believe that at least on these benchmark datasets, our algorithm tends to perform most stably across all tasks, as measured by the mean averaged rank reported in each of our performance tables. This is a single metric that is computed based on the performance of all baselines across all tasks. The end-user can make a judgment based on such metrics. In practice, by looking at how existing baselines perform overall on a set of benchmark tasks that are similar to a target task, one can decide empirically which baseline is most likely to be best for the target task.

3.How does the discretization parameter affect the performance of the algorithm? Our empirical inspections suggest that a discretization $\kappa = 5$ in Eq. 12 is sufficient to get good performance. Further increasing $\kappa$ only results in negligible changes in performance. In addition, increasing kappa will increase the time complexity linearly, as detailed in our Complexity Analysis paragraph following Eq. (13).

We will address your remaining questions in the next comment.