BOCS: Bayesian Optimization for Combinatorial Problems

• BOCS is a surrogate-based Bayesian optimization method for expensive combinatorial problems, using a sparse quadratic pseudo-Boolean surrogate to model first-order and pairwise interactions.
• It integrates Bayesian linear modeling with Thompson sampling, utilizing methods such as simulated annealing and semidefinite programming for effective acquisition optimization.
• The method achieves strong small-data efficiency via sparsity-inducing priors, although its performance may decline in scenarios dominated by higher-order interactions.

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 (Baptista et al., 2018).

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 (Baptista et al., 2018). 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)})$ (Shikanai et al., 11 May 2026). 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 (Shikanai et al., 11 May 2026). 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 (Deshwal et al., 2020).

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 (Baptista et al., 2018). This parametric structure is repeatedly identified as the reason BOCS works well when only a small amount of data is available (Shikanai et al., 11 May 2026).

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 $y = f(x) + \epsilon$0 (Shikanai et al., 11 May 2026). Equivalent matrix forms such as $y = f(x) + \epsilon$1 also appear, and an affine mapping to Ising variables $y = f(x) + \epsilon$2 is standard (Shikanai et al., 11 May 2026).

In the original BOCS paper, the data likelihood is Gaussian:

$y = f(x) + \epsilon$3

where $y = f(x) + \epsilon$4 is the design matrix built from constant, linear, and pairwise features (Baptista et al., 2018). The number of coefficients is $y = f(x) + \epsilon$5, so BOCS uses a sparsity-inducing horseshoe prior to control variance when $y = f(x) + \epsilon$6 (Baptista et al., 2018). 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 (Shikanai et al., 11 May 2026). 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 (Koshikawa et al., 2021).

The original paper gives a fully specified hierarchical prior:

• $y = f(x) + \epsilon$7,
• $y = f(x) + \epsilon$8,
• $y = f(x) + \epsilon$9,
• $\epsilon \sim \mathcal{N}(0,\sigma^2)$0,

together with the Makalic–Schmidt auxiliary-variable reparameterization and conjugate Gibbs conditionals (Baptista et al., 2018). The same paper emphasizes that exact sampling of $\epsilon \sim \mathcal{N}(0,\sigma^2)$1 can be implemented using the sampler of Bhattacharya et al. (2016), giving $\epsilon \sim \mathcal{N}(0,\sigma^2)$2 complexity when $\epsilon \sim \mathcal{N}(0,\sigma^2)$3 (Baptista et al., 2018).

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 (Baptista et al., 2018).

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

$\epsilon \sim \mathcal{N}(0,\sigma^2)$4

with typical penalties $\epsilon \sim \mathcal{N}(0,\sigma^2)$5 or $\epsilon \sim \mathcal{N}(0,\sigma^2)$6, which coincide on $\epsilon \sim \mathcal{N}(0,\sigma^2)$7 (Baptista et al., 2018). In the minimization form used in later work, one samples $\epsilon \sim \mathcal{N}(0,\sigma^2)$8 from $\epsilon \sim \mathcal{N}(0,\sigma^2)$9 and computes

$x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$0

which is again a QUBO (Shikanai et al., 11 May 2026).

With fixed sampled coefficients, the acquisition reduces to a BQP. In one standard maximization form,

$x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$1

where $x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$2 for $x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$3 and $x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$4 (Baptista et al., 2018). 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^\star \in \arg\max_{x \in \mathcal{D}} f(x)$5 to $x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$6, introducing an auxiliary bias variable, and constructing a lifted matrix $x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$7 so that the acquisition becomes

$x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$8

This is relaxed to the semidefinite program

$x^\star \in \arg\max_{x \in \mathcal{D}} f(x)$9

followed by randomized hyperplane rounding (Baptista et al., 2018). The paper also presents BOCS-SA, which replaces the SDP with simulated annealing over Hamming-1 neighborhoods (Baptista et al., 2018).

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 (Shikanai et al., 11 May 2026). In the vehicle-design study, the acquisition objective is explicitly written as

$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$0

and optimized by simulated annealing (Koshikawa et al., 2021).

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 (Shikanai et al., 11 May 2026).

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 (Shikanai et al., 11 May 2026). At the same time, the number of coefficients is still $x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$1, and acquisition-function optimization remains a hard combinatorial problem (Baptista et al., 2018).

The original SDP-based acquisition has substantial computational cost. Later analyses summarize generic SDP solvers for this relaxation as having worst-case time $x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$2, which becomes prohibitive as $x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$3 grows, especially when categorical variables are represented by one-hot encoding and the effective binary dimension increases (Deshwal et al., 2020). 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 (Deshwal et al., 2020).

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 (Shikanai et al., 11 May 2026). The 2026 hybrid study uses fully connected HUBO black-box functions,

