Combining Constraint Solving and Bayesian Techniques for System Optimization

1 Introduction

Bayesian optimization (BO) [Mockus, 1975; Frazier, 2018] is a popular technique for optimizing an objective function f , which is done through searching for input points x such that f (x) approximates the maximum max f . These solutions are presumed to be near-optimal but BO solvers do not provide formal guarantees on the accuracy and it may well happen that the solution is arbitrary far from the real solution. BO does not require an explicit representation of f, and instead search for near-optimal points is based on sampling the values of f at a limited number of input points. BO is therefore used mainly for optimizing functions whose evaluation is expensive. In particular, BO is often used for hyper-parameter optimization when training machine learning models [Snoek et al., 2012]. During the last few decades, BO has been used extensively for designing engineering systems [Mockus, 1989].

BO iteratively builds a statistical model of f(x) usually from a prior distribution defined by a Gaussian process [Rasmussen and Williams, 2006]. At each iteration i the current model is used to select the most promising candidate point xi to evaluate f(xi). This evaluation is used to update the posterior belief of the model. This process is repeated until some bound on the number of iterations is reached.

BO behaves well in practice and usually can find near-optimal points in just a few iterations. Nevertheless BO has its limitations. In particular, since BO is based on statistical approximations, there are no formal guarantees that the result achieved after a finite number of iterations is actually optimal or even close to optimal. Another limitation of BO is that even if the found point is near optimal, the solution may not be stable: there may be close points on which the value of the objective function is very different from the optimum. Such solutions are undesirable in many applications which require regions around solutions to be also near optimal. This motivates the combination of BO with SMT-based constraint solving to achieve safe and stable near-optimal solutions which we develop in this paper.

In this work we develop optimization techniques for systems modeled using real-valued functions. In many cases such functions are approximated using neural networks which are also covered in our framework. In this context, the assumption that the objective function is very expensive to sample is no longer valid. Instead, achieving tighter approximations to the maximum maxf becomes more important even if computing near-optimal points that further improve the approximation accuracy becomes significantly more expensive. Besides, in this context, it is important to guarantee that the computed near-optimal point x satisfies the design constraints and it is robust [Beyer and Sendhoff, 2007] in the sense that small perturbations to x still yield legal nearoptimal points for the objective function f.

We propose an optimization algorithm GEAROPTδ-BO which combines SMT-based constraint solving and BO and has the following properties:

• Safety: The computed near-optimal points are feasible (satisfy design constraints). In BO the feasibility of computed near-optimal solution is achieved by sampling the objective function and evaluating the constraints [Gardner et al., 2014]. In contrast, in our approach we use the explicit representation of the constraints and of the objective function to guide the search for a feasible near-optimal point.
• Accuracy: Our algorithm can find safe near-optimal points for which the value of the objective function f is within a predefined distance from the real maximum, maxf. This is enabled by the fact that in our problem setting the objective function and the constraints are given explicitly, thus we can utilize SMT solving to precisely analyse them, whereas purely probabilistic methods cannot provide full guarantees due to a limited number of sampling of f.
• Stability: Furthermore, safe near-optimal points computed by our algorithm are stable in the sense that after

perturbation of inputs within user-specified regions:

1. the safety constraints remain valid; and
2. the output remains within the near-optimal range.

Stability is a critical requirement from analog devices because the inputs and the output can be perturbed due to uncertainties in the environment such as uncertain operating conditions, design parameter tolerances or actuator imprecisions [Beyer and Sendhoff, 2007]. This concept is studied in the BO setting as robustness [Bogunovic et al., 2018; Sanders et al., 2019; Cardelli et al., 2019], where it can be estimated with high confidence but cannot be proven formally.

The paper is organized as follows. In the next section we introduce the notation. In Section 3 we recall the basics of BO and in Section 4 we define the stability and accuracy requirements for near-optimal solutions. In Section 5, we present our optimization algorithm GEAROPTδ-BO, which combines the strengths of the constraint solving algorithm GEARSATδ [Brauße et al., 2020] and reasoning with probabilities (Bayesian inference). We evaluated GEAROPTδ-BO on industrial examples coming from optimization of microprocessor design at Intel where stability is an essential requirement. We conclude in Section 7.

2 Preliminaries

