Scalable Determination of Penalization Weights for Constrained Optimizations on Approximate Solvers

The gist: This work introduces a novel, broadly applicable algorithm for determining penalization terms for approximate solvers in constrained optimization problems, proving that it can yield an efficient strategy that approximates the minimal penalization constant M with provable guarantees

Problem Formulation and the Big-M Challenge

Combinatorial optimization problems are often cast into a Quadratic Unconstrained Binary Optimization (QUBO) problem by adding quadratic penalty terms weighted by a constant M, which is known as the Big-M approach. The choice of this constant M significantly shapes the energy landscape; if M is too high, the solver prioritizes constraint satisfaction at all costs, returning infeasible states that may be far from optimal with respect to the original objective function. Existing computationally-efficient prescriptions for Big-M often overestimate the required penalty, which degrades solution quality. The paper addresses this by formalizing the problem as finding a minimal penalization constant M∗η exact that ensures an η-reformulation for a solver with guaranteed energy threshold Ef

The Proposed Strategy for Determining Penalization Weights

The core of the algorithm combines analytical considerations with uniform sampling over feasible configurations to derive efficiently evaluable bounds on the solver’s output distribution, from which M is calculated. The strategy involves calculating four specific bounds: (1.) a lower bound B<F on the probability of sampling feasible points with objective below and exceeding Ef, (2.) an upper bound B>F on the probability of feasible points with a high objective value, and (3.) an upper bound BF¯ on the probability of infeasible events. The penalization weight M is then determined as the unique root of a scalar function g(M), which depends on the targeted success probability η

Algorithm and Theoretical Guarantees

Algorithm 1 outlines the systematic strategy for determining M, involving steps such as determining an lower bound ELB from an SDP relaxation, estimating feasible spectral weights n∆(e) via uniform sampling, and computing the bounds BF, and BF¯. The correctness of this algorithm in the limit of infinite samples is established by Theorem 2, which proves that if Algorithm 1 returns M∗ ≠ <>, then (PM) with M = M∗ is an η-reformulation with guaranteed energy threshold Ef for a Gibbs sampler at inverse temperature β Theorem 4 further establishes the finite sample error control, showing that with probability at least 1 − δ, p ≥ η – ϵ when Ns ≥ 2/ϵ2δ2

Practical Applicability and Performance

The method is shown to be efficient for large classes of problems, scaling polynomially in the system size under appropriate assumptions on the problem specification and solver behavior. Numerical demonstrations across diverse problems like the Traveling Salesman Problem (TSP), Multiway Number Partitioning Problem (MNPP), and Portfolio Optimization (PO) show robust performance, including order-of-magnitude speedups over existing heuristics when benchmarking on hardware like Fujitsu’s Digital Annealer The paper demonstrates that using this method can reduce the number of solver calls required for a binary search approach by factors of 10 or more, showing practical advantages in time to solution

Parameter Sensitivity and Robustness

The determination of M is sensitive to hyperparameters like Ef and vcut, which control the energy threshold and the maximal penalization value considered for infeasible points. Analysis shows that for problems like MNPP and TSP, small values of vcut do not change M∗ accounting for most of the probability mass of infeasible points, confirming robustness For PO instances, values vcut ≥ 16 are required to capture all relevant infeasible configurations and achieve an ηeff close to the target Furthermore, the algorithm provides a tool for solving the inverse problem—determining the inverse temperature β that allows a thermal solver to sample desired solutions with probability at least η when M is fixed

Conclusion

The proposed algorithm improves on state-of-the-art general heuristics for penalization constants by providing a tool to precisely set the penalization constant so as to control the probability of success, using knowledge of the problem structure and the solver’s statistical behavior This approach is particularly relevant for quantum solvers, such as quantum annealers or QAOA, where classical computational power for pre-processing is abundant compared to the resource costs of the quantum solver The work concludes by showing that using this method can reduce the time to solution on a solver by an order of magnitude compared to strategies based on binary search for a penalization weight

--- Page 1 ---

The gist: This work introduces a novel, broadly applicable algorithm for determining penalization terms for approximate solvers in constrained optimization problems, proving that it can yield an efficient strategy that approximates the minimal penalization constant M with provable guarantees

How it works

