{'fraction_non_alphanumeric': 0.045112648885899866, 'fraction_numerical': 0.009981770879422331, 'mean_word_length': 4.07764198418404, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 13, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Efficient MIP Techniques for Computing the Relaxation Complexity

The relaxation complexity rc(X) of the set of integer points X contained in a polyhedron is the minimal number of inequalities needed to formulate a linear optimization problem over X without using auxiliary variables. Besides its relevance in integer programming, this concept has interpretations in aspects of social choice, symmetric cryptanalysis, and machine learning. We employ efficient mixed-integer programming techniques to compute a robust and numerically more practical variant of the relaxation complexity. Our proposed models require row or column generation techniques and can be enhanced by symmetry handling and suitable propagation algorithms. Theoretically, we compare the quality of our models in terms of their LP relaxation values. The performance of those models is investigated on a broad test set and is underlined by their ability to solve challenging instances that could not be solved previously.


Introduction
Let X ⊆ Z d be such that X = conv(X) ∩ Z d and let Y ⊆ Z d . A fundamental problem in various fields is to find a polyhedron P with the minimum number of facets such that X ⊆ P and (Y \ X) ∩ P = ∅. We call this quantity the relaxation complexity of X w.r.t. Y , in formulae, rc(X, Y ), and any such polyhedron a relaxation. In case Y = Z d , we write rc(X) instead of rc(X, Z d ). In the theory of social choice, X ⊆ {0, 1} d can be interpreted as the winning strategies of a simple game, see [26,Chap. 8.3]. One is then interested in computing rc(X, {0, 1} d ), i.e., the smallest number of inequalities needed to distinguish winning and loosing strategies. In symmetric cryptanalysis, a subfield of cryptography, rc(X, {0, 1} d ) corresponds to the minimum number of substitutions in symmetric key algorithms [25]. In machine learning, relaxations P correspond to polyhedral classifiers that distinguish two types of data points [1]. The relaxation complexity is then the minimum size of a polyhedral classifier. Finally, of course, rc(X) is the minimum number of inequalities needed to formulate a linear optimization problem over X ⊆ Z d without using auxiliary variables.
Depending on the application, different strategies have been pursued for computing and bounding the relaxation complexity. For example, Kaibel & Weltge [19] introduced the notion of hiding sets for deriving lower bounds on rc(X). Using this technique, they could show that several sets X arising from combinatorial optimization problems have superpolynomial relaxation complexity. Moreover, rc(X, Y ) can be found by computing the chromatic number of a suitably defined hypergraph; deriving lower bounds on the chromatic number allowed Kurz & Napel [20] to find a lower bound on rc(X, {0, 1} d ) in the context of social choice. In machine learning, algorithms have been devised to construct polyhedral classifiers and thus providing upper bounds on rc(X, Y ), see [1,8,21,22]. To find the exact value of rc(X, {0, 1} d ) in the context of symmetric cryptanalysis, mixed-integer programming models have been investigated. For higher dimensions, however, many of these models cannot compute rc(X, {0, 1} d ) efficiently in practice.
In this article, we follow the latter line of research. Given the relevance of knowing the exact value of rc(X, Y ), our aim is to develop efficient mixedinteger programming (MIP) techniques for computing rc(X, Y ), if both X and Y are finite. More precisely, we investigate methods to compute rc ε (X, Y ), a more robust variant of rc(X, Y ) that is numerically more practical as we discuss below. To this end, we propose in Section 2 three different MIP models that allow to compute rc ε (X, Y ): a compact model as well as two more sophisticated models that require row or column generation techniques. Section 3 compares the quality of the three models in terms of their LP relaxation value, and we discuss several enhancements of the basic models in Section 4. These enhancements include tailored symmetry handling and propagation techniques as well as cutting planes. Finally, we compare the performance of the three different models on a broad test set comprised of instances with different geometric properties and instances arising in symmetric cryptanalysis (Section 5). Our novel methods allow to solve many challenging instances efficiently, which was not possible using the basic models.
We remark that the basic versions of two models have already been used by us in [2] to find rc ε (X, Y ) for X being the integer points in low-dimensional cubes and crosspolytopes. These experiments helped us to prove general formulae for rc(X) in these cases. For this reason, we believe that the more sophisticated algorithms described in this article are not only of relevance for practical applications, but also to develop hypotheses for theoretical results. Our code is publicly available at github 1 .

Related Literature
One of the earliest references on the relaxation complexity goes back to Jeroslow [15] who showed the tight bound rc(X, {0, 1} d ) ≤ 2 d−1 , for any X ⊆ {0, 1} d . This result has been complemented by Weltge [28] who showed that most X ⊆ {0, 1} d have rc(X, {0, 1} d ) ≥ 2 d c·d 3 , for some absolute constant c > 0. Moreover, hiding sets proposed by Kaibel & Weltge [19] provide a lower bound on rc(X). The bound given by hiding sets can be improved by computing the chromatic number of a graph derived from hiding sets, see [2]. Regarding the computability of rc(X), it has been shown in [3] that there exists a proper subset Obs(X) of Z d \ X such that rc(X) = rc(X, Obs(X)). If Obs(X) is finite, they show that rc(X, Obs(X)), and thus rc(X), can be computed by solving a mixed-integer program. They also provide sufficient conditions on X that guarantee Obs(X) to be finite. Moreover, they establish that rc(X) is computable if d ≤ 3; for d = 2, a polynomial time algorithm to compute rc(X) is discussed in [2]. In general, however, it is an open question whether rc(X) is computable.
One drawback of relaxations of X as defined above is that they might be sensitive to numerical errors. If a ⊺ x ≤ β is a facet defining inequality of a relaxation of X that separates y ∈ Z d \ X, then we only know a ⊺ y > β. Thus, slightly perturbing a might not separate y anymore. To take care of this, we suggested in [2] to add a safety margin ε > 0 to the separation condition. That is, if a ⊺ x ≤ β is a facet defining inequality of a relaxation of X with a ∞ = 1 that separates y, then we require a ⊺ y ≥ β+ε. In this case, we say that y is ε-separated from X. Then, rc ε (X) denotes the smallest number of facets of any relaxation of X that satisfies the safety margin condition 2 . We call such a relaxation an εrelaxation of X. Analogously to rc(X, Y ), we define rc ε (X, Y ) to be the smallest number of inequalities needed to ε-separate X and Y \ X. As ε-relaxations are more restrictive than relaxations, rc ε (X) ≥ rc(X) for each ε > 0. In contrast to rc(X), however, we show in [2] that for every finite and full-dimensional Thus, rc ε (X) is computable and the aim of this article is to develop MIP techniques that allow to find rc ε (X, Y ) efficiently. In particular, if ε approaches 0, then rc ε (X) converges towards rc Q (X), a variant of the relaxation complexity which requires the relaxations to be rational. Further variations of rc(X) in which the size of coefficients in facet defining inequalities are bounded are discussed in [10,11].
Besides finding relaxations of X, another field of research aims to find outer descriptions of P = conv(X) to be able to use linear programming techniques to solve optimization problems over X. Since P might have exponentially many facets, the concept of extended formulations has been introduced. Extended formulations are polyhedra Q ⊆ R d+k whose projection onto R d yields P . The smallest number of facets of an extended formulation of P is its extension complexity xc(P ). We refer the reader to the surveys of Conforti et al. [7] and Kaibel [18] as well as the references therein. Extended formulations that allow to use integer variables have been discussed, e.g., by Bader et al. [4], Cevallos et al. [6], and Weltge [28,Chap. 7.1]. A combination of rc(X, {0, 1} d ) and xc(conv(X)) has been studied by Hrubeš & Talebanfard [14].  [19] proved that the cardinality of any hiding set is a lower bound on rc(X). The maximum size of a hiding set is denoted by

Basic Definitions and Notation
Analogously to H(X), H(X, Y ) denotes the maximum size of a Y -hiding set.

Mixed-Integer Programming Models to Compute
In this section, we discuss three different mixed-integer programming models to compute rc ε (X, Y ). The three different MIP formulations that we discuss differ in the way how they model rc ε (X, Y ). The first model uses only polynomially many variables and inequalities, the second model needs exponentially many inequalities while the number of variables is still polynomial, and the third model requires exponentially many variables but only polynomially many inequalities. For this reason, we refer to these three models as the compact, the cutting plane, and the column generation model, respectively. In preliminary experiments with our code, we have already used the compact and column generation model [2]. Nevertheless, we provide the full details of these models to make the article self-contained and to be able to explain the model enhancements. For the sake of convenience, we assume for the remainder of this article that X and Y are disjoint. This is without loss of generality, because we can replace Y by Y \ X, which does not change the value of rc ε (X, Y ). We also refer to X as the set of feasible points, whereas the points in Y are called infeasible.

Compact Model
Observe that lattice-convex sets are exactly those subsets of Z d that admit a relaxation. In [3], a mixed-integer programming formulation has been proposed to check whether a finite lattice-convex set X admits a relaxation with k inequalities, and we have explained in [2] how to adapt the model to be able to compute rc ε (X, Y ).
Given an upper bound k on the number of inequalities needed to separate X and Y , the model's idea is to introduce variables a ij and b i , to model the k potential inequalities needed in a relaxation. Moreover, for each y ∈ Y and i ∈ [k], a binary variable s yi is introduced that indicates whether the i-th inequality is violated by y; additional binary variables u i , i ∈ [k], indicate whether the i-th inequality is needed in a relaxation. Using a big-M term with M ≥ d(ρ X + ρ Y ) + ε, with ρ X = max{ x ∞ : x ∈ X} and ρ Y = max{ y ∞ : y ∈ Y }, the mixed-integer programming formulation for rc ε (X, Y ) is as follows: Inequalities (1b) ensure that the k inequalities are valid for X and Inequalities (1c) guarantee that each y ∈ Y is cut off by at least one inequality. If an inequality is selected to separate y ∈ Y and X, Inequalities (1d) ensure that this is consistent with the k inequalities defined by the model.

Cutting Plane Model
To be able to find rc ε (X, Y ), Model (1) introduces two classes of variables: variables u and s model which inequalities are used and subsets of Y that are separated by the selected inequalities, respectively, whereas variables a and b guarantee that the subsets defined by s can be cut by valid inequalities for X. The problem of computing rc ε (X, Y ) can thus be interpreted as a two stage problem, where the first stage selects a set of subsets of Y and the second stage checks whether the selected subsets correspond to feasible cut patterns. Since the first stage variables are binary and the second stage problem is a feasibility problem, logic-based Benders decomposition can be used to compute rc ε (X, Y ), see [13]. While classical Benders decomposition requires the subproblem to be a linear programming problem, logic-based Benders decomposition allows the subproblem to be an arbitrary optimization problem. Let C = {C ⊆ Y : C and X are not linearly ε-separable}. We refer to C as the conflict set. For all (C, i) ∈ C × [k], the conflict inequality y∈C s yi ≤ |C| − 1 models that not all points in C can be cut by an inequality valid for X. Conse-quently, is an alternative model for computing rc ε (X, Y ).

Column Generation Model
Let I = {I ⊆ Y : I and X are linearly ε-separable}. Then, rc ε (X, Y ) is the smallest number ℓ of sets I 1 , . . . , I ℓ ∈ I such that Y = ℓ i=1 I i . Thus, instead of using the matrix s ∈ {0, 1} Y ×[k] to encode which inequality cuts which points from Y , we can introduce for every I ∈ I a binary variable z I ∈ {0, 1} that encodes whether an inequality separates I or not: where I y = {I ∈ I : y ∈ I}.
Remark 2.1. In contrast to Model (1), Models (2) and (3) do not directly provide an ε-relaxation of X w.r.t. Y . To find such a relaxation, rc ε (X, Y ) many linear programs need to be solved in a post-processing step.

Comparison of Basic Models
While the compact model (1) can be immediately handed to an MIP solver due to the relatively small number of variables and constraints, the cutting plane model (2) and column generation model (3) require to implement separation and pricing routines, respectively. At least for the column generation model, this additional computational effort comes with the benefit of a stronger LP relaxation in comparison with the compact model. To make this precise, we denote by v ⋆ com , v ⋆ cut , and v ⋆ CG the optimal LP relaxation value of the compact, cutting plane, and column generation model, respectively. Proposition 3.1. Let X ⊆ Z d be finite and lattice-convex, let Y ⊆ Z d \ X be finite, let ε > 0 such that rc ε (X, Y ) exists, and suppose both Models (1) and (2) are feasible.
Note that 2 is a technical assumption that is almost always satisfied in practice, e.g., to approximate rc(X, Y ) by rc ε (X, Y ), one selects Proof. First we show v ⋆ cut ≥ 1 and v ⋆ com ≥ 1. Observe that we get for every (partial) feasible solution (s, u) and everyȳ ∈ Y the estimation where k is the upper bound used in Model (1) of (2). Hence, v ⋆ cut ≥ 1 and v ⋆ com ≥ 1. If the upper bound k = 1, we thus have necessarily v ⋆ cut = 1. If k ≥ 2, we construct a feasible solution for (2) with objective value 1 by assigning all variables value 0 except for s yi , (y, i) ∈ Y × [2], u 1 , and u 2 , which get value 1 2 . Indeed, the left-hand side of each conflict inequality evaluates to |C| 2 , while the right-hand side is |C| − 1. Thus, because |C| ≥ 2 for any conflict as X is lattice-convex, the find |C| 2 ≤ |C| − 1, i.e., all conflict inequalities are satisfied. Since the remaining inequalities hold trivially, For the second statement, we assume k ≥ 2, because otherwise v ⋆ com = 1 follows as above. We define a feasible solution with objective value 1 of Model (1) by assigning all variables value 0 except for The inequalities a i· ⊺ x ≤ b i defined this way are either 0 ≤ 0 or x 1 ≤ ρ X , which are valid for X. Moreover, the Inequalities (1d) are satisfied, because for i = 1 and every y ∈ Y , we have and for the remaining i ≥ 2, we get ε − M (1 − s y2 ) = 0. Since one can easily check that the remaining inequalities of (1) are also satisfied, v ⋆ com ≤ 1 follows, concluding the proof using the first part of the assertion.
The value of the LP relaxations thus does not indicate whether the compact or cutting plane model performs better in practice. An advantage of the latter is that the conflict inequalities encode a hypergraph coloring problem, which is a structure appearing frequently in practice. Hence, there might be a chance that a solver can exploit this structure if sufficiently many inequalities have been separated. The compact model, however, might have the advantage that the aand b-variables guide the solver in the right direction when branching on sor u-variables, because feasibility is already encoded in the model and does not need to be added to the model by separating cutting planes. (1) and (2) are feasible. Let k be the number of inequalities encoded in Model (1).

and suppose both Models
1. If there exists an optimal solution of the LP relaxation of (3) that assigns at most k variables a positive value, , and this can be strict.
recall that for each I ∈ I there exists an inequality a(I) ⊺ x ≤ b(I) + ε separating I and X. Due to rescaling, we may assume that a(I) If we are given a solution z ∈ R I + of (3) with at most k non-zero entries, we define a solution of the LP relaxation of (1) with the same objective value as follows. Let I 1 , . . . , I ℓ ∈ I be the indices of non-zero entries in z. For each i ∈ [ℓ] and y ∈ Y , define For i ∈ {ℓ + 1, . . . , k} and y ∈ Y , we define s yi = 0 and Indeed, this solution adheres to (1b) since (a, b) defines valid inequalities, and also (1e)-(1g) hold trivially. By definition, s and u also satisfy the box constraints corresponding to (1h). To see that (1c) holds, note that for each y ∈ Y , since z is feasible for the LP relaxation of (3). For the last constraint (1d), note that the constraint is trivially satisfied if s yi = 0. If s yi > 0, then a i· ⊺ x ≤ b i corresponds to an inequality separating X and y, which finally shows that the newly defined solution is feasible for the LP relaxation of (1). To conclude, note com and the remaining estimations hold by Proposition 3.1.
For the second part, let H ⊆ Y be a hiding set for X and let z ∈ R I + be an optimal solution of the LP relaxation of (3). Then, for distinct y 1 , y 2 ∈ H, we have . To see that the inequality can be strict, consider X = {0, 1} 2 and let Y be all infeasible points in Z 2 with ℓ ∞ -distance 1 from X. One can readily verify that a maximum hiding set for X has size 2, while the LP relaxation of (3) has value 8 3 . If Y contains a hiding set of size at least 2, the column generation model is thus strictly stronger than the compact and cutting plane model. In particular, the gap between v ⋆ CG and v ⋆ cut (and v ⋆ com ) can be arbitrarily large: if d = 2 and Y = Obs(X), there is always a hiding set of size rc(X, Y ) − 1, see [2,Thm. 23].

Enhancements of Basic Models and Algorithmic Aspects
In their basic versions, the compact and cutting plane model are rather difficult to solve for a standard MIP solver, e.g., because not enough structural properties of rc ε (X, Y ) are encoded in the models that are helpful for a solver. Moreover, the cutting plane and column generation model require to solve a separation and pricing problem, respectively, to be used in practice. In this section, we discuss these aspects and suggest model improvements.

Incorporation of Structural Properties
In the following, we describe cutting planes, propagation algorithms, and techniques to handle symmetries and redundancies in the compact and cutting plane model.

Cutting Planes
In both the compact and cutting plane model, variable s yi encodes whether a point y ∈ Y is separated by inequality i ∈ [k]. To strengthen the compact model and the initial LP without separated inequalities in the cutting plane model, we can add inequalities that rule out combinations of points from Y that cannot be separated simultaneously. For any hiding set H ⊆ Y , the hiding set cut encodes that each inequality i ∈ [k] can separate at most one element from a hiding set. Although these cuts are the stronger the bigger the underlying hiding set, we add these inequalities just for hiding sets of size 2. The reason for this is that such hiding sets can be found easily by iterating over all pairs (y 1 , y 2 ) of distinct points in Y and checking whether the line segment conv({y 1 , y 2 }) intersects conv(X) non-trivially. In our implementation, we insert the expression λy 1 +(1−λ)y 2 in each facet defining inequality of conv(X) to derive bounds on the parameter λ. Then, the final bounds on λ are within [0, 1] if and only if {y 1 , y 2 } is a hiding set. For hiding sets of arbitrary cardinality, the task is more difficult, because there might exist exponentially many hiding sets. Thus, we are relying on a separation routine for hiding set cuts. The separation problem for hiding set cuts, however, is at least as difficult as finding a maximum hiding set for X, and the complexity of the latter is open.
Propagation Suppose we are solving the compact and cutting plane model using branch-and-bound. At each node of the branch-and-bound tree, there might exist some binary variables that are fixed to 0 or 1, e.g., by branching decisions. The aim of propagation is to find further variable fixings based on the already existing ones.
Our first propagation algorithm is based on the following observation.
The convexity propagation algorithm computes the sets , and fixes s yi to 1 for all y ∈ F ′ i . If there is y ′ ∈ F ′ i such that s y ′ i is already fixed to 0, then the algorithm prunes the node of the branch-and-bound tree. This is indeed a valid operation, because Inequalities (1c) and (2b) allow each point y ∈ Y to be separated by several inequalities.
The second propagation algorithm exploits that F i ∩ conv(X) needs to be empty in each feasible solution. The intersection propagation algorithm thus iterates over all y ∈ Y \ F i and checks whether conv If the check evaluates positively, s yi is fixed to 0.
Comparing both propagation algorithms, the convexity propagator requires to compute only a single convex hull per set F i , whereas the intersection propagator needs to compute O(|Y |) convex hulls per set F i , which can be rather expensive. In our experiments, we will investigate whether the additional effort pays off in reducing the running time drastically. To avoid computing unnecessary convex hulls, we call both propagation algorithms in our implementation only if the branching decision at the parent node is based on a variable s yi , and in this case only for this particular inequality index i and no further i ′ ∈ [k]\ {i}.

Symmetry Handling
It is well-known that the presence of symmetries slows down MIP solvers, because symmetric solutions are found repeatedly during the solving process leading to an exploration of unnecessary parts of the search space. In a solution of the compact and cutting plane model, e.g., we can permute the inequality labels i ∈ [k] without changing the structure of the solution. For this reason, one can enforce that only one representative solution per set of equivalent solutions is computed without losing optimal solutions.
One way of handling symmetric relabelings of inequalities is to require that the columns of the matrix s ∈ {0, 1} Y ×[k] are sorted lexicographically nonincreasingly. To enforce sorted columns, we use a separation routine for orbisack minimal cover inequalities as suggested in [12] and the propagation algorithm orbitopal fixing by Bendotti et al. [5]. Both algorithms' running time is in O(|Y | · k). Moreover, sorting the columns of s implies that we can also require the u-variables to be sorted, i.e., the first rc ε (X, Y ) inequalities are the inequalities defining an ε-relaxation, which can be enforced by adding to the problem. Besides the symmetries of relabeling inequalities, we might also be able to relabel points in Y without changing the structure of the problem. This is the case if we find a permutation π of [d] such that π(X) = X and π(Y ) = Y , where for a set T ⊆ R d we define π(T ) = {π(t) : t ∈ T } and π(t) = (t π −1 (1) , . . . , t π −1 (d) ). The permutation π gives rise to a permutation φ of Y and ψ of X, where φ(y) := π(y) and ψ(x) := π(x).
If there exists a permutation π of [d] such that π(X) = X and π(Y ) = Y , then also (s ′ , u) is a (partial) solution, where s ′ arises from s by reordering the rows of s according to φ.
Proof. Suppose (s, u) is a solution of Model (2). Then, (s, u) can be extended to a solution of Model (1) , such that the i-th inequality separates the points in F i = {y ∈ Y : s yi = 1} from X. If we apply permutation π to X and Y , we do not change the structure of the problem, that is, , defines also a relaxation of X w.r.t. Y . Thus, if the original i-th inequality separated point y ∈ Y , the permuted inequality separates φ(y). Consequently, if we define s ′ by relabeling the rows of s according to φ, (π(a), b, s ′ , u) is a solution of Model (1) and thus (s ′ , u) is a solution of Model (2).
If Π = {π ∈ S d : π(X) = X, π(Y ) = Y } and Φ is the group containing all φ associated with the permutations π ∈ Π, Lemma 4.2 tells us that we can also force the rows of s to be sorted lexicographically non-increasingly w.r.t. permutations from Φ. In our implementation, we compute a set Γ of generators of the group Φ and enforce for each γ ∈ Γ that matrix s is lexicographically not smaller than the reordering of s w.r.t. γ. We enforce this property by separating minimal cover inequalities for symresacks and a propagation algorithm, see [12]. Both run in O(k) time per γ ∈ Γ.
To detect the symmetries Φ, we construct a colored bipartite graph G = (V, E). The left side of the bipartition is given by X ∪ Y and the right side is defined as There is an edge between z ∈ X ∪ Y and (v, j) ∈ R if and only if z j = v. Moreover, each node gets a color uniquely determining its type: all nodes in X are colored equally with color "X", all nodes in Y are colored equally by color "Y ", and node (v, j) ∈ R is colored by color "v". Then, the restriction of every automorphism σ of G to R corresponds to a permutation in Π, and thus, restricting σ to Y is a permutation in Φ.
Note that the graph G defined above might not allow to detect symmetries if a symmetric arrangement of X and Y is translated asymmetrically. For example, if X = t + ∆ 2 , Y = t + (∆ 2 + ♦ 2 ) \ ∆ 2 , and t = 1 2 , then there is no permutation keeping X invariant. For this reason, we use in the construction of G relative coordinates. That is, for each coordinate j ∈ [d], we compute µ j = min z∈X∪Y z j and translate X ∪ Y by −µ before building G.
Another way of handling symmetries for the compact model (1) is to handle symmetries of the inequalities d j=1 a ij x j ≤ b i defined in the model. We can reorder the inequalities d j=1 a ij x j ≤ b i , i ∈ [k] that are (not) used in the relaxation, to obtain another solution with the same objective value. To handle these symmetries, we can add the inequalities Inequalities (5) sort the inequalities (not) present in a relaxation by their first coefficient. The inequalities are compatible with Inequalities (4), but not necessarily with the lexicographic ordering constraints. The latter is the case because cutting the point y ∈ Y associated with the first row of matrix s might require a very small first coefficient in any separating inequality, whereas other points might require a very large first coefficient. In our experiments, we will investigate which symmetry handling method works best for the compact and cutting plane model. (1) can be handled by enforcing that d j=1 a ij x j ≤ b i becomes the trivial inequality 0 ⊺ x ≤ dρ X if it is not used in a relaxation of X (i.e., u i = 0). This removes infinitely many equivalent solutions from the search space, and can be modeled by replacing (1f) by

Finally, additional redundancies in Model
and the lower bound constraint in (1g) by This method is compatible with both the lexicographic symmetry handling approach and Inequalities (5).

Algorithmic Aspects of the Cutting Plane Model
To be able to deal with the exponentially many conflict inequalities (2c) in the cutting plane Model (2), we are relying on a separation routine. We start by discussing the case that the point s ⋆ to be separated is contained in {0, 1} Y ×[k] , i.e., for each of the k inequalities we already know which points it is supposed to separate.
To check whether s ⋆ ∈ {0, 1} Y ×[k] satisfies all conflict inequalities, we can compute for each i ∈ [k] the set F i = {y ∈ Y : s ⋆ yi = 1}, and build a linear program similar to Model (1) that decides whether X and F i are ε-separable. If the answer is yes, we know s ⋆ is feasible. Otherwise, we have found a violated conflict inequality, namely y∈Fi s yi ≤ |F i | − 1. Of course, this inequality will be rather weak in practice, because it excludes only the single assignment F i .
One way to strengthen the inequality is to search for a minimum cardinality subset F min of F i , which cannot be separated from X. The corresponding inequality y∈Fmin s yi ≤ |F min | − 1 then does not only cut off s ⋆ , but every solution that assigns inequality i all points from F min . However, we do not expect that F min can be computed efficiently, because detecting a minimum cardinality set of inequalities whose removal leads to a feasible LP is NP-hard, see Sankaran [24]. Instead, we compute a minimal cardinality subset F ⊆ F i by initializing F = ∅, adding points y ∈ F i to F until F and X are no longer separable, and then iterating over all points y ′ in F and checking whether their removal leads to a separable set. In the latter case, we keep y ′ in F ; otherwise, we remove it. Although this procedure is costly as it requires to solve Θ(|F i |) LPs to find F , preliminary experiments revealed that the running time of the cutting plane model can be reduced drastically when using the sparsified inequalities.
Since we expect the separation problem of (2c) to be difficult even for integer points, we only heuristically separate non-integral points s ⋆ ∈ [0, 1] Y ×[k] in our implementation. To this end, for each i ∈ [k], we again initialize an empty set F and iteratively add y ∈ Y in non-increasing order w.r.t. s ⋆ yi until F ∈ C and s ⋆ violates the inequality (or we know that such an inequality cannot be violated).

Algorithmic Aspects of the Column Generation Model
In contrast to the compact model (1), the number of variables in (3) grows exponentially in |Y |, which makes it already challenging to solve the LP relaxation of (3). In our implementation, we thus use a branch-and-price procedure for solving (3), i.e., we use a branch-and-bound procedure in which each LP relaxation is solved by column generation. In the following, we discuss the different components of the branch-and-price procedure.
Solving the Root Relaxation At the root node of the branch-and-bound tree, we are given a subset I ′ of all possible variables in I and solve the LP relaxation of (3) restricted to the variables in I ′ . To check whether the solution obtained for the variables in I ′ is indeed an optimal solution of the LP relaxation, we need to solve the pricing problem, i.e., to check whether all variables in I have nonnegative reduced costs. Since the pricing problem is equivalent to the separation problem for the dual, we determine the dual of the root node LP relaxation of (3), which is given by The pricing problem at the root node is thus to decide, for given dual weights α y , y ∈ Y , whether there exists a set I ∈ I with y∈I α y > 1. Unfortunately, we cannot expect to solve this problem efficiently in general.

Proposition 4.3.
Let X ⊆ Z d be finite and lattice-convex, let Y ⊆ Z d \ X be finite, and let α y ≥ 0 be a rational weight for y ∈ Y . Then, the pricing problem for the LP relaxation of (3), i.e., deciding whether there exists I ∈ I(X, Y ) with y∈Y α y > 1, is NP-hard.
Proof. Note that the pricing problem is equivalent to finding a set I ∈ I(X, Y ) that maximizes the value y∈I α y . If all weights α y , y ∈ Y , have the same value α > 0, the problem reduces to find a set I ∈ I of maximum cardinality. The latter problem is NP-hard even if X consists of a single point, in which case it reduces to the open hemisphere problem, see Johnson & Preparata [16].
To solve the pricing problem, we use a mixed-integer program that is a variant of (1) with k = 1. The only difference is that instead of minimizing the number of needed inequalities, we maximize the expression y∈Y α y s y1 . If this value is at most 1, we have found an optimal solution of the LP relaxation. Otherwise, we have found a variable z I with negative reduced cost, add I to I ′ , and iterate this procedure until all reduced costs are non-negative. In our implementation, we initialize the set I ′ by Branching Strategy Let u be a node of the branch-and-bound tree and denote by z u an optimal solution of the LP relaxation at node u. A classical branching strategy is to select a variable z I with z u I / ∈ Z and to create two child nodes u 0 and u 1 by enforcing z I = 0 in u 0 and z I = 1 in u 1 . While the branching decision z I = 1 has strong implications for computing rc ε (X, Y ) (we basically fix an inequality used in the relaxation), branching z I = 0 only rules out one of the exponentially many choices in I for a separated set.
To obtain a more balanced branching rule, we use the branching rule suggested by Ryan & Foster [23]. We are looking for two distinct variables z I and z J with z u I , z u J / ∈ Z such that both the intersection I ∩ J and symmetric difference I∆J of I and J are non-empty. Let y 1 ∈ I ∩ J and y 2 ∈ I∆J. Then, two child nodes u 0 and u 1 of u are created as follows. In u 0 , variables z I ′ are fixed to 0 if I ′ contains both y 1 and y 2 . In u 1 , we fix z I ′ to 0 if I ′ contains either y 1 or y 2 . That is, u 0 enforces y 1 and y 2 to be contained in different sets, and u 1 forces them to be contained in the same set I ′ . This branching rule obviously partitions the integer solutions feasible at node u. To show its validity it is thus sufficient to show that for every non-integral solution z u the sets I and J exist. Proof. Let I ∈ I ′ be such that z u I / ∈ Z. Then, z u I ∈ (0, 1), since z u is an optimal solution of the LP relaxation. Due to (3b), for every y ∈ I, there exists J y ∈ I ′ \ {I} with y ∈ J y such that z u J y > 0. For at least one J y we have z u J y ∈ (0, 1), because otherwise, we could improve the objective value of z u by setting z u I to 0 and still satisfying all constraints. Such a set J y together with I satisfy the properties in the statement of the lemma: Since y is contained in both I and J y , we have I ∩ J y = ∅. Moreover, as I = J y , I∆J y = ∅.
In our implementation, we compute for each variable z u I its fractionality θ(I) = 1 2 − min{z u I , 1 − z u I }. Then, we select I and J such that θ(I) + θ(J) is maximized; the branching candidates y 1 ∈ I ∩ J and y 2 ∈ I∆J are selected arbitrarily.

Solving LP Relaxations in the Tree
To not re-generate variables that have been fixed to 0 by the branching rule, we need to incorporate the branching decisions active at a node of the branch-and-bound tree into the pricing problem. This can easily be done by adding linear constraints to the root node formulation of the pricing problem. If a branching decision was that y 1 and y 2 shall be contained in different sets, we add s y11 + s y21 ≤ 1 to the pricing problem. The branching decision that y 1 and y 2 have to be contained in the same set can be enforced by the constraint s y11 = s y21 .

Numerical Experiments
The aim of this section is to compare the practical performance of the three models for computing rc ε (X, Y ) as well as their enhancements. To this end, we have implemented all three models in C/C++ using SCIP 7.0.3 as modeling and branch-and-bound framework and SoPlex 5.0.2 to solve all LP relaxations. All branching, propagation, separation, and pricing methods are implemented using the corresponding plug-in types of SCIP. Since we are not aware of an alternative separation routine for hiding set cuts, we compute all hiding sets of size two in a straightforward fashion before starting the branch-and-bound process. During the solving process, we separate these inequalities if the corresponding cuts are violated. To handle symmetries via lexicographic orderings, we use SCIP's internal plug-ins cons orbitope, cons orbisack, and cons symresack that implement the methods discussed in Section 4; the branching and pricing plug-ins for the column generation model strongly build up on the corresponding plug-ins of the binpacking example provided in the SCIP Optimization Suite. All convex hull computations have been carried out using cdd 0.94m [9] and graph symmetries are detected using bliss 0.73 [17].
Our implementation is available online at github 3 .

Implementation Details
All models admit some degrees of freedom that we detail in the following. Both the compact model and the cut model require an upper bound on the relaxation complexity. In both models, we impose the trivial upper bound which is given by the number of facets of conv(X). We also use the facet description to derive an initial solution for both models. In the column generation model, we need to select a subset of I to define initial variables. We use the sets I ∈ I that are defined by the facet defining inequalities of conv(X), i.e., the sets of points in Y that are separated from X by the facet defining inequalities. Moreover, we include the singleton sets {y}, for y ∈ Y , to make sure that the LP relaxation remains feasible after branching.
Settings To encode the different settings that we have tested, we make use of the following abbreviations: hiding Whether hiding set cuts are added (1) or not (0).
sym. Which symmetry method is used: none (0), simple (s), or advanced (a), where simple is (4) and (5), and advanced uses (4) and additionally enforcing lexicographically maximal solutions based on symmetries of X and Y .
Note that we do not report on results for the intersection propagation algorithm. This is because, in preliminary experiments, we have seen that its running time is very high, in particular, because it needs to compute in each iteration O(|Y |) convex hulls. As a result, we could hardly solve any instance, not even small ones.

Test Sets
In our experiments, we have used three different test sets: basic The sets X are the vertices of the 0/1 cube, the crosspolytope, or the standard simplex in dimensions d ∈ {3, 4, 5}. For X ⊆ Z d , the sets Y consist of all points in Z d \ X whose ℓ 1 -distance to X is at most k, where 1 ≤ k ≤ 10 − d. The reason for smaller distance in higher dimension is that the problems get considerably more difficult to solve with increasing k.
downcld This test set consists of 99 full-dimensional subsets X of {0, 1} 5 that correspond to down-closed subsets (or abstract simplicial complexes) of the Boolean lattice on 5 elements. The corresponding sets Y are the points in Z 5 \ X whose ℓ 1 -distance to X is at most k ∈ {1, 2, 3}.
The sets X have been generated by the natural one-to-one correspondence between inclusion-maximal sets in a down-closed family and antichains in the Boolean lattice.
sboxes The test set comprises 18 instances modeling 4-bit (12 instances) and 5-bit (6 instances) S-boxes, which are certain non-sparse Boolean functions arising in symmetric-key cryptography. The derived sets X are contained in {0, 1} 8 and {0, 1} 10 , respectively, and Y are the complementary binary points. These instances have also been used by Udovenko [27] who solved the full model (3), i.e., without column generation.
The basic instances feature various aspects that might be relevant for computing rc(X) via computing a series of values rc ε (X, Y ) for different Y and ε according to [3]: The cube is parity complete, thus there exists a small set Y such that rc(X) = rc ε (X, Y ) (in fact, this set is {−1, 0, 1, 2} d \ X); the crosspolytope has an interior integer point and thus there exists a (potentially large) finite set Y with rc(X) = rc ε (X, Y ); for the simplex ∆ 4 in R 4 , no finite set Y exists with rc(∆ 4 ) = rc(∆ 4 , Y ); see [2]. That is, Since the standard simplex ∆ d is a down-closed subset of {0, 1} d , the smallsized downcld instances might be good candidates for further examples X such that rc ε (X, Y ) < rc Q (X), for every finite set Y ⊆ Z d and for ε > 0 small enough. Our aim for selecting these instances is thus to identify whether there are potentially further candidates for sets X whose relaxation complexity cannot be computed via finite sets Y .
Finally, the sboxes instances are used to investigate whether our techniques are suited to compute rc ε (X, Y ) also in higher dimensions. This is relevant, among others, in the field of social choice or symmetric cryptanalysis, where the aim is to find rc(X, {0, 1} d ) for sets X ⊆ {0, 1} d .

Computational Setup
All experiments have been run on a Linux cluster with Intel Xeon E5 3.5 GHz quad core processors and 32 GB memory. The code was executed using a single thread and the time limit for all computations was 4 h per instance.
All mean numbers are reported in shifted geometric mean n − s to reduce the impact of outliers. For mean running times, a shift of s = 10 is used; for nodes of the branch-and-bound tree, we use s = 100. The value of ε in computing rc ε (·, ·) is set to 0.001. The upper bound on the number of inequalities needed in the compact and cutting plane model is given by the number of facets of conv(X). We also provide an initial primal solution corresponding to a facet description of conv(X).

Results for Test Set basic
Due to our choice of the sets X and Y , the basic test set comprises 18 cube, crosspolytope, and simplex instances, respectively. Table 1 shows the results for the compact model. For the plain compact model, we observe that SCIP can already solve quite some instances, but, in comparison to the enhanced variants, the running times are rather high. Checking each of the enhancements separately, handling symmetries is most important to reduce running time and to increase the number of instances solvable within the time limit. Interestingly, handling symmetries on the a-variables modeling the inequalities in a relaxation performs better than handling the symmetries of the points to separate. Adding hiding set cuts to the problem formulation is also beneficial, whereas the convexity propagator seems to harm the solving process in particular for cube instances. The worse performance for enabled propagation cannot be explained on the running time of the propagator: For cube instances, e.g., the maximum running time per instance of the propagator was 27 s, which is much smaller than the increase of mean running time. Thus, it seems that the found reductions guide the branch-and-bound search into the wrong direction or make it more difficult for SCIP to find other reductions. Combining simple symmetry handling and hiding set cuts leads consistently to the best results, reducing mean running time for cube instances by 87 %, for crosspolytope instances by 89 %, and simplex instances by 74 %. In particular, the combined setting can solve all cube instances and almost all crosspolytope and simplex instances within the time limit.
Next, we discuss the column generation model for which we only compare two variants: we either disable or enable hiding set cuts in the pricing problem. Since the convexity propagator does not seem to be helpful for the compact model, we do not enable it when solving the pricing problem. Moreover, symmetry handling is not important, because there is only one inequality to be identified by the pricing model.
Comparing the column generation model with disabled hiding set cuts, we can see that it performs for cube and crosspolytope instances much better than the plain compact model: the running time for cubes reduces by 82 % and for cross polytopes by 74 %. For cubes, all solvable instances are solved within the root node which is, on the one hand, because of the strong dual bound as described in Proposition 3.2. On the other hand, the generated sets I ∈ I allow heuristics to find high quality solutions yielding a matching primal bound. For  crosspolytopes, all instances of 3-dimensional sets X can be solved within the root node; for 4-and 5-dimensional sets, however, SCIP needs to start branching to find an optimal solution. Looking onto results on a per-instance basis reveals that the pricing problems become considerably harder if d and k increases. For example, SCIP is only able to process 2 nodes of the branch-and-bound tree for d = k = 5. For the simplex instances, the column generation model needs approximately twice as much time as the plain compact model, which is again explained by the very high running time of the pricing problem.
Enabling also hiding set cuts helps to solve the pricing problems more efficiently. In comparison with the enhanced compact model, however, the enhanced column generation model is only competitive on the cube instances. On the crosspolytope and simplex instances, it is much slower.
Finally, we consider the cutting plane model. In the plain version, this model can hardly solve any instance efficiently. Comparing the different enhancements with each other, we can see, analogously to the compact model, that adding hiding set cuts and handling symmetries is beneficial. Interestingly, the convexity propagator helps to improve the running time if both the previous enhancements are enabled by 90 %-95 %, leading to the best setting for this model. But even this winner setting cannot compete with the enhanced compact model.
From the results using the compact and cutting plane model, we draw the following conclusion regarding the convexity propagator. In principle, this method models the important aspect that the points being cut by an inequality form a lattice-convex set. The cutting plane method can thus benefit from the propagator as this property is not encoded in the model. The compact model, however, makes use of additional variables modeling the inequalities of a relaxation. Since the convexity propagator does not improve SCIP's performance, we conclude that these additional variables already sufficiently encode the lattice-convexity of cut points.
In summary, the column generation model provides very good primal and dual bounds. If these bounds match, rc ε (X, Y ) can be computed rather effi-ciently if not too many pricing problems need to be solved. However, if the bounds do not match, the NP-hardness of the pricing problem strikes back and solving many further pricing problems is too expensive. In this case, the compact model is a rather effective alternative that also allows to compute rc ε (X, Y ) for d = 5 in many cases.

Results for Test Set downcld
In this section, we turn the focus on 5-dimensional 0/1 down-closed sets. On the one hand, our aim is to investigate whether the findings of the previous section carry over to a much broader test set in dimension 5. On the other hand, we are interested in identifying further sets X ⊆ {0, 1} 5 with rc ε (X, Y ) < rc Q (X) for every choice of a finite set Y ⊆ Z d and ε > 0 small enough. Because of our results on the basic test set, we did not run any experiments using the cutting plane model as we expect that it can hardly solve any instance. Instead, we consider a hybrid version of the compact model and the column generation model: We only solve the column generation model's LP relaxation to derive a strong lower bound on rc ε (X, Y ) and to find good primal solutions. Both are transferred to the compact model with the hope to drastically reduce solving time. The running times and number of nodes reported for the hybrid model are means of the total running time and total number of nodes for solving the LP relaxation in the column generation model and the resulting compact model. Table 4 shows aggregated results for the 99 instances of the downcld test set for different ℓ 1 -neighborhoods Y of X (radius 1-3). While the plain compact model is able to solve two third of all instances for radius 1, computing rc ε (X, Y ) for larger radii becomes much harder. As the plain model can hardly solve any instance for radius at least 2, there is definitively a need for model enhancements. In general, the same observations as in the previous section can be made: symmetry handling and adding hiding set cuts improve the solution process a lot. The biggest impact is achieved by symmetry handling; the convexity propagator is not helpful in the best setting. However, sometimes it can improve the running time, e.g., if the "wrong" symmetry handling method is used.
For radius 2 and 3, we find that the simple symmetry handling methods perform much better than the advanced methods. Using hiding set cuts and simple symmetry handling is 59 % faster than the corresponding setting with advanced symmetry handling if the radius is 2; for radius 3, it is 45 % faster. Moreover, simple symmetry handling can solve all 99 instances for radius 2 (resp. 69 instances for radius 3), whereas the advanced setting can only solve 75 (resp. 14) instances. Interestingly, for radius 1, the advanced setting is 26 % faster than the simple setting. A possible explanation is based on the nature of the advanced setting: Each inequality defining a relaxation of X w.r.t. Y defines a pattern on the points from Y that are cut by this inequality. The advanced method enforces that the cut patterns of the inequalities are sorted lexicographically based on a sorting of the elements of Y . Since the results of the lexicographic comparison is determined by the first position in which two vectors differ, it is unlikely that points having a late position in the ordering of Y are very relevant for the lexicographic constraint. Thus, the symmetries are in a certain sense mostly handled for the early points in this ordering. In contrast to this, the simple method takes the geometry of the inequalities in a relaxation into account by sorting inequalities based on their first coefficients. Together with other components of the solver, this seems to have more implications on the cut points from Y if the radius becomes larger. In comparison to the enhanced compact model, the column generation model is again inferior. For radius at least 2, it can hardly solve any instance and, as already discussed in the previous section, the reason for this is the long running time of the pricing models that need to be solved often at each node of the tree. This is reflected by the number of processed nodes during the branchand-price procedure that drops drastically (as the number of solved instances) if the radius is getting larger. However, we can again observe that the root node can be solved relatively efficiently and that the obtained primal and dual bounds are rather strong. This is reflected in the hybrid model, which solves most instances and reduces the running time (in comparison to the best compact model) by 52-56 % for radius 2 and 3. For radius 1, the running times are comparable.
Regarding the usefulness of hiding set cuts in the hybrid model, we observe that they are essential for solving the downcld instances efficiently. They allow to solve all instances for radius 1 and 2 and improve on the hybrid setting without cuts by 74 % and 88 %, respectively. This effect is even more dominant for radius 3, where it significantly increases the number of solvable instances, reducing the running time by 72 %. It is also noteworthy that the hybrid setting with hiding set cuts is the only setting allowing to solve 80 instances, which improves the running time of the compact model by 51 %. In summary, based on our experiments, the hybrid model is the best choice for computing rc ε (X, Y ) as it combines the strong bounds from the column generation model with the ability of the compact model to quickly solve LP relaxations within the branch-andbound tree. In particular, it benefits from hiding set cuts since their implications are very difficult to be found by SCIP.
Finally, concerning our goal to identify candidates for sets X ⊆ {0, 1} 5 such that rc ε (X, Y ) < rc Q (X) for all finite Y ⊆ Z 5 and ε > 0 small enough, our experiments for radius 3 revealed the following: If Y (X) are the integer points in the ℓ 1 -neighborhood of X with radius 3, then there are three sets X such that rc ε (X, Y (X)) = 4. These sets are ∆ 5 , ∆ 5 ∪ {e 1 + e 2 } and ∆ 4 × {0, 1}. Moreover, there are 16 sets X with rc ε (X, Y (X)) = 5. It is left open for future research to identify which other sets than ∆ 5 satisfy rc ε (X, Y (X)) < rc Q (X). Note that rc Q (X) ≥ 6, whenever X ⊆ {0, 1} 5 is full-dimensional, because rational relaxations must be bounded.

Results for Test Set sboxes
The results for the sboxes test set are summarized in Table 5. Note that we do not report on results for the 10-dimensional instances in the compact model with enabled hiding set cuts, because all these experiments hit a memory limit of 20 GB. The reason is that these models grow very large even without any enhancements as we use the number of facets of conv(X) to upper bound rc ε (X, Y ); the number of facets for these instances ranges between 888 and 2395. For the largest instances, even the basic compact model hits the memory limit. Adding hiding set cuts for the remaining instances causes that all instances hit the memory limit. But also for the smaller instances, SCIP is   hardly able to solve any of these instances even if model enhancements are enabled due to huge number of variables and constraints. In contrast to this, we see that the column generation model performs extremely well for the problems in dimension 8. It can solve all twelve 8-dimensional instances within the time limit, on average in 185.7 s if hiding set cuts are disabled and in roughly twice this amount of time with enabled hiding set cuts. An explanation for the worse behavior with enabled cuts is that the number of hiding set cuts increases drastically in comparison with lower dimensional problems. Thus, creating and separating these cuts is a non-trivial task. For dimension 10, the column generation model is also able to solve 2 out of 6 instances within the time limit.
Finally, the hybrid model performs worse than the column generation model. Although the derived bounds from solving the column generation model's LP relaxation yield again very good bounds on the relaxation complexity, the value of rc(X, {0, 1} d ) can still be large if d ∈ {8, 10}. Thus, also the compact model embedded in the hybrid model is struggling with the number of variables and constraints. For this reason, computing rc(X, {0, 1} d ) via the column generation model is most competitive.

Conclusions
Being able to compute the exact value of the quantity rc ε (X, Y ) is highly relevant in many areas, such as, social choice, symmetric cryptanalysis, or machine learning. For this reason, we have proposed three different models that allow to compute rc ε (X, Y ) using mixed-integer programming techniques. As our experiments reveal, each of these models comes with advantages and disadvantages. The compact model, for example, works well in small dimensions as the number of variables and inequalities is small and it encapsulates all essential information about rc ε (X, Y ). In higher dimensions, however, the dual bounds of the compact model become weaker. In this case, the column generation model provides very good bounds that can be transferred to the compact model to still compute rc ε (X, Y ) rather efficiently if d = 5. But if the dimension d grows even larger, only the column generation model seems to be competitive as it does not scale as badly as the compact model when rc ε (X, Y ) increases. The main reason is that the compact model is relying on a good upper bound on rc ε (X, Y ) to be indeed compact.
These findings thus open the following directions for future research. Since the compact model requires a good upper bound on rc ε (X, Y ), it is natural to investigate heuristic approaches for finding ε-relaxations of X or to develop approximation algorithms. Moreover, since the column generation model becomes more relevant if d is large, it is essential that the pricing problem can be solved efficiently. Since the pricing problem is NP-hard, also here a possible future direction could be to develop heuristics or approximation algorithms for solving it. For both the compact and column generation model, hiding set cuts turned out to be useful. However, we are not aware of an efficient routine for generating these cutting planes. Thus, it is natural to devise an efficient scheme for generating hiding set cuts on the fly. Finally, as additional inequalities such as hiding set cuts and symmetry handling inequalities drastically improved the performance of the compact model, the development of further inequalities modeling structural properties of relaxation complexity might allow to solve the compact model even more efficiently.