IJOC-0001-2024.00

Barbato, Ceselli and Messana Mathematical programming algorithms for convex
hull approximation with a hyperplane budget

Mathematical Programming Algorithms for Convex Hull Approximation with a
Hyperplane Budget

Michele Barbato, Alberto Ceselli, Rosario Messana Department of Computer
Science, Università degli Studi di Milano, Via Giovanni Celoria 18, 20133
Milan, Italy
Contact<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
We consider the following problem in computational geometry: given, in the
$d$-dimensional real space, a set of points marked as positive and a set of
points marked as negative, such that the convex hull of the positive set does
not intersect the negative set, find $\mathcal{K}$ hyperplanes that separate,
if possible, all the positive points from the negative ones. That is, we
search for a convex polyhedron with at most $\mathcal{K}$ faces, containing
all the positive points and no negative point. The problem is known in the
literature for pure convex polyhedral approximation; our interest, stems from
its possible applications in constraint learning, where points are feasible or
infeasible solutions of a Mixed Integer Program, and the $\mathcal{K}$
hyperplanes are linear constraints to be found.

We cast the problem as an optimization one, minimizing the number of negative
points inside the convex polyhedron, whenever exact separation cannot be
achieved. We introduce models inspired by support vector machines and we
design two mathematical programming formulations with binary variables. We
exploit Dantzig-Wolfe decomposition to obtain extended formulations, and we
devise column generation algorithms with ad-hoc pricing routines.

We compare computing time and separation error values obtained by all our
approaches on synthetic datasets, with number of points from hundreds up to a
few thousands, showing our approaches to perform better than existing ones
from the literature. Furthermore, we observe that key computational
differences arise, depending on whether the budget $\mathcal{K}$ is sufficient
to completely separate the positive points from the negative ones or not. On
$8$-dimensional instances (and over), existing convex hull algorithms become
computational inapplicable, while our algorithms allow to identify good convex
hull approximations in minutes of computation. They can also be used for
parametric analyses, being the hyperplane budget a hyperparameter trading
approximation quality for model complexity.

Convex Hull, Mathematical Programming, Column Generation

##  1 Introduction

We consider the following geometric problem. Two finite sets of points in
$\mathbb{R}^{d}$ are given, labelled as positive and negative respectively. We
assume that the convex hull of the positive points does not contain any
negative point. Our aim is to find a convex polyhedron having at most
$\mathcal{K}$ faces, which contains all the positive points and, if possible,
no negative points.

Besides its applications in fields like image processing (Jayaram and Fleyeh
2016), and the interest into the combinatorial problem on its own, we are
motivated by its perspective use in constraint learning (Lombardi and Milano
2018). Let us assume that a constrained optimization problem is given,
together with a suitable sampling of solutions (encoded as points in
$\mathbb{R}^{d}$), some of them being feasible (the positive ones) and others
being infeasible (the negative ones). Let us also assume to know that the
region of feasible solutions can be described as a polyhedron, defined by at
most $\mathcal{K}$ linear constraints. The problem of finding them is
equivalent to our problem described above: finding a convex polyhedron having
at most $\mathcal{K}$ faces, enclosing all the points that encode feasible
solutions, and excluding all those encoding infeasible ones. Hardness results
date back to Megiddo (1988).

Without the requirement of limiting the number of faces to $\mathcal{K}$, the
problem could be solved by producing the convex hull of the feasible points
using well established techniques (Bradford Barber, Dobkin, and Huhdanpaa
1996), and in particular by finding the equations describing the facets of the
convex hull. The set of negative points becomes irrelevant, but the
combinatorial complexity of these approaches makes them computationally viable
only for positive sets of small size. Faster alternatives exist for convex
hull approximation (Balestriero, Wang, and Baraniuk 2022), adjusting the
solution in post processing to contain all the positive points, when
necessary.

Since we assume that the budget $\mathcal{K}$ is given as input, it is
possible that $\mathcal{K}$ hyperplanes are not sufficient to separate all the
positive points from the negative ones. Our goal, in this case, is to produce
a _small_ approximation of the convex hull. We therefore map the problem to an
optimization one, in which the number of negative points falling inside the
convex polyhedron is minimized (and no positive points falls outside the
polyhedron). The intuition behind this choice is the following: if the points
in the positive and negative sets came from a uniform random drawing (and were
dense enough), then the $\mathcal{K}$ hyperplanes define a convex polyhedron
that over-approximates the convex hull of the feasible points, and whose
volume is minimized.

(a) Dataset

(b) Convex hull

(c) Hull approximation

Figure 1: Given a dataset of positive and negative points (Fig. 1(a)), one can
alternatively compute the convex hull of the positive points (Fig. 1(b)) or
fix a budget on the number of hyperplanes and compute an over-approximation of
the convex hull that satisfies the budget requirement and does not contain any
of the negative points (Fig. 1(c)).

#### Contributions.

We formally introduce the problem and we discuss similarities and differences
with related applications from the literature (Section 2). Then, we
investigate the use of mixed-integer linear programming (MILP) to design new
resolution algorithms. In particular, we propose two mixed-integer linear
programming models (Section 3). We derive extended formulations from them,
through the use of Dantzig-Wolfe decomposition. We design column generation
algorithms, employing techniques such as early stopping, and pricing solutions
polishing. Additionally, we design a heuristic parallel algorithm for the
pricing problem that proves to be empirically advantageous. We also introduce
fast upper bounding and master initialization greedy methods (Section 4). We
compare with both state-of-the-art convex hull algorithms, and to alternative
polyhedral separation techniques from the literature on a wide range of
instances (Section 5), focusing on how dimensionality, number of points and
hyperplane budget impact the behaviour of each method. Some conclusions are
finally drawn (Section 6), highlighting how our methods improve on existing
ones.

##  2 Notation and Background

We call our problem _polyhedral convex hull approximation with hyperplane
budget_ (PCAB). It can be formalized as follows. We are given two indexed sets
$\mathcal{D}^{+}=\\{\bm{a}^{+}_{i}\\}_{i=1}^{m}$ and
$\mathcal{D}^{-}=\\{\bm{a}^{-}_{i}\\}_{i=1}^{n}$ of points in $\mathbb{R}^{d}$
and a positive integer $\mathcal{K}$, the _hyperplane budget_.

A hyperplane $h\in\mathbb{R}^{d}$ is called _valid_ if it admits a valid
equation, i.e. it exist $b\in\mathbb{R}$ and $\bm{w}\in\mathbb{R}^{d}$ such
that $h=\\{\bm{x}\in\mathbb{R}^{d}\colon b+\bm{w}\cdot\bm{x}=0\\}$ and for
each point $a_{i}\in\mathcal{D}^{+}$ it holds that
$b+\sum_{j=1}^{d}w_{j}a_{ij}\geq 0$. For the sake of brevity, we identify an
equation with its coefficients $(b,\bm{w})$.

A valid solution to PCAB is a set of $\mathcal{K}$ valid equations
$s=\left\\{(b^{k},\bm{w}^{k})\right\\}_{k=1}^{\mathcal{K}}$. The _separation
error_ associated with $s$ is:

$err(s)=\left|\left\\{\bm{a}^{-}_{i}\ |\ i\in\\{1,\ldots,n\\}\wedge
b^{k}+\sum_{j=1}^{d}w_{j}^{k}a_{ij}\geq 0\ \forall
k\in\\{1,\ldots,\mathcal{K}\\}\right\\}\right|$.

A valid solution $s^{*}$ is optimal if $err(s^{*})$ is minimum over all valid
solutions. PCAB consists in the problem of finding one of these optimal
solutions.

Terminology.  We define the following index sets: $I^{+}:=\\{1,\ldots,m\\}$,
$I^{-}:=\\{1,\ldots,n\\}$, $I:=\\{1,\ldots,m+n\\}$, $J:=\\{1,\ldots,d\\}$,
$K:=\\{1,\ldots,\mathcal{K}\\}$.

Let $\ell:\mathcal{D}^{+}\cup\mathcal{D}^{-}\rightarrow\\{-1,1\\}$ be the
labeling function such that $\ell(\mathcal{D}^{+})=1$ and
$\ell(\mathcal{D}^{-})=-1$. We call _positive_ the points labeled 1 and
_negative_ the remaining ones. When we need a uniform naming convention for
positive and negative points we use that $\bm{a}_{i}:=\bm{a}^{+}_{i}$ for
every $i\in I^{+}$ and $\bm{a}_{i}:=\bm{a}^{-}_{i-m}$ for every $i\in
I\setminus I^{+}$.

We say that a hyperplane separates the positive points from a certain subset
of negative points, indexed by $\bar{I}^{-}\subseteq I^{-}$, if it admits a
valid equation $(b,\bm{w})$ such that $b+\sum_{j=1}^{d}w_{j}a_{ij}<0$ for
every $i\in\bar{I}^{-}$. The _indicator vector_ of $\bar{I}^{-}$ is the vector
$\bm{x}\in\\{0,1\\}^{n}$ such that for each $i\in I^{-}$, $x_{i}=1$ if and
only if $i\in\bar{I}^{-}$.

A solution $s$ of PCAB separates the positive points from all the negative
points if and only if $err(s)=0$. The minimum hyperplane budget $\mathcal{K}$
such that PCAB admits a solution with separation error 0 is called _minimum
separation budget_.

The concatenation of two real-valued vectors $\bm{x}$ and $\bm{y}$ is
indicated as $\bm{x}\mathbin{\|}\bm{y}$. Their element-wise product is written
$\bm{x}\circ\bm{y}$.

Related Work.  Despite its simplicity, we are not aware of algorithmic
approaches specifically targeting our PCAB . At the same time, it links to
several relevant problems from the literature, not only from computational
geometry, but also from machine learning and optimization.

Concerning computational geometry, our PCAB is related to the problem of
computing the convex hull of a set of points in terms of the equations of the
delimiting hyperplanes. It is well know in the literature: dedicated software
packages were also conceived, engineered and maintained by the community
(Assarf, Gawrilow, Herr, Joswig, Lorenz, Paffenholz, and Rehn 2017). Computing
the convex hull in arbitrary dimension can be done exactly applying the
algorithm QuickHull (Bradford Barber, Dobkin, and Huhdanpaa 1996, Bykat 1978,
Eddy 1977) or approximately using, for example, more recent algorithms
(Balestriero, Wang, and Baraniuk 2022). Solving PCAB gives an over-
approximation of the convex hull of the positive points but, on the other
hand, PCAB focuses on finding much simpler separation boundaries by setting a
hyperplane budget. It is worth noting, for instance, that while convex hull
algorithms always produce polytopes, our PCAB models also allow unbounded
convex polyhedra as valid solutions.

Another related problem is learning intersections of halfspaces using points
labeled as positive or negative. As in the PCAB , the focus in on identifying
the structure of a specific unknown polyhedron containing the positive points
but none of the negative ones. However, unlike the PCAB , an implicit
representation of the polyhedron is normally given, for example, by means of
an oracle that augments the learning process. In this setting, polynomial-time
algorithms w.r.t the dimension and the number of faces of the polyhedron were
given by Blum et al. (1995), Baum (1990) and Kwek and Pitt (1996). Learning
from random examples without access to an oracle, as in our case, is proved to
be a challenging task whose related complexity results can be independent of
the distribution from which the data is drawn (Klivans, O’Donnell, and
Servedio 2004) or distribution-dependent, as in previous works by Baum (1991)
and Vempala (1997).

A close connection also exists with convex polyhedral classification (which is
sometimes named convex polyhedral separation in the literature), a supervised
learning task that consists in finding a binary classifier whose positive
region is a convex polyhedron. It is a particular case of binary
classification with piece-wise linear separation boundary. A foundational work
is (Astorino and Gaudioso 2002), in which a weighted sum of the classification
error on the positive class and on the negative one is minimized solving, at
each iteration, a linear program providing a descent direction. The problem
has received a fair amount of attention over the years, with solution
approaches that cover mixed-integer programming (Orsenigo and Vercellis 2007),
statistical learning (Manwani and Sastry 2010) local SVM-like classification
(Leng, He, Liu, ping Qin, and Li 2020), and many others (Astorino, Avolio, and
Fuduli 2023, 2021, Astorino and Fuduli 2015, Gottlieb, Kaufman, Kontorovich,
and Nivasch 2018, Kantchelian, Tschantz, Huang, Bartlett, Joseph, and Tygar
2014, Raviv, Hazan, and Osadchy 2018, Strekalovsky, Gruzdeva, and Orlov 2015).
However, the focus of convex polyhedral classification is on finding a
classifier with good generalization properties. Unlike our PCAB , the positive
and negative sets are known only partially at training time. Therefore, they
typically aim at separation boundaries with margins as large as possible,
while PCAB tries to find convex hull approximations which are as tight as
possible. Still some works in this area can be adapted to our case, as
discussed in the remainder.