Given a function f : A → B, by domf we denote its domain A. Vectors (x1, . . . , xn) may occur in the abbreviated form x. Given a, b ∈ Rn with ai ≤ bi for all i, by [[a, b]] we denote the Cartesian product×i [ai, bi] of their componentwise closed intervals [ai, bi].

In this paper we consider formulas over ⟨R, 0, 1,F, P⟩, where P are the usual order predicates <,≤,=, etc. and F contains addition, multiplication with rational constants and can also contain non-linear functions supported by SMT solvers including polynomials, transcendental functions such as combinations of sine, cosine, exponentials, solutions of differential equations and more generally computable functions [de Moura and Bjørner, 2008; Brauße et al., 2019; Brauße et al., 2021; Cimatti et al., 2013; Gao et al., 2012]. We extend functions F by functions definable by formulas: FD, i.e., we assume f ∈ FD is represented by a formula F(x1, . . . , xn, y) over variables x1, . . . , xn corresponding to the n inputs and y corresponding to the output f(x1, . . . , xn). We assume that satisfiability of such formulas is decidable or more generally δ-decidable [Gao et al., 2012]. Let us note that even when basic functions F contain just linear functions, FD will contain, e.g., functions represented by neural networks with ReLU activation functions. Optimization of functions represented by neural networks is one of the main motivations behind this work.

Throughout, x, y denote variables in formulas while a, b, c, d, e, T stand for rational constants; both forms may be indexed. Whenever we use a norm ∥·∥, we refer to the Chebyshev norm (x1, . . . , xn) ↦→ max{|x1|, . . . , |xn|}.

3 Bayesian Optimization

A basic form of BO algorithm involves two primary components: a method for statistical inference, typically Gaussian process regression [Rasmussen and Williams, 2006]; and an acquisition function for deciding where to sample next. The most popular choices for an acquisition function include expected improvement towards the optimum, knowledge gradient, and entropy search. BO iteratively builds a statistical model of an objective function f(x) from a prior distribution defined by a Gaussian process. At each iteration i the current model is used to select the most promising candidate point xi to evaluate the objective function f(xi). This evaluation is used to update the posterior belief of the model. This process is repeated until bound MaxIter on the number of iterations is reached.

The basic optimization algorithm for maximizing a blackbox function f : Rn → R, depicted in Algorithm 1, which shows a BO solver Amax with the following interface.

INIT Inputs: a, b defining bounds [[a, b]] ⊆ domf, and a vector of initial points (xi, yi)i satisfying yi = f(xi) for all i; Output: an initialized Bayesian optimizer for f with the posterior probability distribution on f updated using the available samples (xi, yi)i.

SUGGEST Inputs: none. Output: A maximizer z ∈ [[a, b]] of the acquisition function over x, where the acquisition function is computed using the current posterior distribution.

OBSERVE Inputs: (x, y). Output: the Bayesian optimizer with updated posterior probability based on y = f(x).

Algorithm 1 Basic optimization algorithm using the BO solver Amax for maximizing f(x).

MaxIter - a bound on the number of sampling iterations Sample f at m input points based on a heuristic: Let xi ∈ [[a, b]] and yi = f(xi) for i = 1, . . . ,m A ← Amax .INIT(a, b, (xi, yi)i) for j = m+ 1, . . . ,MaxIter do xj ← A.SUGGEST Sample yj = f(xj) A ← A.OBSERVE(xj , yj ) end yk ← max{y1, . . . , yMaxIter} return a near-optimal solution: (xk, yk)

We will also use a minimizing BO solver Bmin that has the same interface. This basic BO algorithm does not take as input any constraints on the inputs of f thus it is not concerned with feasibility of the computed near-optimal point. Furthermore, this algorithm is not concerned with robustness of the computed near-optimal point either and does not give any estimates of how far from the real maximum the computed near-optimal solution is. We use this basic BO algorithm, complemented with an SMT procedure where needed, as a building block of a more comprehensive, hybrid optimization algorithm that gives formal guarantees that the computed near-optimal points are feasible, stable, with values close to the real maximum up to arbitrary given accuracy. We will define stability and accuracy formally in the next section.

4 Stability and Accuracy