The problem is formulated by converting linearly constrained binary optimization (LCBO) problems into a QUBO problem by promoting constraints A x = b to quadratic penalty terms weighted by a penalization constant M, resulting in the function minimize x∈G(x) = x tQx + M(Ax − b) 2 The goal is to find the minimal penalization constant M∗η exact that ensures an η-reformulation for a solver with guaranteed energy threshold Ef

A Big-M strategy

The general idea for an efficient strategy is to calculate bounds on the probability of three distinct events: (i) observing a feasible point with objective smaller or equal than Ef, (ii) a feasible point with objective larger than Ef and (iii) an infeasible point. These three bounds are combined to determine M∗ ≥ M∗η which will be closer to M∗η the tighter the bounds are

Algorithm 1

Algorithm 1 combines these estimates to compute the bounds BF, and BF¯, that bound the probabilities of observing feasible points with a low objective value, feasible points with a high objective value, and infeasible points, respectively. The final penalization weight M is determined as the unique root of the scalar function g(M), which depends on the targeted success probability η

Validation and Numerical Benchmarks

The proposed strategy was benchmarked on three classes of constrained problems: Traveling Salesman Problem (TSP), Multiway Number Partitioning Problem (MNPP) and Portfolio Optimization (PO) Experiments using an ideal Gibbs sampler, simulated annealing algorithm, and the Digital Annealer were conducted across various system sizes up to several thousand bits The results show that for the Digital Annealer, ηeff is consistently higher than the targeted η in all settings

Discussion

The algorithm improves on state-of-the-art general heuristics for penalization constants by providing a tool to precisely set the penalization constant so as to control the probability of success, using knowledge of the problem structure and the solver’s statistical behavior It provides a tool to precisely set the penalization constant so as to control the probability of success, using knowledge of the problem structure and the solver’s statistical behavior This method can yield practical advantages in time to solution by reducing overhead proportional to log2(Ml1 /M∗)

References

[1] A common approach to constrained combinatorial optimization problems is to enforce the constraints in the QUBO formulation by adding penalization terms. Penalization introduces an additional hyperparameter that significantly affects the solver’s efficacy: the relative weight between the objective terms and the penalization terms

[2] Combinatorial optimization problems arise widely in both theoretical research and real-world applications, aiming to identify the optimal configuration of discrete decision variables that minimizes a given objective function. The corresponding search space typically has a size scaling exponentially with the number of variables, rendering the optimization a computationally hard problem

[3] Practical solvers tend to be heuristic algorithms that yield approximate solutions. For instance, simulated annealing (or any solver based on Gibbs sampling) outputs a configuration sampled from a low-temperature thermal distribution over the combinatorial space, with an objective value close (but in general not equal) to the optimum

[4] Importantly, general combinatorial optimization problems, arising in applications, usually include constraints that restrict the search space. Yet, a standard approach to tackle such problems with QUBO solvers is to convert constraints into penalization terms weighted by a constant, commonly referred to as Big-M

[5] In this work, we introduce a novel, broadly applicable algorithm for a priori determining the penalization term for a given constrained optimization problem and a specified approximate solver.

1. Bold M determination for approximate solvers: The algorithm combines analytical considerations with uniform sampling over feasible configurations, to derive efficiently evaluable bounds on the solver’s output distribution, from which the penalization weight M is calculated. This allows a solver like Fujitsu’s Digital Annealer (version 3) to determine a controllable, guaranteed minimum probability of observing feasible solutions with an energy below a threshold.

2. Efficient runtime scaling: The method demonstrates that for large classes of constrained optimization problems, the algorithm’s runtime and memory scales polynomially in the system size, which is crucial for handling large-scale instances on Fujitsu’s Digital Annealer.

3. Reduced solver overhead: By using the derived optimal penalty weight, the system can achieve a reduction in time to solution by a factor of 10 or more compared to binary search methods, quantified as log2(Ml1 /M∗).

4. Adaptive constraint satisfaction control: The system can be tuned by setting the energy threshold Ef such that it is not impeding the success probability, allowing for precise control over the trade-off between solution quality and feasibility guarantee.

5. Inverse temperature determination: An extension allows determining an optimal inverse temperature β∗ as the root of β 7→ g(M) for given M, enabling thermal solvers to meet sampling requirements even when penalty weights are fixed due to hardware limitations.