We finally mention the link to linear constraint learning. Some algorithms
have been proposed that provide explicit representations of linear constraints
based on examples of feasible and infeasible solutions, e.g. (Pawlak and
Krawiec 2017, Schede, Kolb, and Teso 2019). These works represent just a
sketch of the broader area of constraint learning. The interested reader can
refer to (Fajemisin, Maragno, and den Hertog 2024) for a general survey and
the description of a conceptual framework for development of constraint
learning techniques. Indeed, our PCAB algorithms are designed with a
perspective application in a constraint learning framework. On one side as
observed in (Messana, Chen, and Lodi 2024), pure optimization algorithms as
ours are hardly effective when used _alone_ for constraint separation, because
there is no guarantee that available feasible and infeasible solutions are
uniformly distributed. On the other side, when optimization is coupled with
proper _solutions sampling_ , approaches like PCAB become promising.

##  3 Models

We introduce two compact mixed-integer programming models for PCAB . The first
one is directly inspired by Support Vector Machines (SVMs), borrowing from
earlier works on polyhedral separation. The second one incorporates key ideas
to improve the approximation performance.

A classical mathematical programming model for the training of SVMs is
reported in Appendix 7. It uses one continuous variable $w_{j}$ for each
dimension, encoding coefficients of a hyperplane, and one continuous (non-
negative) variable $\xi_{i}$ to track the absolute misclassification error on
point $i\in I$.

In fact, the SVM model can be adapted to solve PCAB as follows. We name this
formulation as Model A, which is reported in full in Appendix 7. Model A
prescribes to find $\mathcal{K}$ hyperplane equations, encoded by variables
$w_{j}^{k}$ with $j\in J$ and $k\in K$. A set of non-negative variables
$\xi_{i}^{k}$ encode the distance from point $i\in I^{-}$ to hyperplane $k$,
if $i$ is on the wrong side of the hyperplane $k$; they take value $0$
otherwise. Since PCAB requires every hyperplane in a solution to be valid, we
neglect the violation variables corresponding to positive points, imposing
them to be always on the correct side of each hyperplane.

We introduce binary variables $x_{i}^{k}$ for $i\in I^{-}$ and $k\in K$,
having the following meaning: when $x_{i}^{k}=1$, it means that point $i$ is
_separated_ by the $k^{th}$ hyperplane. That is, we impose that $x_{i}^{k}$
can be set to $1$ only if the point $i\in I^{-}$ is on the negative side of
hyperplane $k$. Each negative point must be separated by one hyperplane. We
still allow negative points to be on the positive side of all hyperplanes
(thereby being misplaced), at the price of a violation penalty. The objective
function sums up the weighted violations over the negative points. For each
negative point, the violation cost is considered only with respect to a single
hyperplane (that yielding minimum violation). The weights are expressed as
values of a function $\hat{c}:\mathbb{R}^{d}\to[0,+\infty)$.

The potential and limits of using binary variables in SVM-like models are
discussed in the literature (Carrizosa, Martin-Barragan, and Romero Morales
2010). Additionally, as we show in Section 5.2, the solutions to PCAB obtained
solving Model A using Gurobi are weak in terms of approximation of the convex
hull. It is not surprising, since the objective function of Model A targets a
margin maximization, as in SVMs. Furthermore, it is quadratic, involving the
product of variables $x_{i}^{k}$ and $\xi_{i}^{k}$. As final formulation step,
therefore, we introduce the following mixed-integer programming model, denoted
Model B.

min $\displaystyle\sum_{i\in I^{-}}c(\bm{a}^{-}_{i})e_{i}$ (1) s.t.
$\displaystyle\sum_{k\in K}x_{i}^{k}=1-e_{i}$ $\displaystyle\forall i\in
I^{-}$ (2) $\displaystyle b^{k}+\sum_{j\in J}w_{j}^{k}a_{ij}\geq 1$
$\displaystyle\forall i\in I^{+},\forall k\in K$ (3) $\displaystyle
b^{k}+\sum_{j\in J}w_{j}^{k}a_{ij}\leq-1+\xi_{i}^{k}$ $\displaystyle\forall
i\in I^{-},\forall k\in K$ (4) $\displaystyle\xi_{i}^{k}\leq
M\cdot(1-x_{i}^{k})$ $\displaystyle\forall i\in I^{-},\forall k\in K$ (5)
$\displaystyle e_{i}\in\\{0,1\\}$ $\displaystyle\forall i\in I^{-}$ (6)
$\displaystyle\xi_{i}^{k}\geq 0$ $\displaystyle\forall i\in I^{-},\forall k\in
K$ (7) $\displaystyle x_{i}^{k}\in\\{0,1\\}$ $\displaystyle\forall i\in
I^{-},\forall k\in K$ (8) $\displaystyle b^{k}\in\mathbb{R}$
$\displaystyle\forall k\in K$ (9) $\displaystyle w_{j}^{k}\in\mathbb{R}$
$\displaystyle\forall j\in J,\forall k\in K$ (10)

Each binary variable $e_{i}$, with $i\in I^{-}$, represents the condition of
negative point $i$ to be incorrectly included in the convex hull. When
$e_{i}=0$, constraints (2) impose that exactly one $x_{i}^{k}=1$; in turn,
when $x_{i}^{k}=1$, constraints (5) impose that $\xi_{i}^{k}=0$, and
constraints (4) impose $i$ to be on the correct side of hyperplane $k$.
Constraints (3) impose all positive points $i\in I^{+}$ to be on the correct
side. The objective function (1) counts the number of misplaced negative
points (possibly weighting each of them differently). As technical remarks,
the weights in the objective function do not need to be non-negative (while
they do in Model A). Therefore, we can replace the weight function $\hat{c}$
with a real valued function $c:\mathbb{R}^{d}\to\mathbb{R}$. Furthermore, to
enhance the effect of counting the misplaced points, saving one
hyperparameter, no regularization term on $w$ was placed in the objective
function. Accordingly, the right hand sides of constraints (3) and (4) contain
constants ($1$ and $-1$) to avoid trivial solutions with $w=0$. Finally, the
big-M terms in constraints (5) have been introduced to avoid quadratic terms
in the objective function.

In order to improve the bounds that can be obtained from Model B we employed
Dantzig-Wolfe decomposition. Let

$\mathcal{F}^{k}=\left\\{\left(x_{i}^{k}\right)_{i\in
I^{-}}\mathbin{\|}\left(w^{k}_{j}\right)_{j\in
J}:\eqref{st:model_b_positive},\eqref{st:model_b_negative},\eqref{st:model_b_big_M},\eqref{st:model_b_domain_xi},\eqref{st:model_b_domain_x},\eqref{st:model_b_domain_b},\eqref{st:model_b_domain_w}\right\\}$

and let $\mathcal{G}^{k}=\left\\{\left(x_{i}^{k}\right)_{i\in I^{-}}\ |\
\exists\left(w^{k}_{j}\right)_{j\in J}\ s.t.\ \left(x_{i}^{k}\right)_{i\in
I^{-}}\mathbin{\|}\left(w^{k}_{j}\right)_{j\in J}\in\mathcal{F}^{k}\right\\}$
for each $k\in K$. Since $\mathcal{G}^{k}$ is a set of integer points, as
follows:

min $\displaystyle\sum_{i\in I^{-}}c(\bm{a}^{-}_{i})e_{i}$ s.t.
$\displaystyle\sum_{k\in K}x_{i}^{k}=1-e_{i}$ $\displaystyle\forall i\in
I^{-}$ $\displaystyle\left(x_{i}^{k}\right)_{i\in I^{-}}\in
conv\left\\{\mathcal{G}^{k}\right\\}$ $\displaystyle\forall k\in K$
$\displaystyle e_{i}\in\\{0,1\\}$ $\displaystyle\forall i\in I^{-}$

We relax the integrality constraints on the error variables obtaining $0\leq
e_{i}\leq 1\ \forall i\in I^{-}$.

Let $\Omega^{k}$ be the set of extreme (integer) points of
conv$\\{\mathcal{G}^{k}\\}$, and we denote (the coordinates of) each of them
as $\bar{x}_{i}^{kp}$ for every $p\in\Omega^{k}$. Valid values for $x_{i}^{k}$
can be obtained as a convex combination of these extreme points. Let us
introduce multipliers $z^{kp}$ representing such a convex combination, i.e.

$x_{i}^{k}=\sum_{p\in\Omega^{k}}\bar{x}_{i}^{kp}z^{kp}\ \ \ \forall i\in
I,\forall k\in K$ (11)

with $\sum_{p\in\Omega^{k}}z^{kp}=1$ for each $k\in K$. Constraints (11) can
be used to project out the $x_{i}^{k}$ variables:

min $\displaystyle\sum_{i\in I^{-}}c(\bm{a}^{-}_{i})e_{i}$ s.t.
$\displaystyle\sum_{k\in
K}\sum_{p\in\Omega^{k}}\bar{x}_{i}^{kp}z^{kp}=1-e_{i}$ $\displaystyle\forall
i\in I^{-}$ $\displaystyle\sum_{p\in\Omega^{k}}z^{kp}=1$ $\displaystyle\forall
k\in K$ $\displaystyle 0\leq z^{kp}\leq 1$ $\displaystyle\forall k\in
K,\forall p\in\Omega^{k}$ $\displaystyle 0\leq e_{i}\leq 1$
$\displaystyle\forall i\in I^{-}$

By observing the sets $\mathcal{G}^{k}$, one can notice all of them to contain
the same set of points, call it $\mathcal{G}$ (and for consistency, denote as
$\Omega$ the set of extreme points of the corresponding convex hull). Our
model can be further simplified as follows:

min $\displaystyle\sum_{i\in I^{-}}c(\bm{a}^{-}_{i})e_{i}$ (12) s.t.
$\displaystyle\sum_{p\in\Omega}\bar{x}_{i}^{p}\cdot z^{p}\geq 1-e_{i}$
$\displaystyle(\lambda_{i})\ \ $ $\displaystyle\forall i\in I^{-}$ (13)
$\displaystyle-\sum_{p\in\Omega}z^{p}\geq-k$ $\displaystyle(\mu)\ \ $ (14)
$\displaystyle 0\leq z^{p}\leq 1$ $\displaystyle\forall p\in\Omega$ (15)
$\displaystyle 0\leq e_{i}\leq 1$ $\displaystyle\forall i\in I^{-}$ (16)

where $\lambda_{i}$ for all $i\in I^{-}$ and $\mu$ represent the dual
variables associated to constraints (13) and (14) respectively. Constraints
(14) are relaxed without loss of quality in the bound, to improve numerical
stability.

##  4 Algorithms

The formulation (12) - (16) contains a combinatorial number of variables
$z^{p}$. In the following, we build an algorithm that encompasses the standard
column generation approach (Amaldi, Dhyani, and Ceselli 2012), using such an
extended formulation as a starting point. In addition, we propose an algorithm
to solve the pricing that is specific to the problem, as it exploits its
combinatorial structure. The reduced cost of a variable $z^{p}$ can be written
as $\bar{z}^{p}=0-\sum_{i\in I^{-}}\bar{x}_{i}^{p}\cdot\lambda_{i}+\mu$. The
pricing model (PM) is therefore the following:

$\displaystyle\text{min }\bar{z}^{p}=$ $\displaystyle-\sum_{i\in
I^{-}}x_{i}\cdot\lambda_{i}+\mu$ (17) s.t. $\displaystyle b+\sum_{j\in
J}w_{j}a_{ij}\geq 1$ $\displaystyle\forall i\in I^{+}$ (18) $\displaystyle
b+\sum_{j\in J}w_{j}a_{ij}\leq-1+\xi_{i}$ $\displaystyle\forall i\in I^{-}$
(19) $\displaystyle\xi_{i}\leq M\cdot(1-x_{i})$ $\displaystyle\forall i\in
I^{-}$ (20) $\displaystyle\xi_{i}\in[0,+\infty)$ $\displaystyle\forall i\in
I^{-}$ (21) $\displaystyle x_{i}\in\\{0,1\\}$ $\displaystyle\forall i\in
I^{-}$ (22) $\displaystyle b\in\mathbb{R}$ (23) $\displaystyle
w_{j}\in\mathbb{R}$ $\displaystyle\forall j\in J$ (24)

We explicitly indicate the dependence on the specific dual variables by
writing PM($\bm{\lambda}$, $\mu)$.

A restriction of model (12) - (16) obtained by replacing $\Omega$ with a
subset $\Omega^{\prime}\subseteq\Omega$, is called Restricted Master Model
(RMM). Given a RMM, we call Hyperplane Choice Model (HCM) the same model with
0-1 variables instead of continuous ones.

Our version of column generation with hyperplane choice is detailed in Alg. 4
of Appendix 8. At each iteration, we solve the restricted master model and the
pricing problem. We address the pricing problem applying one of two different
approaches, solving PM with general purpose solvers, or applying the algorithm
that we show later in Section 4.1. Subsequently, RMM is updated adding the
columns obtained from the solution of the pricing problem and HCP is solved to
find the best possible solution to PCAB using hyperplanes obtained from
solving pricing. This last step allows us to monitor the error value at each
iteration. The algorithm returns with the last available PCAB solution when
all the columns generated in a certain iteration have non-negative reduced
cost or when the error value of the PCAB solution is 0.

Notice that, when solving PM using general purpose solvers, we get by default
a single column. However, Alg. 4 takes into account the possibility to employ
the algorithm AHP of Section 4.1 instead, which can generate more than one
column at each iteration.

###  4.1 Ad Hoc Heuristic Algorithm for the Pricing Problem

Solving the pricing problem implies placing a single valid hyperplane. For the
solution to be optimal, we want to maximize the sum of the dual variables
associated to the negative points that are separated from the positive points
by means of that hyperplane. A subset of negative points is separable from the
set of positive points if and only if their respective convex hulls are
disjoint. Therefore, in order to find a good heuristic solution to the pricing
problem, it is sufficient to find a subset of negative points such that:

  1. 1.

the sum of the corresponding dual variables is “high”;

  2. 2.

its convex hull does not intersect the convex hull of the positive points.

The main intuition behind our ad hoc heuristic algorithm for the pricing
problem is that such a set of negative points can be found in a sequential
fashion, considering the negative points one at a time and keeping only those
that maintain the second condition true. We try to ensure the first condition
considering the negative points accordingly to their dual variable values, in
decreasing order. This approach is appealing in that it lends itself to
parallelization. In our actual implementation, we run twin parallel processes
of Alg. 1 by changing the starting point only; this increases the chances of
obtaining a good solution.

Given a subset $\bar{I}^{-}\subseteq I^{-}$ of negative points, it is possible
to check if the conxex hull of the points indexed by $\bar{I}^{-}$ intersects
the convex hull of the positive points after solving the following linear
programming model, called linear separability verification model and denoted
LSVM($\bar{I}^{-}$).

min $\displaystyle\sum_{j\in J}\left(\xi_{j}+\eta_{j}\right)$ (25) s.t.
$\displaystyle\sum_{i\in I^{+}}\alpha_{i}\cdot
a^{+}_{ij}+\xi_{j}=\sum_{i\in\bar{I}^{-}}\beta_{i}\cdot a^{-}_{ij}+\eta_{j}$
$\displaystyle\forall j\in J$ (26) $\displaystyle\sum_{i\in
I^{+}}\alpha_{i}=1$ (27) $\displaystyle\sum_{i\in\bar{I}^{-}}\beta_{i}=1$ (28)
$\displaystyle\alpha_{i}\in[0,\ 1]$ $\displaystyle\forall i\in I^{+}$ (29)
$\displaystyle\beta_{i}\in[0,\ 1]$ $\displaystyle\forall i\in\bar{I}^{-}$ (30)
$\displaystyle\xi_{j}\geq 0$ $\displaystyle\forall j\in J$ (31)
$\displaystyle\eta_{j}\geq 0$ $\displaystyle\forall j\in J$ (32)

The model admits a solution with all zero slacks, i.e. $\xi_{j}=0$ and
$\eta_{j}=0$ for every $j\in J$, if and only if it exists a convex combination
of the positive points equal to a convex combination of the given negative
points or, in other words, if their respective convex hulls intersect.

The pseudocode of our ad hoc algorithm is given in Alg. 1 and Alg. 2. Alg. 1
shows the process called _heuristic partial linear separation_ (HPLS), whose
role is to find the desired set of negative points and to derive a single
pricing solution, given by an indicator vector $\bm{x}$ and hyperplane
parameters $b$ and $\bm{w}$. It requires as input the index $i^{start}\in
I^{-}$ of a starting negative point, the maximum number $thr$ of negative
points to consider and a permutation $\sigma:I^{-}\rightarrow I^{-}$ such that
$\lambda_{\sigma(i)}\geq\lambda_{\sigma(i^{\prime})}$ for every
$i,i^{\prime}\in I^{-}$ with $i<i^{\prime}$.

Algorithm 1 Heuristic Partial Linear Separation

Input: permutation $\sigma$, starting index $i^{start}$, point threshold $thr$

$\bar{I}^{-}\leftarrow\\{i^{start}\\}$

$counter\leftarrow 0$

while $counter<n$ and $counter<thr$ do

if $\sigma(counter)=i^{start}$ then

$counter\leftarrow counter+1$

continue

end if

Solve LSVM($\bar{I}^{-}$) and get slack variable values $\xi_{j}$ and
$\eta_{j}$ for every $j\in J$

if $\exists j\in J$ such that $\left(\xi_{j}>0\ \textbf{or}\
\eta_{j}>0\right)$ then

$\bar{I}^{-}\leftarrow\bar{I}^{-}\cup\\{\sigma(counter)\\}$

end if

end while

Solve a linear separation problem between all the positive points and the
negative points indexed by $\bar{I}^{-}$ and get hyperplane coefficients $b$
and $\bm{w}$

Shift the hyperplane as close as possible to the set of positive points and
get new coefficient b

$\bm{x}\leftarrow\bm{0}\in\mathbb{R}^{n}$

for $i\in I^{-}$ do

if $b+\langle\bm{w},\ \bm{a}^{-}_{i}\rangle\ \leq-1$ then

$x_{i}\leftarrow 1$

end if

end for

return $\bm{x}$, $b$, $\bm{w}$

Until all the negative points have been considered or the point threshold
$thr$ has been reached, a new point is added to $\bar{I}^{-}$ and the model
LSVM($\bar{I}^{-}$) is solved to determine whether to keep or to neglect the
point. After the end of the loop, we solve a linear separation problem to find
an equation $(b,\bm{w})$ of a hyperplane that separates the negative points
indexed by the final $\bar{I}^{-}$ from the positive points. This hyperplane
is then moved as close as possible to the set of positive points, in order to
possibly increase the number of negative points that it classifies correctly.
The indicator vector $x$ is obtained assigning 1 to every negative point in
$\bar{I}^{-}$ and 0 to the other ones.

Algorithm 2 Ad Hoc Heuristic Algorithm for the Pricing Problem

Input: dual variables $\bm{\lambda}=(\lambda_{i})_{i\in I^{-}}$ and $\mu$,
point threshold $thr$, vector $\bm{x}^{prev}\in\\{0,1\\}^{n}$, maximum number
of processes $N^{max}$

Find a permutation $\sigma:I^{-}\rightarrow I^{-}$ such that
$\lambda_{\sigma(i)}\geq\lambda_{\sigma(i^{\prime})}$ for every
$i,i^{\prime}\in I^{-}$ with $i<i^{\prime}$

Set $\tilde{x}_{i}=0$ for every $i\in I^{-}$

$\bm{p}\leftarrow\bm{x}^{prev}\circ\bm{\lambda}$

$p_{i}\leftarrow p_{i}/\sum_{j\in I^{-}}p_{j}$ for every $i\in I^{-}$

Draw $i^{start}$ from $I^{-}$ according to distribution $\bm{p}$

if $\exists i\in I^{-}$ such that $\bar{x}_{i}^{prev}=0$ then

$\bar{y}_{i}^{prev}\leftarrow 1-\bar{x}_{i}^{prev}$ for every $i\in I^{-}$

else

$\bar{y}_{i}^{prev}\leftarrow 1$ for every $i\in I^{-}$

end if

$\bm{p}\leftarrow\bar{\bm{y}}^{prev}\circ\bm{\lambda}$

$p_{i}\leftarrow p_{i}/\sum_{j\in I^{-}}p_{j}$ for every $i\in I^{-}$

$N\leftarrow 0$

while $N<N^{max}$ and $\exists i\in I^{-}$ such that $\tilde{x}_{i}=0$ do

Run parallel HPLS($\sigma$, $i^{start}$, $thr$) process

$N\leftarrow N+1$

if $N=N^{max}$ then

break

end if

Draw $i^{start}$ from $I^{-}$ according to distribution $\bm{p}$

end while

Wait for the running HPLS processes to terminate and collect $\bm{x}^{k}$,
$\bm{w}^{k}$ and $b^{k}$ for every process $k\in\bar{N}$

return $N$ and $\bm{x}^{k}$, $b^{k}$, $\bm{w}^{k}$ for every $k\in\bar{N}$

Alg. 2 shows the general structure of the ad hoc heuristic for pricing (AHP),
whose task is to obtain multiple pricing solutions by running parallel HPLS
processes. The following input parameters are required: $\bm{\lambda}$ and
$\mu$ are the dual variables obtained from the solution of the most recent
restricted master model, $thr$ is the point threshold to be used in HPLS,
$\bm{x}^{prev}$ is the indicator vector corresponding to the pricing solution
with lowest reduced cost found in the previous iteration of column generation
and finally, $N^{max}$ is the maximum allowed number of parallel instances of
HPLS. Notice that, according to Alg. 4, at the first iteration of column
generation, $\bm{x}^{prev}=\bm{1}$.

As already mentioned, every HPLS process requires a starting index
$i^{start}$. Such starting indices are chosen at random with probability
proportional to the corresponding values of the dual variables $\bm{\lambda}$,
in order to favour points with higher dual variable values. For the first
parallel process, only indices with a corresponding 1 in the indicator vector
$\bm{x}^{prev}$ can be chosen. For the subsequent processes, it holds the
opposite: an index can be chosen only if the corresponding value in
$\bm{x}^{prev}$ is 0. In other words, for the first process we select the
index of a point that, in the previous iteration of column generation, has
been separated from the positive points by the best hyperplane (one with
lowest corresponding reduced cost among the ones obtained); for the other
processes, we choose only indices of negative points that were erroneously
classified as positive by the same hyperplane. This differentiated behaviour
allows us to try an improvement mimicking the best pricing solution of the
previous iteration and, at the same time, to put effort in exploring different
solutions.

###  4.2 A Greedy Matheuritic for PCAB

We devised also a greedy matheuristic approach to solve PCAB , for both
obtaining an initial solution quickly, and to populate the initial RMP. It
consists in placing a single hyperplane at a time and ignoring every negative
point that is already separated correctly by at least one of the hyperplanes
created so far. To find a hyperplane that, given a certain subset
$\bar{I}^{-}\subseteq I^{-}$ of negative points, separates the positive points
from the greatest number of negative points in $\bar{I}^{-}$, it is possible
to solve the model below, called _Partial Linear Separation Model_.