Given a real-valued function f on a bounded space X ⊂ Rn , we are interested in regions of X defined by the stability guard θ which we assume is a reflexive binary relation over domf. In particular, we address the problem of maximizing the minimum value of f in these regions, as stated formally in Definition 1.

The region defined by θ(x, ·) could for example be a ball with a fixed radius around x ∈ X or a radius relative to ∥x∥. It could also take a more complicated shape and we do not impose any restrictions on it except that of being definable by a quantifier-free formula θ. In the following, we assume f is in FD, i.e., defined by a formula F(x1, . . . , xn, y). We also assume that X can be defined by safety constraints.

Define the minimal value of f in the θ region of x ∈ X as

Definition 1. Given a stability guard θ for a definable and continuous f : Rn → R and a compact set X ⊆ domf, the problem

is called max-min optimization under stability guard θ. Consider ε > 0. We refer to ˜x ∈ X as a solution to this problem with accuracy ε, or ε-solution, if ˜ y ≤ y∗ < ˜ y+ε holds where y∗ = �x∗�θ,f for a solution x∗ to (1), and ˜ y = �˜x�θ,f .

In a nutshell, the value ˜ y in the definition of an ε-solution ˜x to the optimization problem under a stability guard θ is an approximation of the exact solution with guaranteed accuracy ε. It provides a lower bound for all the function's values in the stability region around ˜x that is no further away from the maximum value y∗ than ε. When f, X and θ are clear from the context, we will call ˜x just a (θ-)stable ε-solution for short. In the next section, we present algorithms based on combination of BO and SMT that solve this problem.

5 Optimisation Procedure

The problem (1) is equivalent to:

where

Let T denote a candidate bound on the value y∗ of a solution x∗ to (1). Formulas of the form G(T) with a fixed value T are also known to be in the GEAR-fragment of ∃∗∀∗ formulas [Brauße et al., 2020]. In the case we have safety constraints we can conjunctively adjoin them to the condition y ≤ y′ in Equation (3). Such additional constraints do not affect our algorithms and we will omit them for clarity of the exposition.

We first recall the decision procedure GEARSATδ [Brauße et al., 2020] shown in Algorithm 2, which we enhance with BO solvers below. It takes a potential bound T and either verifies T to be a lower bound on the optimum y∗ or proves it to be an upper bound on y∗ . GEARSATδ alternates two phases: search for candidate solution and search for counterexample for stability around the candidate solution. It does that by first finding a point x for which f(x) ≥ T holds - a candidate for a stable, safe, and accurate solution. If there is none, clearly T is an upper bound on y∗ . Otherwise, it checks whether when x is seen as the center of the stability region θ(x, ·), there is a counter-example x′ , that is, θ(x, x′ ) holds but f(x′ ) ≥ T does not (represented by Di(x′ , y′ ) constraint in Algorithm 2). In the case when there are no counter-examples with this property, we can be sure that θ(x, x′ ) implies f(x′ ) ≥ T for all x′ ∈ X. Otherwise, x′ is a counter-example and the algorithm excludes the region θ(·, x′ ) around it from the search for the next candidate for the same bound T. The stability condition θ guides the proof search by generating lemmas excluding regions around counter-examples.

Here, δ ≥ 0 refers to a constant which is used to ensure that solutions have a distance of at least δ from regions defined by θ containing counter-examples for bound T. This is done by learning lemmas of the form ¬θδ(x, d) for each counterexample d, where θδ is the δ-relaxation of θ defined by

Termination and δ-completeness of GEARSATδ was shown in [Brauße et al., 2020]. We can use GEARSATδ to find an optimal value with accuracy ε by a binary search of lower and upper bounds on the optimal value of the objective function until T ≤ y∗ < T + ε. Then ˜ y = T is the minimal value �˜x�θ,f of an ε-solution ˜x.

The key search parts in GEARSATδ rely on finding points (either candidates or counter-examples). One way to achieve this is by using an SMT solver to find points satisfying corresponding constraints (as done in [Brauße et al., 2020]) but this can be computationally expensive. In this work we propose to delegate the search part to BO and the certification part to SMT checks. In this way we take the best from both worlds: efficient search in complex spaces from BO and formal guarantees on the optimization results from SMT.

