Omnipresent Yet Overlooked: Heat Kernels in Combinatorial Bayesian Optimization

We are gratified to see the reviewer found the paper to be well-written, and our unifying framework to be "very interesting and important". We thank the reviewer for this thoughtful response, and reply to the mentioned points below.

I do not think no-ordinal variables is a standard in combinatorial BO papers (line 350). In my opinion, many recent papers either treat ordinal variables as categorical or avoid addressing them due to the additional complexity involved, especially when their methods already demonstrate strong performance on categorical inputs. Hence the ordinal variables have received less attention. In fact, many real-world tasks require the ordinal variables, notably in HPO (number of trees or max depth in random forests, number of layers and units in neural networks, etc.). Examples in other fields include [1-3], which have all been used in BO works [4-5].

We believe a slight misunderstanding took place, and thank the reviewer for bringing it to our attention. Since different sources rely on different definitions for ordinal variables, we now explicitly specify ours: "Ordinal data is a categorical, statistical data type where the variables have natural, ordered categories and the distances between the categories are not known." This is the main definition on Wikipedia [6], which cites [7] as being the source for this definition. We have carefully examined all the variables in the suggested benchmarks [1-3], and conclude, to the best of our knowledge, that none of them are ordinal variables (according to our definition). As a general example, we consider the number of layers in neural networks. Here, although the variable is discrete, there is a clear distance between the different categories (namely the Euclidean distance). Therefore, according to our definition, this is not an ordinal variable. For lack of a better word, we will refer to these discrete variables (containing an inherent notion of distance) as "discrete-quantitative variables".

In light of the above clarification, we believe our claim holds true: "our experiments do not contain ordinal variables, but we note that this is standard in recent combinatorial BO papers" (line 350). Following this insightful discussion with the reviewer, we have updated the camera-ready version of the paper: we have included a more explicit categorization of the different possible discrete variables, as well as an expanded discussion of how existing works fit into this categorization.

As pointed out by the reviewer, combinatorial BO papers do often contain discrete-quantitative variables. For such variables, we believe the associated kernels are rather straightforward: use any (continuous) isotropic kernel, along with an appropriate distance function (usually Euclidean). As such, these kernels can simply be seen as continuous kernels, which is why we did not explicitly address them in our unifying framework. Additionally, although discrete-quantitative variables can borrow their kernels from continuous BO, we note that this is not true for their acquisition-function optimizers, which need to be discrete (e.g. genetic algorithm). We are grateful to the reviewer for highlighting this additional type of discrete variable, and have included an additional paragraph in the camera-ready version that discusses this relevant topic.

Lastly, for ordinal variables (according to our definition of the term), we would like to mention our answer provided to Reviewer b9FW. In this answer, we explain how to model ordinal variables using path graphs $\mathcal{G}_i$, therefore connecting them to our existing heat-kernel framework.

[6] https://en.wikipedia.org/wiki/Ordinal_data

[7] Agresti, Alan. Categorical Data Analysis. John Wiley & Sons (2002).

Given the simplicity of heat kernel pipeline, can the authors implement the discussed new ideas in Hamming-based approach in Sec. 3.4.1(via rational quadratic kernel), and similarly for Graph-based kernels in Sec. 3.4.2 (via graph Matern kernel and sum-of-inverse polynomial kernel)?

We have implemented two of the above three kernels: the Hamming-based rational-quadratic kernel and the graph-based Matérn kernel. To the best of our understanding, the sum-of-inverse polynomial kernel is fundamentally an ARD kernel, and turning this into a non-ARD kernel would not be sensible. Since our experimental evaluation considers only non-ARD kernels, we believe it would be confusing for the reader to add an ARD kernel, which is why we have omitted it.

Unfortunately, the NeurIPS policy this year does not allow us to show or link images in any way whatsoever. Fortunately, in this case, it should not be a problem, as the results are easy to describe: on all seven benchmarks (including the additional two suggested by Reviewer opDa), the Hamming-based rational-quadratic kernel and the graph-based Matérn kernel achieve near-identical performance compared to the five other Hamming kernels evaluated in the paper. As a result, we are now left with one single rainbow-colored line, which includes seven different kernels that all obtain the exact same empirical performance. To the best of our knowledge, we believe we are the first to demonstrate the empirical (and partly theoretical) equivalence between these seven kernels for combinatorial BO, which were previously linked to distinct methods.