max $\displaystyle\sum_{i\in\bar{I}^{-}}c(\bm{a}^{-}_{i})x_{i}$ (33) s.t.
$\displaystyle b+\sum_{j\in J}w_{j}a_{ij}\geq 1$ $\displaystyle\forall i\in
I^{+}$ (34) $\displaystyle b+\sum_{j\in J}w_{j}a_{ij}\leq-1+\xi_{i}$
$\displaystyle\forall i\in\tilde{I}^{-}$ (35) $\displaystyle\xi_{i}\leq
M\cdot(1-x_{i})$ $\displaystyle\forall i\in\tilde{I}^{-}$ (36)
$\displaystyle\xi_{i}\in[0,+\infty)$ $\displaystyle\forall i\in\tilde{I}^{-}$
(37) $\displaystyle x_{i}\in\\{0,1\\}$ $\displaystyle\forall
i\in\tilde{I}^{-}$ (38) $\displaystyle b\in\mathbb{R}$ (39) $\displaystyle
w_{j}\in\mathbb{R}$ $\displaystyle\forall j\in J$ (40)

We write PLSM($\bar{I}^{-}$) to denote the model and its dependence on the set
$\bar{I}^{-}$. This model is similar to (17) - (24), the pricing model of our
column generation. The differences are that not necessarily all the negative
points are involved and that the objective is a particular case of (17), in
which $\lambda_{i}=1$ for all $i$ and $\mu=0$.

In Alg. 3 we show a pseudocode for the greedy algorithm.

Algorithm 3 Greedy Algorithm

Set $x_{i}=0$ for every $i\in I^{-}$

$\bar{I}^{-}\leftarrow I^{-}$

$t\leftarrow 0$

while $t<\mathcal{K}$ and $\exists i\in I^{-}$ such that $x_{i}=0$ do

Solve PLSP($\bar{I}^{-}$) and get $\\{x_{i}^{t}\\}_{i\in\bar{I}^{-}}$, $b^{t}$
and $\bm{w}^{t}=(w_{j}^{t})_{j\in J}$

Set $x_{i}=1$ for every $i\in\bar{I}^{-}$ such that $x_{i}^{t}=1$

$\bar{I}^{-}\leftarrow\\{i\in\bar{I}^{-}\ |\ x_{i}=0\\}$

$t\leftarrow t+1$

end while

return $\\{(\bm{w}^{\tau},\ b^{\tau})\ |\ \tau\in\\{1,\ldots,t\\}\\}$

Binary variables $x_{i}$ indicate whether the negative point $i$ is classified
correctly by at least one chosen hyperplane or not. The set $\bar{I}^{-}$ is
initialized to contain all the negative points. As long as less that
$\mathcal{K}$ hyperplanes have been placed and there is still some
misclassified negative points, a hyperplane is chosen solving
PLSM($\bar{I}^{-}$) and all the points classified correctly by that hyperplane
are removed from $\bar{I}^{-}$.

##  5 Experimental Setting and Results

We create two datasets (D1 and D2), each composed by 3 sets of 100 instances
each, with number of dimensions $d$ equal to 2, 4 and 8 respectively, assuming
the unitary hypercube as the feasible region and generating an exponential
number of points with respect to the dimension (ranging from about 200 points
for $d=2$ to 10000 points for $d=8$). D1 and D2 are generated similarly, but
instances in D1 have a known minimum separation budget, whereas instances in
D2 do not. Details on the generations of our instances are reported in
Appendix 10; both datasets are available online (Messana 2024). The
experimental analysis is carried out in two steps. In the first part, we use
dataset D1 and we consider our compact models (47) - (54) (named MODEL A) and
(1) - (10) (MODEL B), column generation with original pricing formulation
(COLGEN 1) and with ad hoc pricing algorithm (COLGEN 2), and finally the
greedy algorithm of Section 4.2 (GREEDY). We also adapted a compact model
proposed in (Orsenigo and Vercellis 2007) for convex polyhedral
classification; the original model was designed for a problem closely matching
ours, and was therefore easy to adapt. Such an adaptation, that we name OV
2007, is detailed in Appendix 9.

In the second part of the analysis, we investigate the capability of a column
generation approach to approximate the convex hull of the positive points,
using dataset D2. In particular we consider COLGEN 2 and we compare it against
a procedure customized from the state-of-the art software polymake, Assarf et
al. (2017). For this second part of the experiments, we use very similar
datasets as in the first part, but with unknown minimum separation budget.

All the programming models, including the ones involved in column generation
and in the greedy algorithm, have been tackled using Gurobi. The code was
written and executed in Python 3.9 on a machine with Ubuntu 20.04.6 LTS, 32 GB
DDR4 RAM and processor Intel(R) Xeon(R) W-1250P CPU @ 4.10GHz.

###  5.1 Implementation Choices

For Model B, we fixed $M=10000$, which proved to be a sufficiently large value
not to interfere negatively with the solution process. For OV 2007, we tried
different values for $\epsilon$, including machine precision as suggested by
the authors in the original paper. However, despite several attempts, we have
not been able to find a value that always ensures a correct behaviour in the
solution process and in the end we set $\epsilon=0.001$. For further details,
see Section 5.2.

Both the column generation algorithms considered are heuristic. The first one
is Alg. 4 with PM($\bm{\lambda}$, $\mu$) solved using Gurobi and interrupting
the solver inside the branch-and-bound algorithm immediately after the root
node. The reason for this choice is twofold: preliminary tests showed us that
the solution process of PM terminates within a reasonable time only on very
small datasets; moreover, the value of the objective function at the root node
is often close to the values obtained going on in the branch-and-bound. The
second column generation algorithm has the same structure of the first one,
but AHP (Alg. 2) is executed instead of solving PM. We have chosen the maximum
number of parallel processes $N^{max}$ equal to 8. The point threshold $thr$
has been set equal to $d$.

In the greedy algorithm, like in the first column generation, the solution
process of PLSM is interrupted right after the root node of the branch-and-
bound tree for the same reasons.

For all the methods, the weight function (that is, $\hat{c}$ for Model A and
$c$ for the others) has been kept constant and equal to 1 on every
$a\in\mathbb{R}^{d}$.

###  5.2 Comparing Models and Algorithms

We carried out preliminary experiments using 8 out of 100 instances for each
of the three dimensions using large time limits, in order to have feedback
about the performance of the models and the algorithms tested on the long run.
Subsequently, we ran tests using all the 100 instances with lower time limits.

#### Preliminary experiments

For each $d\in\\{2,4,8\\}$, we tested the six methods using three different
hyperplane budgets: $d$, $2d-1$ and $2d$. As already mentioned in Section 10,
$\mathcal{K}=2d$ represents the minimum separation budget. Each method has
been run with a time limit equal to $6\cdot 10^{log_{2}(d)}$ seconds (i.e.
$10^{log_{2}(d)-1}$ minutes).

Figure 2: Error value with 8-dimensional hypercubic data: $\mathcal{K}=8$
(sample instance). Figure 3: Error value with 8-dimensional hypercubic data:
$\mathcal{K}=15$ (sample instance). Figure 4: Error value with 8-dimensional
hypercubic data: $\mathcal{K}=16$ (sample instance).

In order to give a representative example, Figures 2 \- 4 show the error value
over time for a single instance in dimension 8. Analogous figures are present
in Appendix 11 for dimensions 2 and 4. Both error and time are reported in
log-10 scale. In the legends, next to the method names, their respective final
error values are indicated in brackets. In the following, we state the main
observations that we can derive from the figures, but very similar comments
hold for all the instances in the preliminary experiments, according to their
dimensionality.

First of all, we can notice that the performance of Model A in terms of error
value is the worst among all the methods in most of the cases the unique
exception is the experiment whose results are reported Fig. 4). So we can
affirm that generalizing in a straightforward way a popular and vastly applied
model such as SVM does not yield good results for our problem. Also notice
that the lines for Model A are non-monotone, because the objective that it
tries to minimize is different from the error value. The column generation
algorithms yield the best performances overall. This holds especially for
COLGEN 2, which converges very quickly and obtains the best error values in
the majority of the cases. Model B and GREEDY give intermediate results.

For the instance in dimension 2, all the methods converge really quickly. In
dimension 4 and 8, we can notice more differentiated convergence speeds. We
can also see that the extent of the difference range for error value and
convergence speed seems to be directly proportional to the hyperplane budget.
Finally, the reader may notice that OV 2007 is not present in Figures 2 \- 4.
The reason is that the solutions obtained from solving that model are not
valid: the solver outputs $b=0$ and $\bm{w}=\bm{0}$ as a computational
consequence of the fact that we were not able to find a value of the
hyperparameter $\epsilon$ that allows the solver to find a proper solution
when $d=8$.

#### Massive experiments

In the preliminary experiments we were able to verify that one tenth of the
given time limit was always sufficient for the column generation algorithms to
reach approximately the final error value. Therefore, we tested the six
methods on all the 100 instances using the same hyperplane budgets as in the
preliminary experiments, but time limit equal to $6\cdot 10^{log_{2}(d)-1}$
for each $d\in\\{2,\ 4,\ 6\\}$. Tables 1 \- 3 show the results averaged over
all the instances. We report running time, final error value and time needed
to reach such error value (this latter abbreviatd to TE, time up to error),
for each of the three dimensions respectively. Additional tables can be found
in Appendix 11. They contain details about the performance of the two column
generation algorithms.

$\mathcal{K}$ | method | time (seconds) | % error value | TE (seconds)
---|---|---|---|---
2 | COLGEN 1 | 6.044 | 22.28 | 0.625
COLGEN 2 | 6.052 | 22.66 | 2.747
GREEDY | 0.164 | 27.23 | 0.164
MODEL A | 6.072 | 33.87 | 3.843
MODEL B | 6.081 | 21.24 | 0.856
OV 2007 | 6.190 | 22.47 | 4.253
3 | COLGEN 1 | 6.040 | 7.87 | 0.767
COLGEN 2 | 6.052 | 6.52 | 3.304
GREEDY | 0.199 | 9.86 | 0.198
MODEL A | 6.102 | 13.36 | 3.876
MODEL B | 6.109 | 6.07 | 1.774
OV 2007 | 6.243 | 7.40 | 3.358
4 | COLGEN 1 | 1.408 | 0.11 | 0.711
COLGEN 2 | 0.474 | 0.00 | 0.472
GREEDY | 0.208 | 0.92 | 0.207
MODEL A | 5.500 | 1.52 | 4.179
MODEL B | 0.220 | 0.00 | 0.215
OV 2007 | 1.460 | 0.08 | 1.388
Table 1: Average running time, error value and final error hitting time for $d=2$. $\mathcal{K}$ | method | time (seconds) | % error value | TE (seconds)
---|---|---|---|---
4 | COLGEN 1 | 60.108 | 12.48 | 8.409
COLGEN 2 | 60.457 | 12.07 | 29.835
GREEDY | 2.111 | 15.89 | 2.108
MODEL A | 60.410 | 24.66 | 52.789
MODEL B | 60.408 | 13.49 | 25.264
OV 2007 | 60.911 | 22.37 | 51.594
7 | COLGEN 1 | 58.990 | 1.96 | 7.568
COLGEN 2 | 60.142 | 1.44 | 4.159
GREEDY | 2.232 | 2.81 | 2.227
MODEL A | 60.702 | 13.17 | 57.571
MODEL B | 60.686 | 1.76 | 14.369
OV 2007 | 61.501 | 11.14 | 49.668
8 | COLGEN 1 | 32.357 | 0.36 | 7.760
COLGEN 2 | 1.185 | 0.00 | 1.180
GREEDY | 2.267 | 1.22 | 2.260
MODEL A | 60.781 | 15.92 | 55.100
MODEL B | 18.662 | 0.10 | 14.873
OV 2007 | 61.696 | 8.43 | 52.282
Table 2: Average running time, error value and final error hitting time for $d=4$. $\mathcal{K}$ | method | time (seconds) | % error value | TE (seconds)
---|---|---|---|---
8 | COLGEN 1 | 601.465 | 11.20 | 95.704
COLGEN 2 | 746.987 | 11.08 | 517.667
GREEDY | 30.880 | 12.80 | 30.795
MODEL A | 618.026 | 60.33 | 442.587
MODEL B | 617.558 | 17.68 | 581.188
OV 2007 | 637.087 | - | 35.709
15 | COLGEN 1 | 433.135 | 2.87 | 108.615
COLGEN 2 | 603.572 | 1.27 | 31.553
GREEDY | 33.963 | 5.22 | 33.487
MODEL A | 633.393 | 33.83 | 46.513
MODEL B | 632.671 | 10.25 | 455.574
OV 2007 | 666.978 | - | 65.929
16 | COLGEN 1 | 288.734 | 2.66 | 109.489
COLGEN 2 | 26.953 |  $<$ 10-4 | 26.717
GREEDY | 34.245 | 4.78 | 33.687
MODEL A | 635.822 | 30.20 | 337.217
MODEL B | 634.875 | 10.20 | 506.273
OV 2007 | 671.311 | - | 70.104
Table 3: Average running time, error value and final error hitting time for
$d=8$.