Algorithm 2 (GEARSATδ) General procedure for deciding whether a constant T is a lower or upper bound on the solution y∗ of (2).

N ←∅ ▷ known counter-examples and upper bounds

end

Our algorithms do not depend on particular types of BO

and SMT solvers used. We only assume that the SMT solver supports quantifier-free fragment including formulas F and θ.

First we integrate BO into GEARSATδ for searching counter-examples. We note that whenever θδ(ci, x′ ) is equivalent to membership in the Cartesian product [[ai, bi]] (as is always the case in our application), BO lends itself well to implement the search for a counter-example by solving Di(x′ , y), as is shown in Algorithm 3. This works by starting a new minimizing BO search in [[ai, bi]] once Algorithm 2 found a candidate ci. Note that variables N, P, ci, di and T are global across all algorithms in this section.

Algorithm 3 Finding counter-examples: Solving Di(x′ , y′ ) with Bayesian optimization and SMT in the i-th iteration of Algorithm 2.

Let (ai, bi) denote the bounds on x′ implied by θ(ci, x′ ) in Di(x′ , y′ ) and let T denote the bound on y′ in Di(x′ , y′ ) Let e1, . . . , ek ∈ [[ai, bi]] ▷ BO if f(ej) < T for some j ∈ {1, . . . , k} then di ← ej return (di, f(di)) elseBi ← Bmin .INIT(ai, bi, (ej , f (ej ))j ) for j = 1, . . . ,MaxIter do di ← Bi.SUGGEST zi ← f(di) Bi ← Bi.OBSERVE(di, zi) if zi < T then return (di, zi) fi end if Di(x′ , y′ ) is unsat then return unsat fi ▷ SMT return (di, zi) ← solution of Di(x′ , y′ ) fi

Next, we integrate BO for finding candidate solutions. For this we need to guide BO to generate points outside of regions excluded by generated lemmas. Even though most BO solvers do not support constraints like the lemmas ¬θδ(x, dm) learned by Algorithm 2, when X ⊆ domf is equivalent to membership in [[a, b]] for some a, b, it is still possible to use them in the search for a candidate c in our setting. We achieve this by penalizing any suggestion ci that satisfies θδ(x, d), for any counter-example d, with the minimal value of a counter-example found when generating the lemma. This is shown in Algorithm 4. Alongside the SMT solver it maintains the maximizing BO solver Ai that is initialized externally with selected points which are either generated in previous iterations or randomly.

Note that due to penalizing any of Ai's suggestions inside regions containing counter-examples, the function that Ai observes is not necessarily f. Whenever Ai's suggestion ci lies in a region around a previous counter-example d in the set of known ones N, that is, θ(ci, d) holds, we make Ai believe that the value of the function it optimizes is the minimum value of all counter-examples around ci instead of f(ci). Since d was a counter-example, specifically f(d) < T, this has the effect of penalizing the suggestion ci since d has already been proved to be a counter-example in the region defined by θ(ci, ·).

Algorithm 4 Finding candidates: Solving Ci(x, y) with Bayesian optimization and SMT in the i-th iteration of Algorithm 2. Ai maximizes.

if i > 1 then ▷ record previous counter-example Ai ← Ai-1.OBSERVE(ci-1, f (di-1)) fi for j = 1, . . . ,MaxIter do ▷ BO ci ← Ai.SUGGEST U ← {(d, z) ∈ N : θ(ci, d) ∧ z < T} if U ̸= ∅ then ▷ ci is excluded by lemmas yi ←min{z : (d, z) ∈ U} else ▷ ci is an eligible candidate yi ← f(ci) if yi ≥ T then return (ci, yi) fi ▷ f(ci) ≥ T N ← N ∪ {(ci, yi)} fiAi ← Ai.OBSERVE(ci, yi) ▷ yi < T end if Ci(x, y) is unsat then return unsat fi ▷ SMT return (ci, yi) ← solution of Ci(x, y)

We call GEARSATδ with integrated BO for counterexamples and candidate search GEAROPTδ-BO. This procedure is shown in Algorithm 5. The following theorem states the correctness of GEAROPTδ-BO.

Theorem 1. For any accuracy ε > 0, any stability guard θ and δ > 0, GEAROPTδ-BO is a sound, δ-complete and terminating procedure for the problem of finding safe, optimal and θ-stable ε-solutions.

