---
title: 'BOCS: Bayesian Optimization for Combinatorial Problems'
url: https://www.emergentmind.com/topics/bayesian-optimization-of-combinatorial-structures-bocs
type: topic
---

# BOCS: Bayesian Optimization for Combinatorial Problems

Bayesian Optimization of Combinatorial Structures (BOCS) is a surrogate-based Bayesian optimization method for expensive-to-evaluate black-box objectives defined over combinatorial domains, with the canonical setting given by binary decision vectors $x \in \{0,1\}^d$. In BOCS, the objective is approximated by a sparse quadratic pseudo-Boolean surrogate, posterior uncertainty over surrogate coefficients is maintained with a sparsity-inducing prior, and the next query is selected by optimizing a sampled quadratic acquisition over the discrete domain. The original formulation couples Bayesian linear modeling with Thompson sampling and a binary quadratic program (BQP) solver, using semidefinite programming (SDP) or simulated annealing (SA) for acquisition-function optimization [1806.08838].

## 1. Problem formulation and conceptual basis

BOCS addresses the setting in which evaluations are expensive, noisy, and budget-limited. The original formulation assumes observations of the form $y = f(x) + \epsilon$, with $\epsilon \sim \mathcal{N}(0,\sigma^2)$, and seeks
$x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$
using as few black-box evaluations as possible [1806.08838]. A closely related minimization form is used in later presentations,
$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$,
with sequential data
$D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$
and $y^{(i)} = f(x^{(i)})$ [2605.10856]. This suggests that maximization and minimization versions differ mainly by sign convention; the central object is the same: a sequential model-based search over a combinatorial space.

The method is designed for discrete-variable black-box optimization, where the number of candidate solutions grows combinatorially and each evaluation is often expensive [2605.10856]. In the original exposition, the domain is binary, but later discussions of BOCS and closely related methods situate it in broader combinatorial settings including sets, sequences, trees, and graphs, typically through binary encodings or one-hot expansions [2008.08177].

The defining modeling assumption is that much of the structure relevant for optimization can be captured by first-order and pairwise interactions. BOCS therefore uses a quadratic surrogate that is linear in its coefficients, trading expressive power for tractability and data efficiency [1806.08838]. This parametric structure is repeatedly identified as the reason BOCS works well when only a small amount of data is available [2605.10856].

## 2. Surrogate model, sparsity prior, and Bayesian inference

The canonical BOCS surrogate for $x \in \{0,1\}^d$ is
$$
\hat f_t(x \mid \alpha)
= \alpha_0 + \sum_{i=1}^d \alpha_i x_i + \sum_{i<j} \alpha_{ij} x_i x_j
= \phi(x)^\top \alpha,
$$
with
$$
\phi(x) = \big[ 1,\ x_1,\ldots,x_d,\ \{x_i x_j\}_{i<j} \big],\quad
\alpha \in \mathbb{R}^{p},\quad
p = 1 + d + \tbinom{d}{2}.
$$
This is the QUBO form over $x \in \{0,1\}^d$ [2605.10856]. Equivalent matrix forms such as $f(x)=x^\top Q x + s^\top x + c$ also appear, and an affine mapping to Ising variables $z \in \{-1,1\}^d$ is standard [2605.10856].

In the original BOCS paper, the data likelihood is Gaussian:
$$
y \mid X,\alpha,\sigma^2 \sim \mathcal{N}(X\alpha,\sigma^2 I_N),
$$
where $X$ is the design matrix built from constant, linear, and pairwise features [1806.08838]. The number of coefficients is $\Theta(d^2)$, so BOCS uses a sparsity-inducing horseshoe prior to control variance when $N \ll p$ [1806.08838]. Later work also describes BOCS as using a horseshoe prior to shrink unneeded coefficients while retaining important ones, with posterior samples obtained via Gibbs sampling [2605.10856]. The vehicle-design study likewise states that BOCS employs a horseshoe prior for sparse parameters and uses sparse Bayesian linear regression to fit the surrogate [2110.00226].

The original paper gives a fully specified hierarchical prior:
- $\alpha_k \mid \beta_k^2,\tau^2,\sigma^2 \sim \mathcal{N}(0,\beta_k^2 \tau^2 \sigma^2)$,
- $\beta_k \sim \mathcal{C}^+(0,1)$,
- $\tau \sim \mathcal{C}^+(0,1)$,
- $P(\sigma^2) \propto \sigma^{-2}$,

together with the Makalic–Schmidt auxiliary-variable reparameterization and conjugate Gibbs conditionals [1806.08838]. The same paper emphasizes that exact sampling of $\alpha$ can be implemented using the sampler of Bhattacharya et al. (2016), giving $O(N^2 p)$ complexity when $N \ll p$ [1806.08838].