The massive experiments confirm the relative performances showed by the
methods in the preliminary experiments. In dimension 2, the small size of the
instances allows MODEL B and OV 2007 to reach error values lower or comparable
to those obtained by column generation. In dimension 4, MODEL B, COLGEN 1 and
COLGEN 2 all give small error values but we can notice that the time up to
error can be significantly higher for MODEL B. Finally, in dimension 8, COLGEN
2 proves to be the method with best error values, while COLGEN 1 is slightly
less performing. In decreasing order of final error, the other methods are
GREEDY, MODEL B, MODEL A and OV 2007. Notice that greedy is the only method
that always terminates within the given time limit, thanks to the stopping
criterion of its iteration step. For the other methods, reaching optimality
within the time limit turns out to be difficult except when $d$ is not too
large and $\mathcal{K}=2d$. However, both column generation algorithms
terminate within the time limit when $\mathcal{K}=2d$ for every
$d\in\\{2,4,8\\}$.

Summarizing, the experiments highlight COLGEN 2 as the most promising method
overall. However, it is at the same time the most complex. Parallelization
showed to be an essential component. Its effective use require the fine tuning
of the point threshold _thr_ , which might not be easy in all applications.

###  5.3 Producing Polyhedral Approximations

In a second round of experiments we assessed the effectiveness of our
algorithms in producing good convex hull approximations. We considered COLGEN
2, being the most effective algorithm in previous experiments, with a time
limit set as in the previous section. As a benchmark we considered a procedure
adapted from Polymake. Polymake implements convex hull algorithms, possibly
producing way more than $k$ hyperplanes. We considered a greedy approach
(GREEDY POLYMAKE): firstly, it computes the equations of the hyperplanes
describing the facets of the convex hull of the positive points; secondly, it
iteratively selects a hyperplane that correctly classifies the highest number
of negative points and it removes those points from the dataset. An additional
hyperplane gets selected with the same criterion, until the budget
$\mathcal{K}$ is reached. A larger time limit of one hour was given to GREEDY
POLYMAKE.

We compared COLGEN 2 and GREEDY POLYMAKE on all the 100 instances generated
for this second part of the experiments, applying the same hyperplane budget
as in part one. The metrics for the comparison are the error, i.e. the number
of misclassified negative point, and the volume of the convex region delimited
by the resulting hyperplanes. In the cases in which the hyperplanes generated
by an algorithm for a specific instance describe an unbounded polyhedron, we
consider the region bounds for the generated points, i.e. $-1\leq x_{i}\leq 2$
for every $i\in\\{1,\ldots,d\\}$. The computation of the convex hulls for
GREEDY POLYMAKE and the computation of the volume of the resulting polyhedron
for both the methods were achieved using ake (Gawrilow and Joswig 2000,
Assarf, Gawrilow, Herr, Joswig, Lorenz, Paffenholz, and Rehn 2017). Tables 4,
5 and 6 show the results in dimension 2, 4 and 8 respectively, averaged over
the 100 instances for each hyperplane budget and each algorithm. Table 6
doesn’t contain those results that require to use ake, as GREEDY POLYMAKE was
always hitting the time limit (without producing any valid solution) also
computing the volume was extremely computational expensive in dimension 8.

$\mathcal{K}$ | method |  % error value |  volume
---|---|---|---
2 | COLGEN 2 | 22.01 | 3.00
GREEDY POLYMAKE | 21.52 | 2.96
3 | COLGEN 2 | 5.30 | 1.69
GREEDY POLYMAKE | 9.44 | 1.91
4 | COLGEN 2 | 0.00 | 1.00
GREEDY POLYMAKE | 0.06 | 0.98
Table 4: Average error value and volume of the approximate convex hull for $d=2$. $\mathcal{K}$ | method |  % error value |  volume
---|---|---|---
4 | COLGEN 2 | 5.33 | 7.64
GREEDY POLYMAKE | 6.93 | 9.00
7 | COLGEN 2 | 0.05 | 2.11
GREEDY POLYMAKE | 0.35 | 1.98
8 | COLGEN 2 | 0.00 | 1.71
GREEDY POLYMAKE | 0.04 | 1.37
Table 5: Average error value and volume of the approximate convex hull for $d=4$. $\mathcal{K}$ | method |  % error value |  volume
---|---|---|---
8 | COLGEN 2 | 1.07 | -
GREEDY POLYMAKE | - | -
15 | COLGEN 2 | 0.01 | -
GREEDY POLYMAKE | - | -
16 | COLGEN 2 | $<0.003$ | -
GREEDY POLYMAKE | - | -
Table 6: Average error value and volume of the approximate convex hull for
$d=8$.

Regarding dimension 2 and 4, the column generation algorithm obtains lower
average separation error in 5 cases over 6. Because the hyperplanes describing
the convex hull have distance 0 from the set of positive points, we could
expect an advantage of GREEDY POLYMAKE in terms of the volume of the resulting
polytopes. Instead, COLGEN 2 obtains a better result in 2 cases over 6 and its
volume values are close to those of GREEDY POLYMAKE in the remaining cases.

##  6 Conclusions

We proposed a mathematical programming approach to a problem in computational
geometry, which is to find a polyhedron with a limited number of faces,
approximating the convex hull of a set of points (PCAB in the paper).

Ultimately, our approach allows us to design formulations for the PCAB , and
ad-hoc algorithms which exploit column generation techniques, and rely on
special purpose pricing routines.

The first task was to design formulations with a polynomial number of
variables and constraints. Although our study began by inspirations from
support vector machines, experiments indicate a direct adaptation not to be
competitive. By elaborating on them, however, we could propose compact models
which outperformed earlier MIPs from the literature. Key ideas were to
properly change the penalty function, removing both the handling of margins
and the hyperplane norm regularization terms.

The second task was to design more effective resolution algorithms. We used
decomposition methods to obtain extended formulations and a column generation
scheme. When solving pricing problems as MIPs, unfortunately, such a
decomposition process had only a partial impact: even if smaller, pricing MIPs
remain challenging to tackle with general purpose solvers. We designed an ad-
hoc heuristic pricing routine that takes advantage from considering only a
subset of points simultaneously and from parallelization. Its embedding in the
column generation framework led to an effective algorithm for PCAB . In our
experiments, in its best configuration, the column generation algorithm was
able to solve instances up to dimension 8 with number of points up to over 10
thousand, comparing favourably to all the benchmark algorithms we have
considered.

With respect to the target of obtaining convex hull approximations, our
experiments highlight our approach to be effective. We mention three appealing
features. First, our algorithms are flexible tools for numerical resolution;
their structure makes it easy to tune them to be more accurate (at the price
of higher computing times) or more quick (at the price of a few more
separation errors). Standard approaches from the literature do not have such
an option. Second, our models allow to possibly produce unbounded polyhedra as
solutions; convex hull algorithms are forced to produce convex polytopes.
Third, our models assume the number of hyperplanes, to be produced as faces,
as part of data. This may actually become a hyperparameter in applications.
That is, our methods allow parametric analyses to be performed, to evaluate
the convex hull approximation error which is made by progressively allowing
less faces, and therefore producing less complex polyhedra. It is interesting
to note that, in our experiments, very few hyperplanes are enough to get
approximations whose volume closely match that of the actual convex hull.

As a benchmark, we considered adaptations of convex hull algorithms from the
literature, and in particular from the package Polymake; these latter show to
be faster than our mathematical programming approaches when the number of
dimensions is very low (providing approximations of similar quality), but get
orders of magnitude slower as the dimensionality increases. The use of
Polymake becomes inapplicable already in dimension $8$. We report that both
approaches are affected by the number of points in the instance, which may
however increase exponentially (in case the positive and negative points are
produced by a proper sampling) or not (in case the positive and negative
points are feasible and infeasible solutions of a MIP).

The set of In our PCAB , the sets of positive and negative points are assumed
to be given. As perspective research, we plan to consider the case in which
new positive and negative points can be generated by suitable queries, at a
cost, as in Messana et al. (2024). We expect this line of research to yield
methods which are applicable in even wider settings.

Acknowledgments.  The work has been partially supported by the University of
Milan, Piano Sostegno alla Ricerca.

## References

  * Amaldi et al. (2012) Amaldi E, Dhyani K, Ceselli A (2012) Column generation for the minimum hyperplanes clustering problem 25.
  * Assarf et al. (2017) Assarf B, Gawrilow E, Herr K, Joswig M, Lorenz B, Paffenholz A, Rehn T (2017) Computing convex hulls and counting integer points with polymake. _Mathematical Programming Computation_ 9:1–38.
  * Astorino et al. (2023) Astorino A, Avolio M, Fuduli A (2023) Maximum-margin polyhedral separation for binary multiple instance learning. _EURO Journal on Computational Optimization_ .
  * Astorino et al. (2021) Astorino A, Francesco MD, Gaudioso M, Gorgone E, Manca B (2021) Polyhedral separation via difference of convex (dc) programming. _Soft Computing_ 25:12605 – 12613.
  * Astorino and Fuduli (2015) Astorino A, Fuduli A (2015) Support vector machine polyhedral separability in semisupervised learning. _Journal of Optimization Theory and Applications_ 164:1039–1050.
  * Astorino and Gaudioso (2002) Astorino A, Gaudioso M (2002) Polyhedral separability through successive lp. _Journal of Optimization Theory and Applications_ 112:265–293.
  * Balestriero et al. (2022) Balestriero R, Wang Z, Baraniuk R (2022) Deephull: fast convex hull approximation in high dimensions. _ICASSP 2022 - 2022 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ 3888–3892.
  * Baum (1990) Baum EB (1990) A polynomial time algorithm that learns two hidden unit nets. _Neural Computation_ 2(4):510–522.
  * Baum (1991) Baum EB (1991) Neural net algorithms that learn in polynomial time from examples and queries. _IEEE Transactions on Neural Networks_ 2(1):5–19.
  * Blum et al. (1995) Blum A, Chalasani P, Goldman SA, Slonim DK (1995) Learning with unreliable boundary queries. _Proceedings of the eighth annual conference on Computational learning theory_ , 98–107.
  * Bradford Barber et al. (1996) Bradford Barber C, Dobkin DP, Huhdanpaa H (1996) The quickhull algorithm for convex hulls. _ACM Trans. Math. Softw._ 22:469–483.
  * Bykat (1978) Bykat A (1978) Convex hull of a finite set of points in two dimensions. _Inf. Process. Lett._ 7:296–298, URL https://api.semanticscholar.org/CorpusID:31239431.
  * Carrizosa et al. (2010) Carrizosa E, Martin-Barragan B, Romero Morales D (2010) Binarized support vector machines. _INFORMS Journal on Computing_ 22(1):154–167.
  * Eddy (1977) Eddy WF (1977) A new convex hull algorithm for planar sets. _ACM Trans. Math. Softw._ 3:398–403, URL https://api.semanticscholar.org/CorpusID:6464920.
  * Fajemisin et al. (2024) Fajemisin AO, Maragno D, den Hertog D (2024) Optimization with constraint learning: A framework and survey. _European Journal of Operational Research_ 314(1):1–14.
  * Gawrilow and Joswig (2000) Gawrilow E, Joswig M (2000) Polymake: a framework for analyzing convex polytopes. _Polytopes—combinatorics and computation_ , 43–73 (Springer).
  * Gottlieb et al. (2018) Gottlieb LA, Kaufman E, Kontorovich A, Nivasch G (2018) Learning convex polytopes with margin. _Neural Information Processing Systems_.
  * Jayaram and Fleyeh (2016) Jayaram M, Fleyeh H (2016) Convex hulls in image processing: a scoping review. _American Journal of Intelligent Systems_ 6(2):48–58.
  * Kantchelian et al. (2014) Kantchelian A, Tschantz MC, Huang L, Bartlett PL, Joseph AD, Tygar JD (2014) Large-margin convex polytope machine. _Neural Information Processing Systems_.
  * Klivans et al. (2004) Klivans AR, O’Donnell R, Servedio RA (2004) Learning intersections and thresholds of halfspaces. _Journal of Computer and System Sciences_ 68(4):808–840.
  * Kwek and Pitt (1996) Kwek S, Pitt L (1996) Pac learning intersections of halfspaces with membership queries. _Proceedings of the ninth annual conference on Computational learning theory_ , 244–254.
  * Leng et al. (2020) Leng Q, He Z, Liu Y, ping Qin Y, Li Y (2020) A soft-margin convex polyhedron classifier for nonlinear task with noise tolerance. _Applied Intelligence_ 51:453–466.
  * Lombardi and Milano (2018) Lombardi M, Milano M (2018) Boosting combinatorial problem modeling with machine learning. _Proceedings of the Twenty-Seventh International Joint Conference on Artificial Intelligence_ .
  * Manwani and Sastry (2010) Manwani N, Sastry PS (2010) Learning polyhedral classifiers using logistic function. _Asian Conference on Machine Learning_.
  * Megiddo (1988) Megiddo N (1988) On the complexity of polyhedral separability. _Discrete and Computational Geometry_ 3\.
  * Messana (2024) Messana R (2024) Replication Data for: Mathematical programming algorithms for convex hull approximation with a hyperplane budget. URL http://dx.doi.org/10.13130/RD_UNIMI/IE2MX0.
  * Messana et al. (2024) Messana R, Chen R, Lodi A (2024) Actively learning combinatorial optimization using a membership oracle. _ArXiv_ abs/2405.14090.
  * Orsenigo and Vercellis (2007) Orsenigo C, Vercellis C (2007) Accurately learning from few examples with a polyhedral classifier. _Computational Optimization and Applications_ 38:235–247.
  * Pawlak and Krawiec (2017) Pawlak TP, Krawiec K (2017) Automatic synthesis of constraints from examples using mixed integer linear programming. _Eur. J. Oper. Res._ 261:1141–1157.
  * Raviv et al. (2018) Raviv D, Hazan T, Osadchy M (2018) Hinge-minimax learner for the ensemble of hyperplanes. _J. Mach. Learn. Res._ 19:62:1–62:30.
  * Schede et al. (2019) Schede E, Kolb S, Teso S (2019) Learning linear programs from data. _2019 IEEE 31st International Conference on Tools with Artificial Intelligence (ICTAI)_ 1019–1026.
  * Strekalovsky et al. (2015) Strekalovsky AS, Gruzdeva TV, Orlov AV (2015) On the problem polyhedral separability: a numerical solution. _Automation and Remote Control_ 76:1803–1816.
  * Vempala (1997) Vempala S (1997) A random sampling based algorithm for learning the intersection of half-spaces. _Proceedings 38th Annual Symposium on Foundations of Computer Science_ , 508–513 (IEEE).