Proof. (Sketch) The proof follows from the properties of GEARSATδ (Algorithm 2) which was shown to be a sound and δ-complete decision procedure for finding safe and stable solutions [Brauße et al., 2020]. δ-completeness means that if GEARSATδ terminates with unsat, then around any potential solution there is a counter-example in its θδ region. By the δ-completeness, the binary search performed by Algorithm 5 terminates with a solution to the max-min optimization problem. Since all BO enhancements are complemented by SMT checks these can be seen as an additional guidance that does not affect the correctness properties of the algorithm.

One of the distinguishing features of GEAROPTδ-BO is the exchange of candidates and counter-examples between BO and SMT solvers. In particular, when BO finds a counterexample (Algorithm 3), the generated symbolic lemma is used (a) by the SMT solver to exclude regions around this counterexample (Algorithm 2), and (b) in BO for penalizing points in excluded regions (Algorithm 4). In the other direction, candidate solutions found by the SMT solver are used to initialize the BO solver (update of P in Algorithm 5). Also, failed attempts of the BO solver to find candidates are recorded (update of N in Algorithm 4) in order to steer it away from unstable regions. As shown in Section 6, this exchange of candidates and counter-examples between BO and SMT solvers proved to be critical for both achieving tighter bounds and making runtimes feasible in our experiments.

Algorithm 5 (GEAROPTδ-BO) General procedure for optimizing f : Rn → R represented by predicate F under stability conditions θ. Arbitrary candidate lower and upper bounds l0, u0 such that l0 < u0.

P ←∅ ▷ known candidates and lower bounds N ←∅ ▷ known counter-examples and upper bounds l ← -∞ ▷ lower bound u ← +∞ ▷ upper bound loop if u = +∞ then ▷ binary search for upper bound (T, u0) ← (u0, 2u0 - l0) else if l = -∞ then ▷ binary search for lower bound (T, l0) ← (l0, 2l0 - u0) else ▷ binary search for optimum T ← (l + u)/2 fiF1(x, y) ← F (x, y) A1 ← Amax .INIT(a, b, P ∪ N) for i = 1, 2, 3, . . . do Ci(x, y) ← Fi(x, y) ∧ y ≥ T if Algorithm 4 returns unsat on Ci(x, y) then u ← T break fi (ci, yi) ← solution of Ci(x, y) Di(x′ , y′ ) ← θ(ci, x′ ) ∧ Fi(x′ , y′ ) ∧ y′ < T if Algorithm 3 returns unsat on Di(x′ , y′ ) then l ← T ▷ ci is solution to bound T P ← P ∪ {(ci, yi)} break fi (di, zi) ← solution of Di(x′ , y′ ) N ←N ∪ {(di, zi)} ▷ di is a counter-example Fi+1(x, y) ← Fi(x, y) ∧ ¬θδ(x, di) end if l + ε > u then break fi ▷ false if l or u is ±∞ end loop return {(x, l) : (x, y) ∈ P ∧ l ≤ y}

6 Benchmarks