This inference scheme is central to the identity of BOCS. The surrogate is quadratic in the binary variables but Bayesian in its coefficients, so uncertainty enters through posterior sampling rather than through a nonparametric posterior over functions. A common misconception is that combinatorial Bayesian optimization is necessarily Gaussian-process-based. In BOCS, the surrogate is instead a sparse quadratic regression model with Bayesian coefficient uncertainty [1806.08838].

## 3. Thompson sampling and acquisition-function optimization

At each iteration, BOCS draws a sample from the posterior over surrogate coefficients and optimizes the resulting sampled surrogate. In the original maximization form, this is
$$
x_t \in \arg\max_{x \in \mathcal{D}} f_{\tilde\alpha}(x) - \lambda \Pi(x),
$$
with typical penalties $\Pi(x)=\|x\|_1$ or $\Pi(x)=\|x\|_2^2$, which coincide on $\{0,1\}^d$ [1806.08838]. In the minimization form used in later work, one samples $\alpha^{(t)}$ from $p(\alpha \mid D_t)$ and computes
$$
x_{t+1}^{\mathrm{BOCS}} \in \arg\min_{x\in\{0,1\}^d} \hat f_t(x \mid \alpha^{(t)}),
$$
which is again a QUBO [2605.10856].

With fixed sampled coefficients, the acquisition reduces to a BQP. In one standard maximization form,
$$
\max_{x \in \{0,1\}^d} x^\top A x + b^\top x,
$$
where $A_{ij}=\tilde\alpha_{ij}$ for $i \neq j$ and $b_i=\tilde\alpha_i-\lambda$ [1806.08838]. Because such BQPs are NP-hard in general, BOCS relies on specialized optimization procedures.

The original BOCS paper introduces an SDP relaxation by mapping $x \in \{0,1\}^d$ to $y \in \{-1,1\}^d$, introducing an auxiliary bias variable, and constructing a lifted matrix $B$ so that the acquisition becomes
$$
\max_{z \in \{-1,1\}^{d+1}} z^\top B z.
$$
This is relaxed to the semidefinite program
$$
\max_{Z \succeq 0} \operatorname{Tr}(B^\top Z)
\quad \text{subject to} \quad
\operatorname{diag}(Z)=1,
$$
followed by randomized hyperplane rounding [1806.08838]. The paper also presents BOCS-SA, which replaces the SDP with simulated annealing over Hamming-1 neighborhoods [1806.08838].

Subsequent descriptions of the original methodology note that semidefinite relaxations with randomized rounding and other binary-optimization heuristics were used in the initial work, while later implementations often employ SA because it is fast in practice for moderate dimensions [2605.10856]. In the vehicle-design study, the acquisition objective is explicitly written as
$$
x^* \in \arg\max_{x \in \{0,1\}^n}
\alpha_0 + \sum_i \alpha_i x_i + \sum_{i<j} \alpha_{ij} x_i x_j,
$$
and optimized by simulated annealing [2110.00226].

The acquisition mechanism is therefore structurally distinctive: posterior sampling provides exploration, while optimization of a sampled quadratic model provides exploitation. This separation is one reason BOCS is often described as data efficient in small-sample discrete optimization [2605.10856].

## 4. Computational properties, strengths, and limitations

BOCS is built around a tractable statistical model, but its tractability is conditional rather than absolute. The quadratic surrogate keeps the number of degrees of freedom dramatically below that of nonparametric surrogates, and the combination of sparsity priors and Thompson sampling is repeatedly identified as the source of strong small-data behavior [2605.10856]. At the same time, the number of coefficients is still $\Theta(d^2)$, and acquisition-function optimization remains a hard combinatorial problem [1806.08838].

The original SDP-based acquisition has substantial computational cost. Later analyses summarize generic SDP solvers for this relaxation as having worst-case time $O(d^6)$, which becomes prohibitive as $d$ grows, especially when categorical variables are represented by one-hot encoding and the effective binary dimension increases [2008.08177]. The same source notes that the approximation gap of the SDP approach deteriorates with dimensionality and cites logarithmic-factor worst-case approximation bounds for general BQPs [2008.08177].

BOCS also has modeling limitations. Its canonical surrogate includes only linear and pairwise interactions, so objectives dominated by higher-order structure induce model misspecification unless the surrogate is explicitly extended [2605.10856]. The 2026 hybrid study uses fully connected HUBO black-box functions,
$$
f_{\mathrm{HUBO}}(x)= \sum_i \alpha_i x_i + \sum_{i<j}\alpha_{ij}x_i x_j + \sum_{i<j<k}\alpha_{ijk}x_i x_j x_k,
$$
precisely to probe this mismatch, while keeping the BOCS surrogate quadratic [2605.10856].