##  7 Classical SVM Models and Adaptations

In order to introduce our MIPs, we report below a classical soft margin SVM
model shown below.

min $\displaystyle\sum_{i\in I}C_{i}\xi_{i}-\sum_{j\in J}w_{j}^{2}$ (41) s.t.
$\displaystyle b+\sum_{j\in J}w_{j}a_{ij}\geq 1-\xi_{i}$ $\displaystyle\forall
i\in I^{+}$ (42) $\displaystyle b+\sum_{j\in J}w_{j}a_{ij}\leq-1+\xi_{i}$
$\displaystyle\forall i\in I\setminus I^{+}$ (43) $\displaystyle\xi_{i}\geq 0$
$\displaystyle\forall i\in I$ (44) $\displaystyle b\in\mathbb{R}$ (45)
$\displaystyle w_{j}\in\mathbb{R}$ $\displaystyle\forall j\in J$ (46)

Variables $\bm{w}$ encode hyperplane coefficients; variables $\xi_{i}$ encode
misclassification error on point $i\in I$. Constraints (42) and (43) impose
correct classification of the points in both positive and negative classes,
with the possibility of misclassifying some of the points paying in terms of
the corresponding violation variables. The objective is to minimize the sum of
the violations minus a regularization term proportional to the squared 2-norm
of the hyperplane coefficients $\bm{w}$. Each constants $C_{i}\geq 0$ is a
weight that quantifies the importance of correctly classifing the $i^{th}$
point. In practice, all such weights are often set equal to a same constant
$C$.
We designed the following adaptation to PCAB .

min $\displaystyle\sum_{i\in
I^{-}}\hat{c}(\bm{a}^{-}_{i})\xi_{i}^{k}x_{i}^{k}$ (47) s.t.
$\displaystyle\sum_{k\in K}x_{i}^{k}=1$ $\displaystyle\forall i\in I^{-}$ (48)
$\displaystyle b^{k}+\sum_{j\in J}w_{j}^{k}a_{ij}\geq 1$ $\displaystyle\forall
i\in I^{+},\forall k\in K$ (49) $\displaystyle b^{k}+\sum_{j\in
J}w_{j}^{k}a_{ij}\leq-1+\xi_{i}^{k}$ $\displaystyle\forall i\in I^{-},\forall
k\in K$ (50) $\displaystyle\xi_{i}^{k}\geq 0$ $\displaystyle\forall i\in
I^{-},\forall k\in K$ (51) $\displaystyle x_{i}^{k}\in\\{0,1\\}$
$\displaystyle\forall i\in I^{-},\forall k\in K$ (52) $\displaystyle
b^{k}\in\mathbb{R}$ $\displaystyle\forall k\in K$ (53) $\displaystyle
w_{j}^{k}\in\mathbb{R}$ $\displaystyle\forall j\in J,\forall k\in K$ (54)

The model prescribes to find $\mathcal{K}$ hyperplane equations, encoded by
variables $w_{j}^{k}$ with $j\in J$ and $k\in K$. A set of non-negative
variables $\xi_{i}^{k}$ encode the distance from point $i\in I^{-}$ to
hyperplane $k$, if $i$ is on the wrong side of the hyperplane $k$; they take
value $0$ otherwise, according to constraints (50). Since PCAB requires every
hyperplane in a solution to be valid, we neglect the violation variables
corresponding to positive points, imposing them to be always on the correct
side of each hyperplane, according to constraints (49). We introduce binary
variables $x_{i}^{k}$ for $i\in I^{-}$ and $k\in K$. When $x_{i}^{k}=1$, it
means that point $i$ is assigned to the $k^{th}$ hyperplane. According to
constraints (48), each negative point must be assigned to exactly one
hyperplane. Constraints (49) force the hyperplanes to be valid and constraints
(50) require the hyperplanes to keep the negative points on a same side with
possibility of violation. The objective function sums up the weighted
violations over the negative points with respect to the hyperplanes they are
assigned to, where the weights are expressed as values of a function
$\hat{c}:\mathbb{R}^{d}\to[0,+\infty)$. In force of constraints (48), if we
fix the $\mathcal{K}$ hyperplanes, the minimum value of the objective function
is the sum over the negative points of the minimum weighted violation with
respect to the hyperplanes.

##  8 Column Generation and Hyperplane Choice

Algorithm 4 Column Generation and Hyperplane Choice

Input: point limit $thr$ and maximum number of processes $N^{max}$ (only when
using AHP, see Section 4.1)

$\bar{N}\leftarrow 0$

$R\leftarrow-\infty$

$E\leftarrow+\infty$

$\bm{x}^{prev}\leftarrow\bm{1}\in\mathbb{N}^{n}$

Set $\bm{W}$ to empty list

Set $\bm{b}$ to empty list

Initialize RMM without $z$ variables

while $R<0$ and $E>0$ do

Solve RMM and get dual variable values $\bm{\lambda}=(\lambda_{i})_{i\in
I^{-}}$ and $\mu$

Solve PM($\bm{\lambda}$, $\mu$) (or AHP($\bm{\lambda}$, $\mu$, $thr$,
$\bm{x}^{prev}$, $N^{max}$), see Section 4.1) and get:

(a) the number $N>0$ of columns generated

(b) vectors $\bm{x}^{k}\in\\{0,1\\}^{n}$ for each $k\in\\{1,\ldots,N\\}$

(c) scalars $b^{k}\in\mathbb{R}$ for each $k\in\\{1,\ldots,N\\}$

(d) vectors $\bm{w}^{k}\in\mathbb{R}^{d}$ for each $k\in\\{1,\ldots,N\\}$

$\bm{r}\leftarrow\left(-\sum_{i\in
I^{-}}x_{i}^{k}\cdot\lambda_{i}+\mu\right)_{k=1}^{N}$

$k^{min}\leftarrow argmin_{k=1}^{N}r_{k}$

$R\leftarrow r_{k^{min}}$

$\bm{x}^{prev}\leftarrow\bm{x}^{k^{min}}$

Add variable $z^{\bar{N}+k}$ and column
$\bm{c}^{\bar{N}+k}=\left(x^{k}\mathbin{\|}(-1)\right)^{\top}$ to RMM for each
$k\in\\{1,\ldots,N\\}$ such that $r_{k}<0$

Solve HCP and get the objective value, i.e., the error value $E$

Append $b^{k}$ to $\bm{b}$ for every $k\in\\{1,\ldots,N\\}$

Append $\bm{w}^{k}$ to $\bm{W}$ for every $k\in\\{1,\ldots,N\\}$

$\bar{N}\leftarrow\bar{N}+N$

end while

return $\\{(\bm{b}(k),\ \bm{W}(k))\ |\ k\in\\{1,\ldots,\bar{N}\\}\ \text{and}\
z^{k}=1\ \text{in the last HCP solution}\\}$

##  9 Adaptation to PCAB of a MIP from Literature

We show below the adaptation for PCAB of the model proposed in (Orsenigo and
Vercellis 2007) for convex polyhedral separation. We neglected part of its
objective function together with the respective variables and the related
constraints, in order to make the function equal to the separation error of
PCAB . We also added constraints (62) to impose correct separation of all the
positive points. The parameter $\epsilon$ is greater than $0$, while $Q$ is
equal to $|J|D$, where $D$ is greater or equal than the diameter of
$\mathcal{D}^{+}\cup\mathcal{D}^{-}$.

The binary variables take value 1 according to the following conditions:

  1. 1.

$e_{i}$ if and only if point $i$ is misclassified with respect to its true
label;

  2. 2.

$d_{i}^{k}$ if and only if point $i$ satisfies $\bm{w}^{k}\cdot\bm{a}_{i}\geq
w_{0}^{k}$, i.e. hyperplane $k$ classifies point $i$ as “positive”;

  3. 3.

$m_{i}$ if and only if point $i$ lies outside the convex region identified by
the $\mathcal{K}$ hyperplanes, i.e. $d_{i}^{k}=1$ for some $k\in K$;

  4. 4.

$v$ if and only if the points inside the polyhedron identified by the
$\mathcal{K}$ hyperplanes are labeled as “1” by the model;

The continuous variables represent the parameters of the hyperplanes. The
objective function of the model is a weighted sum of the separation errors
over all the points, positive and negative. Constraints (56) and (57)
identify, for each point $i$, which halfspace determined by hyperplane $k$ the
point is in. Constraints (58) ensures that $m_{i}=0$ if $d_{i}^{k}=0$ for
every $k\in K$, while constraints (59) impose that $m_{i}=1$ if $d_{i}^{k}=1$
for some $k\in K$. Constraints (60) and (61) characterize the error values
$e_{i}$ in terms of $m_{i}$, $v$ and the true labels of the points. Finally,
constraints (62) are not present in the original model (Orsenigo and Vercellis
2007) and impose that all the points whose true label is 1 are classified as
positive by all the hyperplanes. These constraints make the model suitable to
solve PCAB, which requires all the positive samples to be classified
correctly.