GEAROPTδ-BO is implemented in the solver called SMLP (https://github.com/fbrausse/smlp). As Bayesian optimizers Amax ,Bmin , we used the SKOPT implementation [Pedregosa et al., 2011] based on Gaussian processes, with acquisition function gp hedge. The SMT part is implemented using the state of the art solver Z3 [de Moura and Bjørner, 2008]. Just like the implementation of GEARSATδ in [Brauße et al., 2020], besides the constraints associated with neural networks, our implementation of GEAROPTδ -BO in SMLP supports constraints over reals, integers and finite sets; this enables SMLP to optimize systems with numerical and categorical variables. Furthermore, SMLP supports optimization of f(x1, . . . , xn) by optimizing over only part of the variables x1, . . . , xk while the remaining variables xk+1, . . . , xn are treated as free inputs; thus safety, accuracy and stability of the selected near-optimal solutions for x1, . . . , xk are guaranteed for all legal input combinations of xk+1, . . . , xn as defined by the safety constraints. Transformation of GEAR formulas enabling this is explained in [Brauße et al., 2020, Section V ].

Let us note that there is a tradeoff between the size of the stability regions and proximity of the value of the solution to the pointwise optimum. One can iteratively shrink the stability regions to get arbitrary close to the pointwise optimum.

We evaluated GEAROPTδ-BO on 6 industrial examples coming from the Electrical Validation Lab at Intel. These are neural network models representing signal integrity of transmitters and receivers of a channel to a peripheral device. This application requires solutions to be safe and stable (see, Section 4), moreover the radii of stability regions are required to be proportional to the value of their respective centers. We evaluated all 4 combinations with and without BO-guided searches for candidates and counter-examples, respectively. The accuracy was set to ε = 0.05. These results are shown in Table 1. In the left-most column i refers to the problem instance, c to whether BO search was used for candidates and d to whether BO search was used for counter-examples.

6.1 Optimality under Stability Guards

Throughout our experiments, the combination of BO with SMT solvers proved to find the best bound ˜ y = 1 to the optimum, whereas SMT alone timeouts in many cases. This can be attributed to the facts that 1) when BO failed to find counter-examples (nsa), none did exist (Nsa), and that 2) the average time taken to find a counter-example is much shorter for BO (tce/nce) than for SMT (Tce/Nce for Nce ̸= 0 results). This suggests that BO constitutes a very good heuristic for finding counter-examples. On the other hand, there are many failed candidates suggested by BO (ncai in *:1:*) which indicates that BO alone is not able to find stable candidates and SMT is required to guarantee stability of the solutions.

The combinations *:0:1 correspond to those where BO tries to refute stability of SMT candidates. Throughout, it manages to do that with on average ncci/nce < 3 iterations in Algorithm 3. On the other hand, for *:1:1 when BO tries to find counter-examples to BO candidates, this average with value ≈ 8.8 is much higher on average in our experiments.

This suggests that the quality of BO candidates when guided by BO counter-examples (*:1:1) compares favorably to that of the SMT candidates with BO counter-examples (*:0:1).

6.2 Stability

Figure 1 shows the dependencies between the stability radius and the bound on the optimum achievable on candidate points. We can observe that initial candidates x found by BO are not necessarily stable (x′ close to x, but objective values differ considerably wrt. the accuracy ε), and therefore usage of SMT is essential to discover stable solutions.

6.3 Performance

From the timings in Table 1 we can see that the best bounds and timings are achieved when SMT is combined with BO for both candidate and counter-example search. When SMT is combined with BO time spent in the SMT solver (T∗) is considerably reduced in many cases by two orders of magnitude. On the other hand, the time spent by BO is also considerable increases in these cases. It turns out that as long as the total number of candidates found by BO-solver A remains low (ncap < 50), the asymptotically cubic complexity of the Gaussian process appears to be negligible compared to the overhead of rigorously solving the existential candidate formula by SMT. On the other hand, this advantage vanishes quickly after that, which motivates employing a kind of restart process similar to that successfully practiced by current SAT and SMT solvers - at least for the BO solver A. Initial experiments suggest that training samples for the restarted BO solver require careful selection.

All in all, we can see BO and SMT complement each other to solve the problem stated in Definition 1. BO's ability to rapidly produce candidates and counter-examples initially allows the combination to proceed to optimal regions quickly in most cases while still providing the formal guarantees on the validity and accuracy of stable solutions.

7 Conclusions

We introduced a hybrid optimization algorithm GEAROPTδBO that uses Bayesian optimization and SMT solvers as its building blocks. SMT solving is used to establish formal guarantees to optimality and stability of the computed solutions; while BO is used to suggest valuable candidates towards stable near-optimal solutions and significantly speeds up the overall search. In this way we combine the strengths of both approaches: the power of statistical inference by BO to guide the search with formal guarantees provided by SMT.

To the best of our knowledge this is the first work that combines BO and SMT solving to overcome basic limitations of the BO, in particular, its inability to give formal guarantees of stability and accuracy of the computed optimum, which becomes possible to resolve in cases when the objective function is given explicitly rather than as a black-box.

We believe that the observation that BO is very good in finding counter-examples in large multi-dimensional spaces, opens up new opportunities for applying BO for counterexample generation and directing the search in multiple areas of automated reasoning and formal verification.

Acknowledgements

This research was supported by a grant from Intel Corporation.