A further practical issue is search stagnation. The 2026 analysis reports that, as the number of observations increases, BOCS tends to repeatedly propose points that have already been evaluated, leading to “learning stagnation” [2605.10856]. There the operational criterion is
$$
\text{declare stagnation at iteration } t
\quad\text{iff}\quad
x_t^{\mathrm{BOCS}} \in \{x^{(i)}\}_{i=1}^{t-1}.
$$
The same work attributes duplicates to posterior concentration, degeneracy among local minima in discrete spaces, and surrogate underfitting of higher-order structure [2605.10856].

These features motivate two standard cautions. First, BOCS is not a universal discrete optimizer; its sample efficiency depends on the adequacy of a sparse quadratic surrogate. Second, the method’s practical success depends as much on acquisition optimization as on Bayesian modeling. Both observations are borne out in the later literature [2008.08177].

## 5. Algorithmic extensions and successor variants

A major line of work improves BOCS by replacing or augmenting its acquisition-function optimizer. “Scalable Combinatorial Bayesian Optimization with Tractable Statistical models” introduces Parametrized Submodular Relaxation (PSR) for the BOCS acquisition step [2008.08177]. Starting from the minimization-form BQP
$$
E(x)=x^\top A x + b^\top x + c,
$$
PSR decomposes $A=A^+ + A^-$, lower-bounds the positive part by an affine function
$$
h_\Lambda(x) =
x^\top (A^+ \circ \Lambda)\mathbf{1}
+ \mathbf{1}^\top (A^+ \circ \Lambda)x
- \mathbf{1}^\top (A^+ \circ \Lambda)\mathbf{1},
$$
with $\Lambda \in [0,1]^{d\times d}$, and then solves
$$
E_{\min}(\Lambda)
=
\min_{x \in \{0,1\}^d}
\left[
h_\Lambda(x) + x^\top A^- x + b^\top x
\right]
$$
exactly by an $s$–$t$ min-cut for fixed $\Lambda$ [2008.08177]. The outer problem maximizes this lower bound over $\Lambda$ using projected proximal gradient updates, and the paper reports that PSR requires dramatically less time than SDP and returns better acquisition values across BQP, contamination, Ising, and LABS benchmarks [2008.08177].

A distinct extension addresses stagnation rather than raw acquisition cost. “Improving search efficiency via adaptive acquisition function selection in discrete black-box optimization” keeps BOCS as the main search framework and activates a Gaussian-process-based module only when a duplicate is proposed [2605.10856]. In that hybrid, the GP uses a Hamming-distance kernel
$$
k(x,x')=\exp\big(-\gamma\, d_H(x,x')\big),
$$
forms a portfolio of Lower Confidence Bound acquisitions
$$
a_m(x)=\mu_{t-1}(x)-\kappa_m \sigma_{t-1}(x), \quad \kappa_m \in \{1,\dots,10\},
$$
and adaptively selects among them with GP-Hedge [2605.10856]. On fully connected QUBO and HUBO benchmarks with $d=50$, the hybrid attains the smallest relative gap in middle and later stages and improves final RelGap by 98.07% and 96.08% over GP-Hedge only and BOCS+Random, respectively, on QUBO, with corresponding HUBO improvements of 86.46% and 79.25% [2605.10856].

Another important comparison is with Bayesian Variational Optimization (BVO), which is presented as a scalable Bayesian optimization method tailored to combinatorial domains [2011.02004]. BVO replaces the quadratic BOCS surrogate with a Bayesian neural network and optimizes acquisition functions by continuous relaxations such as Concrete and Gumbel-Softmax, using pathwise gradients in a relaxed space [2011.02004]. The comparison is explicit: BOCS models up to pairwise interactions unless extended, whereas BVO’s surrogate can represent high-order, nonlinear interactions; BOCS solves a discrete combinatorial program via SDP or SA, whereas BVO performs gradient-based acquisition optimization in a continuous relaxation [2011.02004]. This suggests that later methods have often treated BOCS as the canonical sparse quadratic baseline from which either more scalable solvers or more expressive surrogates depart.

## 6. Empirical behavior and application domains