Moreover, in Section 3.4 of the paper, we show that Hamming- and graph-based kernels can be seen as broader (generalized) versions of heat kernels, and that these two classes also coincide (for $|\mathcal{X}_1| = \ldots = |\mathcal{X}_n|$). With these new experiments, we now have additional empirical evidence for these theoretical results. We believe these experiments meaningfully improve the paper, and sincerely thank the reviewer for suggesting them.

Is there any way to make heat kernel enjoy the same benefit of BODi and BOSS kernels such that the performance is also great in problems with special structure? It would be much better if a kernel can work well in both scenarios.

The short answer is "yes", and we see two ways of doing this: incorporating group-invariances, or including additional features.

When it is known that the problem at hand has a certain type of structure, this structure can often be modeled by incorporating invariance to a certain group $G$ directly into the kernel. As explained in Section 3.5 and Appendix C.1, our heat-kernel framework allows for the incorporation of any group-invariance, and this through a straightforward procedure. In a way, this is also what the SSK (kernel from BOSS) does, since it models certain string-invariances, which seems to be a powerful inductive bias for certain types of optima. In Figures 5 and 7, we empirically demonstrate the advantages of incorporating such group-invariances, resulting in state-of-the-art performance. However, in general, modeling invariances can quickly become expensive, as can be seen by BOSS' high computational cost in Figure 6 (and the fact that it needs highly-vectorized code and high-memory GPUs to get to this speed).

A more lightweight alternative, though perhaps less principled, is to include additional "strategic" features. For example, in the Bounce paper [8], the authors argue that BODi achieves its performance mostly due to its way of doing "uniform" sampling of anchor points, which turns out to be heavily skewed towards points with a certain type of structure. As a consequence, BODi's HED embeddings or features $\phi_{\mathbf{A}}(\mathbf{x})$ contain a clear distance metric to anchor points with a certain structure, which can become a useful signal to navigate towards certain optima. The same idea can easily be integrated within our heat-kernel framework: sample anchor points according to the desired type of structure, and add the relevant features (i.e. the Hamming distance to these anchor points) to the existing heat-kernel. We believe this general idea might be an interesting direction for further research, and have added it to the camera-ready version of the paper. Once again, we appreciate the reviewer's efforts to improve the paper and make it more complete.

[8] Papenmeier, Leonard, Luigi Nardi, and Matthias Poloczek. Bounce: Reliable High-Dimensional Bayesian Optimization for Combinatorial and Mixed Spaces. NeurIPS (2023).

We are pleased to see the reviewer found the paper to be "clearly written" and our framework to be "theoretically sound". We are grateful for the reviewer's input and will be addressing their comments.

It could be interesting to consider more real-world black-boxes (such as the ones from the MCBO library); Please evaluate on those:

In the MCBO library, there a four categorical benchmarks left on which we did not evaluate: two biological benchmarks (RNA Inverse Folding and Antibody Design) and two logic-synthesis optimization benchmarks (AIG Optimization and MIG Optimization). Since Antibody Design and MIG Optimization contain much slower simulators (MIG Optimization can be up to $200.000\times$ slower), we focus our efforts in this limited time-window on RNA Inverse Folding and AIG Optimization. As a result, the camera-ready paper now also contains one biological and one logic-synthesis optimization task. We believe this makes the experimental section of the paper more complete, and are grateful to the reviewer for this suggestion.

Unfortunately, the NeurIPS policy this year does not allow us to show or link images in any way whatsoever. Consequently, we use text to describe the main trends observed on the new plots:

1. On both RNA Inverse Folding and AIG Optimization, all Hamming kernels obtain nearly indistinguishable performance (as is the case for the other datasets presented in the paper).
2. On both datasets, SSK (BOSS) and HED (BODi) do not seem particularly sensitive to relocation of the optima (a similar behavior was observed on LABS in Figures 1 and 2 of the paper).
3. On both datasets, our fast and simple heat-kernel pipeline achieves identical performance to Bounce, which is currently the state-of-the-art method. Interestingly, although both methods (Bounce and ours) are the top performers on AIG Optimization, this is not the case on RNA Inverse Folding. On this last dataset, there is a new trend: additive kernels, namely CoCaBO and RDUCB, are able to outperform Hamming kernels for the first time. This is not necessarily surprising: in related biological tasks, there exists evidence that first- and second-order interactions lead to strong predictive performance [1-2]. The MCBO paper also contains potential evidence for this low-order hypothesis: a genetic algorithm, which only mutates a few variables at a time and does not inherently model high-order interactions, performs remarkably well on RNA Inverse Folding, significantly outperforming all BO methods. Lastly, we note that if such additive structure is known a priori, one can easily incorporate it into the heat-kernel, leading to methods such as (but not limited to) RDUCB and CoCaBO (see Appendix C.2).

[1] Domingo, Júlia, Guillaume Diss, and Ben Lehner. Pairwise and higher-order genetic interactions during the evolution of a tRNA. Nature (2018).

[2] Faure, Andre J., et al. The genetic architecture of protein stability. Nature (2024).

line 229: it should be "by choosing k", not "mathbb{k}"

We acknowledge that our notation in the mentioned paragraph is confusing, and are thankful to the reviewer for bringing this to our attention. However, we do not believe that replacing "mathbb{$k$}" by "$k$" is correct. In the updated camera-ready version, we have now written: "by choosing mathbb{$k$}$(d)$ to be any (continuous) isotropic kernel (with $d = \lVert\mathbf{x} - \mathbf{x}'\rVert^2_2$), we always obtain a valid Hamming kernel $k$." We are happy to address any further questions the reviewer might have.

Have a line style distinct for heat-based and non-heat-based to ease the analysis

We thank the reviewer for this great suggestion. We have made use of it for the camera-ready version of the paper, which now contains improved plots.

Is there a way to see the kernel from BODi as a heat kernel variation in some sense?

We appreciate the reviewer's excellent question, and believe the answer is "no". In fact, we have constructed a short sketch for why BODi is not part of the Hamming-kernel class, and, as a consequence, cannot be seen as a (generalized) heat kernel. This proof has now been included in the appendix, and is referenced in the main part of the paper. We believe this is a meaningful addition to the paper, and thank the reviewer for highlighting this point.