$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$4

precisely to probe this mismatch, while keeping the BOCS surrogate quadratic (Shikanai et al., 11 May 2026).

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” (Shikanai et al., 11 May 2026). There the operational criterion is

$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$5

The same work attributes duplicates to posterior concentration, degeneracy among local minima in discrete spaces, and surrogate underfitting of higher-order structure (Shikanai et al., 11 May 2026).

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 (Deshwal et al., 2020).

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 (Deshwal et al., 2020). Starting from the minimization-form BQP

$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$6

PSR decomposes $x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$7, lower-bounds the positive part by an affine function

$x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$8

with $x^\star \in \arg\min_{x \in \{0,1\}^d} f(x)$9, and then solves

$D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$0

exactly by an $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$1–$D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$2 min-cut for fixed $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$3 (Deshwal et al., 2020). The outer problem maximizes this lower bound over $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$4 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 (Deshwal et al., 2020).

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 (Shikanai et al., 11 May 2026). In that hybrid, the GP uses a Hamming-distance kernel

$D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$5

forms a portfolio of Lower Confidence Bound acquisitions

$D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$6

and adaptively selects among them with GP-Hedge (Shikanai et al., 11 May 2026). On fully connected QUBO and HUBO benchmarks with $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$7, 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% (Shikanai et al., 11 May 2026).

Another important comparison is with Bayesian Variational Optimization (BVO), which is presented as a scalable Bayesian optimization method tailored to combinatorial domains (Wu et al., 2020). 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 (Wu et al., 2020). 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 (Wu et al., 2020). 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 (Baptista et al., 2018). On binary quadratic programming with $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$8, BOCS-SDP and BOCS-SA substantially outperform EI, SA, OLS, Random Search, and an MLE-based variant; with correlation length $D_t = \{(x^{(i)}, y^{(i)})\}_{i=1}^t$9 and $y^{(i)} = f(x^{(i)})$0, the average simple regret after 100 iterations is reported as approximately $y^{(i)} = f(x^{(i)})$1 $y^{(i)} = f(x^{(i)})$2 for EI and $y^{(i)} = f(x^{(i)})$3 $y^{(i)} = f(x^{(i)})$4 for both BOCS-SA and BOCS-SDP (Baptista et al., 2018). On sparsification of Ising models, BOCS-SDP achieves the best average values and lowest variance across several $y^{(i)} = f(x^{(i)})$5 values, and on contamination control it consistently finds the best or near-best schedules (Baptista et al., 2018). 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 (Baptista et al., 2018).

Later work continues to use BOCS as a primary baseline. In the BVO paper’s Ising sparsification benchmark with a $y^{(i)} = f(x^{(i)})$6 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 (Wu et al., 2020). In the pest-control benchmark with 5 categories at each of 21 stages, BVO matches COMBO and is better than RS, TPE, and SA (Wu et al., 2020). 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 (Koshikawa et al., 2021). 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 (Koshikawa et al., 2021). The search space is combinatorial: for $y^{(i)} = f(x^{(i)})$7 and $y^{(i)} = f(x^{(i)})$8, the number of CDMs is $y^{(i)} = f(x^{(i)})$9, using the Bell-number count $x \in \{0,1\}^d$0 (Koshikawa et al., 2021). Across these experiments, “BOCS and its variant slightly outperform the random search,” but the gains are modest (Koshikawa et al., 2021).

The same application study gives a concrete explanation for modest improvements. It introduces the local roughness measure

$x \in \{0,1\}^d$1

where $x \in \{0,1\}^d$2 are Hamming-distance-1 neighbors, and finds that the vehicle-design landscape is flatter than a random SK spin-glass model (Koshikawa et al., 2021). The paper argues that such flatness hinders regression and acquisition optimization, reducing the advantage over random search (Koshikawa et al., 2021). 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 (Baptista et al., 2018). This gives it a characteristic profile: strong small-data behavior, direct compatibility with QUBO and Ising formulations, and explicit leverage of sparsity assumptions (Shikanai et al., 11 May 2026).

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 (Koshikawa et al., 2021). The same study discusses BOCS-SA, BOCS-SDR, and BOCS-QA as solver variants, although only SA is used in the reported experiments (Koshikawa et al., 2021).

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 (Shikanai et al., 11 May 2026). If higher-order interactions dominate, the quadratic surrogate can underfit (Shikanai et al., 11 May 2026). If acquisition optimization is the bottleneck, PSR improves scalability and solution quality over SDP (Deshwal et al., 2020). If posterior concentration causes duplicate proposals, GP-Hedge-based candidate generation can restore progress (Shikanai et al., 11 May 2026). 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 (Wu et al., 2020).

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 (Baptista et al., 2018).