min $\displaystyle\sum_{i\in I}c(\bm{a}_{i})e_{i}$ (55) s.t. $\displaystyle
b^{k}+\sum_{j\in J}w_{j}^{k}a_{ij}\geq-Qd_{i}^{k}$ $\displaystyle\forall i\in
I,\forall k\in K$ (56) $\displaystyle b^{k}+\sum_{j\in
J}w_{j}^{k}a_{ij}\leq\left(1-d_{i}^{k}\right)Q-\epsilon$ $\displaystyle\forall
i\in I,\forall k\in K$ (57) $\displaystyle m_{i}\leq\sum_{k\in K}d_{i}^{k}$
$\displaystyle\forall i\in I$ (58) $\displaystyle|K|m_{i}\geq\sum_{k\in
K}d_{i}^{k}$ $\displaystyle\forall i\in I$ (59) $\displaystyle
2v-1-\ell(\bm{a}_{i})(1-2m_{i})\leq 2e_{i}$ $\displaystyle\forall i\in I$ (60)
$\displaystyle 2v-1-\ell(\bm{a}_{i})(1-2m_{i})\geq-2e_{i}$
$\displaystyle\forall i\in I$ (61) $\displaystyle d_{i}^{k}=0$
$\displaystyle\forall i\in I^{+},\forall k\in K$ (62) $\displaystyle
e_{i}\in\\{0,1\\}$ $\displaystyle\forall i\in I$ (63) $\displaystyle
d_{i}^{k}\in\\{0,1\\}$ $\displaystyle\forall i\in I,\forall k\in K$ (64)
$\displaystyle m_{i}\in\\{0,1\\}$ $\displaystyle\forall i\in I$ (65)
$\displaystyle v\in\\{0,1\\}$ (66) $\displaystyle b^{k}\in\mathbb{R}$
$\displaystyle\forall k\in K$ (67) $\displaystyle w_{j}^{k}\in\mathbb{R}$
$\displaystyle\forall j\in J,\forall k\in K$ (68)

##  10 Data Generation

$d$ |  N. det. positive points |  N. det. negative points |  N. random positive points |  N. random negative points |  Total n. positive points |  Total n. negative points
---|---|---|---|---|---|---
2 | 4 | 8 | 141 | 200 | 145 | 208
4 | 16 | 64 | 200 | 500 | 216 | 564
8 | 256 | 2048 | 282 | 8000 | 538 | 10048
Table 7: Details of the datasets for each $d\in\\{2,4,8\\}$.

For the first part of the experiments, we designed and generated what we call
“hypercubic data”, i.e. a dataset intended with the goal of obtaining
instances whose minimum separation budget is known _a priori_. We refer to
this dataset as D1 in the paper. The polyhedron to be approximated is the
$d$-dimensional hypercube $H_{d}=[0,1]^{d}$. After fixing a border gap
$\gamma=0.04$, for each vertex $v$ of the hypercube, we place:

  1. 1.

a positive point obtained from $v$ by replacing its 0 components with $\gamma$
and its 1 components with $1-\gamma$;

  2. 2.

$d$ negative points, obtained from $v$ by applying two subsequent
replacements. With the first one, we replace every 0 with $-\gamma$ and every
1 with $1+\gamma$. Then, in turn for each $j\in\\{1,2,\ldots,d\\}$, we replace
the $j^{th}$ element. If it is $-\gamma$, it becomes $2\gamma$ and if it is
$1+\gamma$ it becomes $1-2\gamma$.

The result of this generation for $d=2$ is shown in Fig. 5(a). The points are
placed in such a way that the minimum separation budget is $2d$. In order to
enlarge this initial set, random points are added, both positive and negative
(Fig. 5(b)), with the negative ones lying inside $[-1,2]^{d}\setminus H^{d}$.
We generated data in dimension 2, 4 and 8, producing 100 instances with the
same number of points of each type, for each of the three dimensions. The
number of deterministic, random and total points is shown in Table 7.

The dataset for the second part of the experiments, that we call D2, is
generated as D1, but without any border gap ($\gamma=0$) and without imposing
any corner points in the instances. In this way, the minimum separating budget
is possibly less that $2d$, while the number of facets of the convex hull of
the positive points can be greater than $2d$. The number of positive and
negative points in each dataset is the same as in the counterpart with corner
points excluding such points.

(a) Corner points

(b) Example of instance

Figure 5: Hypercubic data in dimension 2: points placed deterministically
(5(a)) and entire instance representation including random points (5(b)).

##  11 Additional Figures and Tables

Figures 6 \- 11 present the results of preliminary experiments on
representative instances with corner points, both in dimension 2 (Figures 6 \-
8) and in dimension 4 (Figures 9 \- 11). Such Figures are analogous to Figures
2 \- 4 presented in Section 5 for dimension 8. In dimension 2, the performance
difference between different methods is scarse. In dimension 4, COLGEN 1,
COLGEN 2 and MODEL B confirm to obtain the best results in terms of final
error value.

Tables 8 \- 10 show numerical results specific to our column generation
algorithms, COLGEN 1 and COLGEN 2. We report number of iterations, number of
columns generated and last plain pricing gap between objective and lower bound
(only for for COLGEN 1). We also show the amount of time each of the two
algorithms spent solving the pricing problem, the restricted master problem
and the hyperplane choice problem respectively. We observe that the number of
iterations that COLGEN 2 is able to complete can be significantly lower than
that reported by COLGEN 1 when $\mathcal{K}<2d$, especially with the datasets
in higher dimensions.

However, the parallelism inherent in COLGEN 2 enables generation of a number
of columns that is often greater with respect to the number of COLGEN 1. The
number of parallel processes is an essential hyperparameter in this sense:
extra preliminary experiments with a single process have shown far worse
results in terms of error values compared to our standard experiments. The
values of the last pricing gap of COLGEN 1 decrease when $\mathcal{K}$
increases, which indicates that finding a solution to the pricing model is
especially harder for the solver when the hyperplane budget is lower. Finally,
we can notice that the majority of the running times of the column generation
algorithms is spent in solving the pricing model, while the restricted master
model takes the smallest fraction of time and the hyperplane choice problem is
in between.

Figure 6: Error value with 2-dimensional hypercubic data: $\mathcal{K}=2$.
Figure 7: Error value with 2-dimensional hypercubic data: $\mathcal{K}=3$.
Figure 8: Error value with 2-dimensional hypercubic data: $\mathcal{K}=4$.

Figure 9: Error value with 4-dimensional hypercubic data: $\mathcal{K}=4$. Figure 10: Error value with 4-dimensional hypercubic data: $\mathcal{K}=7$. Figure 11: Error value with 4-dimensional hypercubic data: $\mathcal{K}=8$. $\mathcal{K}$ | method |  n. iterations |  n. columns |  last pricing gap |  PM time (seconds) |  RMM time (seconds) |  HCM time (seconds)
---|---|---|---|---|---|---|---
2 | COLGEN 1 | 71.69 | 71.69 | 74.53 | 5.797 | 0.078 | 0.130
COLGEN 2 | 65.09 | 520.72 | - | 5.289 | 0.165 | 0.574
3 | COLGEN 1 | 77.09 | 77.09 | 34.21 | 5.759 | 0.097 | 0.145
COLGEN 2 | 65.07 | 520.55 | - | 5.271 | 0.169 | 0.590
4 | COLGEN 1 | 33.91 | 33.91 | 0.46 | 1.242 | 0.057 | 0.073
COLGEN 2 | 5.35 | 42.79 | - | 0.437 | 0.012 | 0.019
Table 8: Column generation results for $d=2$: number of iterations, number of columns generated, last pricing gap and running times (pricing model, restricted master model, hyperplane choice model). $\mathcal{K}$ | method |  n. iterations |  n. columns |  last pricing gap |  PM time (seconds) |  RMM time (seconds) |  HCM time (seconds)
---|---|---|---|---|---|---|---
4 | COLGEN 1 | 334.37 | 334.37 | 114.322 | 53.185 | 1.099 | 5.633
COLGEN 2 | 110.55 | 884.40 | - | 15.109 | 0.818 | 44.449
7 | COLGEN 1 | 391.23 | 391.23 | 54.84 | 53.072 | 1.269 | 4.445
COLGEN 2 | 280.72 | 2245.76 | - | 36.563 | 2.780 | 20.506
8 | COLGEN 1 | 235.33 | 235.33 | 11.85 | 28.732 | 0.787 | 2.668
COLGEN 2 | 7.93 | 63.44 | - | 1.040 | 0.0448 | 0.082
Table 9: Column generation results for $d=4$: number of iterations, number of columns generated, last pricing gap and running times (pricing model, restricted master model, hyperplane choice model). $\mathcal{K}$ | method |  n. iterations |  n. columns |  last pricing gap |  PM time (seconds) |  RMM time (seconds) |  HCM time (seconds)
---|---|---|---|---|---|---|---
8 | COLGEN 1 | 210.88 | 210.88 | 2118.737 | 474.796 | 10.961 | 113.015
COLGEN 2 | 17.24 | 137.92 | - | 7.843 | 5.114 | 733.746
15 | COLGEN 1 | 182.15 | 182.15 | 1030.32 | 375.068 | 8.613 | 46.755
COLGEN 2 | 125.34 | 1002.72 | - | 54.187 | 19.545 | 528.832
16 | COLGEN 1 | 136.49 | 136.49 | 617.52 | 244.839 | 6.478 | 34.795
COLGEN 2 | 21.41 | 171.28 | - | 9.248 | 3.698 | 13.718
Table 10: Column generation results for $d=8$: number of iterations, number of
columns generated, last pricing gap and running times (pricing model,
restricted master model, hyperplane choice model).


11institutetext: Department of EEE, Imperial College London, London, SW7 2AZ,
United Kingdom.
<EMAIL_ADDRESS>
# A Note On Square-free Commuting Probabilities of Finite Rings

Andrew Mendelsohn

###### Abstract

It is shown that the commuting probability of a finite ring cannot be a
fraction with square-free denominator, resolving a conjecture of Buckley and
MacHale.

##  1 Introduction

Let $G$ be a finite group and $U$ denote the uniform distribution. Let the
commuting probability of $G$, $P_{G}$, be denoted

$P_{G}:=Pr_{a,b\leftarrow U(G)}(ab=ba).$

An alternative characterisation is

$P_{G}=\frac{1}{|G|^{2}}|\\{(x,y)\in G^{2}\text{ : }xy=yx\\}|.$

Joseph made the following conjectures, where $\mathcal{G}$ is the set of all
commuting probabilities of finite groups [1]:

  1. $\text{A}.$

All limit points of $\mathcal{G}$ are rational.

  2. $\text{B}.$

The set $\mathcal{G}$ is well ordered by $>$.

  3. $\text{C}.$

The set ${0}\cup\mathcal{G}$ is closed (that is, contains all its accumulation
points).

In [2], Eberhard resolved conjectures A and B, and in [3] conjecture C was
resolved.
One can extend the definition of commuting probability to finite rings: set

$P_{R}=\frac{1}{|R|^{2}}|\\{(x,y)\in R^{2}\text{ :
}xy=yx\\}|=\frac{1}{|R|^{2}}|\\{(x,y)\in R^{2}\text{ : }xy-yx=0\\}|.$

In [4], the following conjectures were made, where $\mathcal{R}$ is the set
$P_{R}$ over all finite rings $R$:

  1. 1.

$1/n\not\in\mathcal{R}$ when $n\in\mathbb{N}$ is square-free.

  2. 2.

$\mathcal{R}\subset\mathcal{G}$.

  3. 3.

$\mathcal{R}$ coincides with the set of values of $P_{G}$ as $G$ ranges over
all finite nilpotent groups of class at most 2.

  4. 4.

All limit points of $\mathcal{R}$ are rational.

  5. 5.

For each $0<t\leq 1$, there exists $\epsilon_{t}>0$ such that
$\mathcal{R}\cap(t-\epsilon_{t},t)=\emptyset$.

  6. 6.

$\mathcal{R}$ does not contain any of its accumulation points.