As outlined in Definition 6: a kernel is a Hamming kernel (HK) if and only if it depends on $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$ through the square root of the Hamming distance $h(\mathbf{x}, \mathbf{x}')$. That is, $k_{`HK`}(\mathbf{x}, \mathbf{x}') := f\left(\sqrt{h(\mathbf{x}, \mathbf{x}')}\right),$ for some $f : \mathbb{Z}_0^+ \rightarrow \mathbb{R}$.

In contrast, the BODi kernel depends on $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$ through the term $k_{`BODi`}(\mathbf{x}, \mathbf{x}') := g\left(\lVert\phi_{\mathbf{A}}(\mathbf{x})-\phi_{\mathbf{A}}(\mathbf{x}')\rVert^2_2\right),$ with $g(d): =\left( 1 + \frac{\sqrt{5} d}{\ell} + \frac{5 d^2}{3 \ell^2} \right) \exp \left( -\frac{\sqrt{5} d}{\ell} \right).$ Here, $g(\cdot)$ is the Matérn-5/2 function, where we note that $g(\cdot)$ constitutes a bijective function. Additionally, $\phi_{\mathbf{A}}(\cdot)$ denotes the mapping towards the HED embedding (see section 4 of [3]), namely: $\left[ \phi_\mathbf{A}(\mathbf{x}) \right]_i := h(\mathbf{a}_i, \mathbf{x}),$ with $\mathbf{a}_i \in \mathbf{A}$ one of the $M$ anchor points in dictionary $\mathbf{A}$.

For BODi to be a Hamming kernel, there thus needs to exist a function, say $u(\cdot)$, such that $u\left(\sqrt{h(\mathbf{x}, \mathbf{x}')}\right) = \lVert\phi_{\mathbf{A}}(\mathbf{x})-\phi_{\mathbf{A}}(\mathbf{x}')\rVert^2_2.$ If this function $u(\cdot)$ exists, we can define $f:= g \circ u$, obtaining

$$k_{`BODi`}(\mathbf{x}, \mathbf{x}') = g\left(\lVert\phi_{\mathbf{A}}(\mathbf{x})-\phi_{\mathbf{A}}(\mathbf{x}')\rVert^2_2\right) = g\left(u\left(\sqrt{h(\mathbf{x}, \mathbf{x}')}\right)\right) = f\left(\sqrt{h(\mathbf{x}, \mathbf{x}')}\right) = k_{`HK`}(\mathbf{x}, \mathbf{x}').$$

However, as becomes clear by looking at the mapping $\phi_{\mathbf{A}}(\cdot)$, there does not exist a function $u(\cdot)$ satisfying the above claim. We provide a small counter-example with $M=1$ below.

In scenario A, we have: $\mathbf{x} = [1,1,1], \text{ } \mathbf{x'} = [1,1,2] \text{ and } \mathbf{a}\_1 = [1,1,1],$ which gives us $\sqrt{h(\mathbf{x}, \mathbf{x}')} = 1 \text{ and } \lVert\phi_{\mathbf{A}}(\mathbf{x})-\phi_{\mathbf{A}}(\mathbf{x}')\rVert^2_2 = 1.$

In scenario B, we have: $\mathbf{x} = [1,1,2], \text{ } \mathbf{x'} = [1,1,3] \text{ and } \mathbf{a}\_1 = [1,1,1],$ which gives us $\sqrt{h(\mathbf{x}, \mathbf{x}')} = 1 \text{ and } \lVert\phi_{\mathbf{A}}(\mathbf{x})-\phi_{\mathbf{A}}(\mathbf{x}')\rVert^2_2 = 0.$

For $u(\cdot)$ to work in both scenarios, we would need to have: $u(1) = 1 \text{ and } u(1) = 0,$ which is clearly not possible. As a result, we conclude that BODi is not a Hamming kernel, and therefore also not a heat kernel (or a direct generalization thereof).

In fact, using the above proof-structure, we can show an even stronger result: BODi is not part of an (even broader) general class of kernels, which includes not only Hamming- and graph-based kernels, but also additive kernels such as CoCaBO and RDUCB (see Appendix C.2). With this, we feel confident in our answer that BODi cannot be seen as a heat-kernel variation in some sense.

[3] Deshwal, Aryan, et al. Bayesian Optimization over High-Dimensional Combinatorial Spaces via Dictionary-Based Embeddings. AISTATS (2023).

We are excited and encouraged to see the reviewer's strong positive response, and are happy to reply to their questions below.

I do not think there are any significant weaknesses, however, I believe the presentation of the paper could be slightly improved. A lot of important results and insights appear only in the appendix. Authors should think carefully about how to best use the additional page in camera-ready version and move the most important parts of the appendix there. It would also be nice to report the number of seeds used for each run in the main text.

We would like to thank the reviewer for highlighting this aspect. During this rebuttal period, we have taken great care to migrate the most important parts of the appendix to the additional page of the camera-ready version. We have also made sure to include the number of seeds of each run in the main paper.

(typo): I believe equation 4 should have a minus inside the exponent

We politely disagree with the reviewer, and provide a short derivation below to explain our intuition. Here, we start with an RBF kernel applied to one-hot encoded data, namely $\mathbf{z}, \mathbf{z}' \in \\{0,1\\}$: $\exp\left(-\frac{\lVert\mathbf{z} - \mathbf{z}'\rVert^2_2}{2\ell}\right).$ As expected, this kernel has a minus inside the exponent. Afterwards, we use the equivalence from Equation 1 in the paper: $h(\mathbf{x}, \mathbf{x}') = \frac{\lVert\mathbf{z} - \mathbf{z}'\rVert_2^2}{2},$ where $h(\mathbf{x}, \mathbf{x}')$ is the Hamming distance between categorical variables $\mathbf{x}, \mathbf{x}' \in \mathcal{X}$. Using this equivalence, and explicitly writing out the Hamming distance in terms of Kronecker-delta functions $\delta(\cdot,\cdot)$, we get:

$$\exp\left(-\frac{\lVert\mathbf{z} - \mathbf{z}'\rVert^2_2}{2\ell}\right) = \exp\left(-\frac{h\left(\mathbf{x}, \mathbf{x}'\right)}{\ell}\right) = \exp\left(-\frac{1}{\ell} \left[n - \sum_{i=1}^n \delta(x_i, x_i')\right] \right) \propto \exp\left(\sum_{i=1}^n \frac{\delta(x_i, x_i')}{\ell}\right),$$

where the minus sign has now disappeared. If we replace $\ell$ by $\ell_i$ to allow for ARD, and re-parameterize with $\ell_i := n/\gamma_i$, we obtain $\exp\left(\frac{1}{n}\sum_{i=1}^n \gamma_i \delta(x_i, x_i')\right),$ which is exactly the kernel proposed in Equation 1 of [1] (and presented in Equation 4 of our paper). We hope that this derivation was able to clarify the reason for the missing minus sign, and are happy to address any remaining questions the reviewer might have.

[1] Wan, Xingchen, et al. Think Global and Act Local: Bayesian Optimisation over High-Dimensional Categorical and Mixed Search Spaces. ICML (2021).

We are delighted to see the reviewer found the paper to be well-written, and the problem to be relevant. We appreciate the reviewer's insights, and are happy to answer their questions.

Why do we need Figures 2 and 4? Aren't they showing mostly the same thing?

Yes, that is correct: Figures 2 and 4 depict the same baselines, only on different datasets. The same is true for Figures 1 and 3. Due to space constraints, we decided to move some datasets to the appendix. In the camera-ready version, we have now more clearly described the relationship between these two pairs of figures (2/4 and 1/3). We are thankful to the reviewer for pointing out this possible point of confusion.

Do you have concrete ideas on how one could extend the analysis to ordinal problems?

On line 139 of our submitted paper, we explain how to transform the discrete domain $\mathcal{X}$ into a graph $\mathcal{G}$. Specifically, we define the graphs $\mathcal{G}_i$ as complete graphs, and argue that this is a natural choice when the underlying variables $x_i \in \mathcal{X}_i$ are categorical variables. However, for other types of discrete variables, the graphs $\mathcal{G}_i$ can (and probably should) have different structures (as hinted at on line 143). For example, for an ordinal variable $x_i \in \mathcal{X}_i$, it is natural to define $\mathcal{G}_i$ as a path graph, where nodes are not equidistant anymore (whereas this is the case for complete graphs). In fact, this path-graph structure for ordinal variables is also suggested in the COMBO paper [1].

Using the above idea, there is a clear way to connect ordinal variables to our unifying framework. First, define the graphs $\mathcal{G}_i$ according to the nature of their underlying variables $x_i \in \mathcal{X}_i$ (complete graph for categorical, path graph for ordinal, etc.). Second, create the graph $\mathcal{G}$ as the graph Cartesian product of all $\mathcal{G}_i$ (similar to what is done in Definition 2 of our paper). Lastly, apply the diffusion kernel from Equation 2 or 5 on this product-graph $\mathcal{G}$. This connection has now been made more explicit in the camera-ready-version, where we have greatly expanded our discussion regarding ordinal variables. We believe this is a meaningful addition, and are grateful to the reviewer for highlighting this aspect.

Unfortunately, the above connection is also where the analysis stops. Since path graphs do not contain the same type of symmetries as complete graphs, the eigenvalues of the corresponding Laplacian do not simplify as much as for categorical variables. Therefore, we do not obtain a simple closed-form solution. In other words, Equation 2 or 5 cannot be simplified to something like Equation 3 or 4.

Lastly, we would like to refer to our answer given to Reviewer oZnX, where we argue that most combinatorial BO papers do not include truly ordinal variables. As a consequence, we cannot unify many well-known kernels into one unifying framework (as we did for categorical variables), since there simply do not exist many kernels designed for ordinal variables. We have now stressed this argument in the updated version of the paper.

[1] Oh, Changyong, et al. Combinatorial Bayesian Optimization using the Graph Cartesian Product. NeurIPS (2019).