The original BOCS paper reports strong benchmark performance relative to alternative combinatorial and Bayesian optimization methods [1806.08838]. On binary quadratic programming with $d=10$, BOCS-SDP and BOCS-SA substantially outperform EI, SA, OLS, Random Search, and an MLE-based variant; with correlation length $L_c=10$ and $\lambda=0$, the average simple regret after 100 iterations is reported as approximately $2.54 \pm 0.51$ $(\times 10)$ for EI and $0.07 \pm 0.05$ $(\times 10)$ for both BOCS-SA and BOCS-SDP [1806.08838]. On sparsification of Ising models, BOCS-SDP achieves the best average values and lowest variance across several $\lambda$ values, and on contamination control it consistently finds the best or near-best schedules [1806.08838]. The same paper reports that, on the Ising benchmark, BOCS-SA and BOCS-SDP reach targets 4–7× faster than EI in wall-clock time [1806.08838].

Later work continues to use BOCS as a primary baseline. In the BVO paper’s Ising sparsification benchmark with a $4\times4$ grid and 24 binary decisions, BVO achieves lower objective values than COMBO, TPE, SA, and BOCS-SDP, while in contamination control BVO is close to COMBO and slightly better than BOCS-SDP [2011.02004]. In the pest-control benchmark with 5 categories at each of 21 stages, BVO matches COMBO and is better than RS, TPE, and SA [2011.02004]. These results are relevant to BOCS because they identify settings in which a more expressive surrogate and relaxed optimization outperform or match the classical sparse quadratic approach.

Application work illustrates both the flexibility and the limits of BOCS. In the vehicle-design study, the original objective is the “solution space size” of a vehicle Commonality Description Matrix (CDM), which cannot be described directly in QUBO form [2110.00226]. The authors therefore introduce three conversion rules from the constrained integer CDM to binary strings and then apply BOCS or Walsh-based variants to the resulting binary domain [2110.00226]. The search space is combinatorial: for $p_e=4$ and $E=4$, the number of CDMs is $15^4 = 50{,}625$, using the Bell-number count $\big(B(p_e)\big)^E$ [2110.00226]. Across these experiments, “BOCS and its variant slightly outperform the random search,” but the gains are modest [2110.00226].

The same application study gives a concrete explanation for modest improvements. It introduces the local roughness measure
$$
\Delta E(x_0)=\frac{1}{N}\sum_{i=1}^N \big(S(x_0)-S(x_1^i)\big)^2,
$$
where $x_1^i$ are Hamming-distance-1 neighbors, and finds that the vehicle-design landscape is flatter than a random SK spin-glass model [2110.00226]. The paper argues that such flatness hinders regression and acquisition optimization, reducing the advantage over random search [2110.00226]. A plausible implication is that BOCS benefits most when pairwise structure is informative and the local landscape is sufficiently non-flat to be learned from sparse samples.

## 7. Position within combinatorial Bayesian optimization

Within the broader literature on discrete black-box optimization, BOCS occupies a specific methodological niche. It is neither a Gaussian-process combinatorial BO method nor a purely heuristic combinatorial optimizer. Its core contribution is a sparse quadratic Bayesian surrogate combined with Thompson sampling and discrete acquisition optimization [1806.08838]. This gives it a characteristic profile: strong small-data behavior, direct compatibility with QUBO and Ising formulations, and explicit leverage of sparsity assumptions [2605.10856].

That QUBO compatibility has practical consequences. The vehicle-design study emphasizes that keeping the surrogate quadratic “opens the possibility” of using D-Wave quantum annealers, because QUBO is a standard form for combinatorial optimization, even when the original engineering objective is not itself quadratic [2110.00226]. The same study discusses BOCS-SA, BOCS-SDR, and BOCS-QA as solver variants, although only SA is used in the reported experiments [2110.00226].

At the same time, later research clarifies where BOCS can fail or require augmentation. If the objective is well approximated by sparse pairwise interactions, the parametric structure can be highly sample efficient [2605.10856]. If higher-order interactions dominate, the quadratic surrogate can underfit [2605.10856]. If acquisition optimization is the bottleneck, PSR improves scalability and solution quality over SDP [2008.08177]. If posterior concentration causes duplicate proposals, GP-Hedge-based candidate generation can restore progress [2605.10856]. If the search space is very large, highly categorical, or sequence-structured, relaxation-based methods such as BVO may offer better scalability and modeling flexibility [2011.02004].

BOCS is therefore best understood not as a generic label for all combinatorial Bayesian optimization, but as a particular family of methods centered on sparse quadratic pseudo-Boolean surrogates. Its continuing role in the literature is twofold: as a baseline for new discrete BO algorithms and as a practically useful method when the pairwise, sparsity-driven inductive bias is appropriate [1806.08838].

Source: https://www.emergentmind.com/topics/bayesian-optimization-of-combinatorial-structures-bocs