Note conjectures 4, 5, and 6 correspond to Joseph’s conjectures. Moreover,
conjectures 4 and 5 would follow from the veracity of conjectures 2 or 3,
since Eberhard showed $\mathcal{G}$ has rational limit points and is well-
ordered. In [5], conjecture 2 was in fact resolved, and thus conjectures 4 and
5. Moreover, conjecture 3 was partially resolved: the authors obtained that
$\mathcal{R}$ is a subset of the set of values of $P_{G}$ as $G$ ranges over
all finite nilpotent groups of class at most 2. We conclude that conjectures
1, 3, and 6 are open.
In this work, we resolve conjecture 1.

##  2 Preliminaries and Prior Results

######  Definition 1.

Two finite groups $G,H$ are called isoclinic if $G/Z(G)\cong H/Z(H)$ and
$G^{\prime}\cong H^{\prime}$, and if the diagram below commutes:

${{G}/{Z(G)}\times{G}/{Z(G)}}$${{H}/{Z(H)}\times{H}/{Z(H)}}$${G^{\prime}}$${H^{\prime}}$

Isoclinism preserves nilpotency class and commuting probability [10]. A stem
group is a group in a given isoclinism class of minimal order. It is well
known that if $G$ is a stem group, then $Z(G)\leq G^{\prime}$. For more on
isoclinism, see [11].
Below we state existing results in the literature we will need below.

######  Lemma 1.

[5] $\mathcal{R}\subset\mathcal{G}_{n,2}$, where $\mathcal{G}_{n,2}$ is the
set of commuting probabilities of all finite nilpotent groups of class at most
2.

This statement is proved as follows: let $R$ be a finite ring. We can turn
$R\oplus R$ into a nilpotent ring of class 3 by endowing it with the
multiplication rule $(a,x)(b,y)=(0,ab)$. This ring can be turned into a
nilpotent group $G_{R}$ of class at most 2 by endowing it with the binary
operation $a\circ b=a+b+ab$. Both of these transformations preserve the
commuting probability. Thus the values of $\mathcal{R}\setminus\\{1\\}$ are a
subset of the values $P_{G}$, running over nilpotent groups $G$ of class equal
to 2. Note that if $R$ has size $n$, then the resulting group $G_{R}$ has
order $n^{2}$, and if $R$ is noncommutative then the resulting group is
nonabelian.

######  Lemma 2.

[7]
$\operatorname{Pr}(G)=\frac{1}{\left|G^{\prime}\right|}\left(1+\frac{\left|G^{\prime}\right|-1}{|G:Z(G)|}\right)$
if and only if $G$ is nilpotent.

######  Lemma 3.

[6] If $G$ is a nilpotent group, then $P_{G}\neq\frac{1}{p}$.

##  3 Results

######  Theorem 3.1.

$\frac{1}{p}\not\in\mathcal{R}$ for all $p\in\mathbb{N}_{\geq 2}$.

###### Proof.

By Lemma 1, $\mathcal{R}$ is contained within the set of commuting
probabilities of finite nilpotent groups of class at most 2. By Lemma 3, this
latter set does not contain $\frac{1}{p}$ for any prime $p$. ∎

Denote the set of commuting probabilities of rings of prime power order for
some prime $p$ by $\mathcal{R}_{p}$.

######  Proposition 1.

$\frac{1}{n}\not\in\mathcal{R}_{p}$ for any prime $p\in\mathbb{N}_{\geq 2}$
and $n\in\mathbb{N}_{>1}$.

###### Proof.

By Lemma 1, we need only consider commuting probabilities of finite nilpotent
groups of class at most 2. By Lemma 2, we know a fortiori the commuting
probability of finite nilpotent groups of class at most 2 in terms of derived
subgroups and centers. Suppose that for some $n\in\mathbb{N}_{\geq 2}$ we have

$\frac{1}{\left|G^{\prime}\right|}\left(1+\frac{\left|G^{\prime}\right|-1}{|G:Z(G)|}\right)=\frac{1}{n}.$

By the construction of [5] considered above, wlog let $|G|=p^{e}$ for some
even positive integer $e$. Then $Z(G)=p^{f}$ and $G^{\prime}=p^{g}$ with
$0<g\leq f<e$ (since $G$ is at most class 2). Then

$\displaystyle\frac{1}{n}$
$\displaystyle=p^{-g}\left(1+\frac{p^{g}-1}{p^{e-f}}\right)=p^{-g}+\frac{p^{g}-1}{p^{e-f+g}}$
(1) $\displaystyle=\frac{p^{e-f}+p^{g}-1}{p^{e-f+g}}.$ (2)

For this to hold, some power $p^{h}$ must divide the numerator, with $h>0$;
but this cannot hold; for if so, then one must have $p^{e-f}+p^{g}-1=kp$, for
some integer $k\neq 0$. But then $-1\equiv 0\bmod p$, a contradiction. ∎

######  Theorem 3.2.

$\frac{\ell}{n}\not\in\mathcal{R}$ for any squarefree $n\in\mathbb{N}_{>1}$
and $\ell<n$ with $\operatorname{gcd}(\ell,n)=1$.

###### Proof.

Any finite ring can be turned into a nilpotent group of class at most 2, such
that the commuting probability of the ring is equal to the commuting
probability of the group. The construction (outlined above) turns a
commutative ring into a group of class 1, and a noncommutative ring into a
nonabelian group of class at most 2, therefore of class equal to 2. The order
of the group is the square of the order of the ring, so the Sylow subgroups of
the group have order at least the square of a prime. Since the group is
nilpotent, it can be written as a product of its Sylow subgroups, which are
all of class at most 2, and the commuting probability of the group is the
product of the commuting probabilities of its Sylow subgroups. Thus it remains
to analyse the equation

$\frac{\ell}{n}=\prod_{i=1}^{m}\frac{p_{i}^{e_{i}-f_{i}}+p_{i}^{g_{i}}-1}{p_{i}^{e_{i}-f_{i}+g_{i}}},$

for $m>1$ and the $p_{i}$ distinct, where the $e_{i},f_{i}$, and $g_{i}$ are
as before. Via isoclinism, we may replace $G_{R}$ by a class two nilpotent
(stem) group $G$ with identical commuting probability and minimal order. Thus
we may assume that $Z(G_{R})=G_{R}^{\prime}$ (note isoclinism preserves
nilpotency class and $G_{R}$ is class two), and moreover that none of the
Sylow subgroups are abelian. The above equality simplifies to

$\frac{\ell}{n}=\prod_{i=1}^{m}\frac{p_{i}^{e^{\prime}_{i}-f^{\prime}_{i}}+p_{i}^{f^{\prime}_{i}}-1}{p_{i}^{e^{\prime}_{i}}},$

where the exponents $e_{i}^{\prime},f_{i}^{\prime}$ correspond to the group
$G$. We now proceed by induction on the number of prime factors of $|G_{R}|$,
denoted $m$.
By Lemma 14 of [9], if $P_{G_{R}}=\frac{\ell}{n}$ in lowest terms, the prime
factors of $n$ are precisely the prime factors of $|G_{R}|$. If $m=1$, it is
known that $\frac{\ell}{n}\not\in\mathcal{R}_{q}$ for any prime $q$ and
square-free $n$ (in fact, we know this to hold for $m\leq 69$ by Theorem 9 of
[4]). Suppose the statement is true up to $n=k-1$, and consider the case
$n=k$; suppose, for a contradiction, that the commuting probability is equal
to $\frac{\ell}{n}$, for some square-free integer $n$ and $\ell\leq n$ with
$\operatorname{gcd}(\ell,n)=1$, and without loss of generality that $n$ has
prime factors equal to the set of $p_{i}$, $i=1,...,k$:

$\frac{\ell}{n}=\prod_{i=1}^{k}\frac{p_{i}^{e^{\prime}_{i}-f^{\prime}_{i}}+p_{i}^{f^{\prime}_{i}}-1}{p_{i}^{e^{\prime}_{i}}}.$

Rearrange for the following:

$\frac{\ell\cdot
p_{k}^{e_{k}^{\prime}}}{n\cdot(p_{k}^{e^{\prime}_{k}-f^{\prime}_{k}}+p_{k}^{f^{\prime}_{k}}-1)}=\prod_{i=1}^{k-1}\frac{p_{i}^{e^{\prime}_{i}-f^{\prime}_{i}}+p_{i}^{f^{\prime}_{i}}-1}{p_{i}^{e^{\prime}_{i}}}.$

Writing the left hand side in lowest terms, we have

$\frac{\ell\cdot
p_{k}^{e_{k}^{\prime\prime}}}{n^{\prime}\cdot(p_{k}^{e^{\prime}_{k}-f^{\prime}_{k}}+p_{k}^{f^{\prime}_{k}}-1)}=\prod_{i=1}^{k-1}\frac{p_{i}^{e^{\prime}_{i}-f^{\prime}_{i}}+p_{i}^{f^{\prime}_{i}}-1}{p_{i}^{e^{\prime}_{i}}},$

where $n^{\prime}$ is not divisible by $p_{k}$. We have a commuting
probability on the right hand side with $k-1$ prime factors; so by the
induction hypothesis, the denominator of the left hand side has no square
factors. But we also have
$p_{k}^{e^{\prime}_{k}-f^{\prime}_{k}}+p_{k}^{f^{\prime}_{k}}-1$ on the
denominator of the left hand side, which is not divisible by $p_{k}$, and by
Lemma 14 of [9] must have prime factors equal to the set of $p_{i}$, for
$i=1,...,k-1$; moreover, there can be no cancellation between these factors
and $\ell$, by assumption on $\ell$. But then for at least one index $j$,
$n^{\prime}\cdot(p_{k}^{e^{\prime}_{k}-f^{\prime}_{k}}+p_{k}^{f^{\prime}_{k}}-1)$
has a prime factor $p_{j}$ with multiplicity at least two, which is a
contradiction. ∎

######  Remark 1.

Since $\frac{1}{p}$ is an accumulation point of $\mathcal{R}_{p}$,
$\frac{1}{n}$ is an accumulation point of $\mathcal{R}$ for all $n$. The above
result thus means that many accumulation points of $\mathcal{R}$ are not
contained in $\mathcal{R}$. As well as resolving the first conjecture stated
at the beginning of this note, the result also makes progress on the sixth
conjecture stated.

## References

  * [1] Joseph, K. Several Conjectures on Commutativity in Algebraic Structures. The American Mathematical Monthly. 84, 550-551 (1977), https://doi.org/10.1080/00029890.1977.11994411.
  * [2] Eberhard, S. Commuting probabilities of finite groups. Bulletin Of The London Mathematical Society. 47, 796-808 (2015).
  * [3] Browning, T. Limit Points of Commuting Probabilities of Finite Groups. (arXiv, 2022), https://arxiv.org/abs/2201.09402.
  * [4] Buckley, S. & MacHale, D. Contrasting the commuting probabilities of groups and rings. (arXiv, 2014), https://archive.maths.nuim.ie/staff/sbuckley/Papers/bm_g-vs-r.pdf.
  * [5] Jurás̆, M. & Ursul, M. On commuting probabilities in finite groups and rings. Journal Of Algebra Combinatorics Discrete Structures And Applications. (2022), https://jacodesmath.com/index.php/jacodesmath/article/view/148.
  * [6] Castelaz, A. Commutativity Degree of Finite Groups. (Wake Forest University, 2010), https://books.google.co.uk/books?id=QYxBnQAACAAJ.
  * [7] Nath, R. & Das, A. On a lower bound of commutativity degree. Rendiconti Del Circolo Matematico Di Palermo. 59 pp. 137-142 (2010, 4).
  * [8] Buckley, S., MacHale, D. & Sh, A. Finite rings with many commuting pairs of elements. (2014), https://archive.maths.nuim.ie/staff/sbuckley/Papers/bms.pdf.
  * [9] Buckley, S. & MacHale, D. Groups with $\operatorname{Pr}$(G) = 1/3, https://archive.maths.nuim.ie/staff/sbuckley/Papers/bm_GpCP_1_3.pdf.
  * [10] Lescot, P. Isoclinism Classes and Commutativity Degrees of Finite Groups. Journal Of Algebra. 177, 847-869 (1995).
  * [11] Berkovich, Y. Groups of Prime Power Order, Volume 1. (De Gruyter, 2008), https://doi.org/10.1515/9783110208221.