Under consideration for publication in Theory and Practice of Logic Programming

1

arXiv:cs/0311002v1 [cs.PL] 4 Nov 2003

PROGRAMMING PEARL
Computing Convex Hulls with a Linear Solver
FLORENCE BENOY and ANDY KING
Computing Laboratory, University of Kent, UK.
email: {p.m.benoy, a.m.king}@kent.ac.uk

FRED MESNARD
Iremia, Université de La Réunion, France.
email: fred@univ-reunion.fr
submitted 16 October 2002; revised 7 August 2003; accepted 4 November 2003

Abstract
A programming tactic involving polyhedra is reported that has been widely applied in the
polyhedral analysis of (constraint) logic programs. The method enables the computations
of convex hulls that are required for polyhedral analysis to be coded with linear constraint
solving machinery that is available in many Prolog systems.
Keywords: convex hull, polyhedra, abstract interpretation, linear constraints.

1 Introduction
Polyhedra have been widely applied in program analysis (Cousot and Halbwachs 1978)
particularly for reasoning about logic and constraint logic programs. In this context
polyhedra have been used in binding-time analysis (Vanhoof and Bruynooghe 2001),
cdr-coded list analysis (Horspool 1990), argument-size analysis (Benoy and King 1996),
time-complexity analysis (King et al. 1997), high-precision groundness analysis (Codish et al. 2001),
type analysis (Sağlam and Gallagher 1997), termination checking (Codish and Taboch 1999)
and termination inference (Mesnard and Neumerkel 2001; Genaim and Codish 2001).
All these techniques use polyhedra to describe relevant properties of the program and manipulate polyhedra using operations that include projection, emptiness
checking, inclusion testing for polyhedra, intersection of polyhedra (meet) and the
convex hull (join). The classic approach to polyhedral analysis (Cousot and Halbwachs 1978)
uses two representations: (i) frames and rays and (ii) systems of (non-strict) linear inequalities and employs the Chernikova algorithm to convert between them
(Le Verge 1992). The rationale for this dual representation is that the convex hull
can be computed straightforwardly with frames and rays whereas intersection is
more simply computed over systems of linear inequalities. A simpler tactic that has
been widely adopted in the analysis of logic programs is to use only the linear inequality representation and compute the convex hull by adapting (Benoy and King 1996)
a relaxation technique proposed in (De Backer and Beringer 1993). The elegance

2

Florence Benoy, Andy King and Fred Mesnard

of this approach is that it enables the convex hull to be computed without recourse to a dual representation: the problem is recast as a projection problem that
can be subcontracted to standard linear constraint solving machinery with minimal coding effort. Moreover, the performance is acceptable for many applications.
In fact this technique has been widely applied in the analysis of logic programs
(Codish and Taboch 1999; Genaim and Codish 2001; King et al. 1997; Mesnard and Neumerkel 2001;
Sağlam and Gallagher 1997). The next section outlines the method and the following section, an example implementation. The final section presents the concluding
discussion.

2 Method
Consider two arbitrary polyhedra, P1 and P2 , represented in standard form:
~ 1}
P1 = {~x ∈ Qn | A1 ~x ≤ B

~ 2}
P2 = {~x ∈ Qn | A2 ~x ≤ B

~i
such that P1 6= ∅ and P2 6= ∅ so that the problem is non-trivial. Note that Ai ~x ≤ B
are non-strict and therefore P1 and P2 are both closed. The problem in essence is to
compute the smallest polyhedron that includes P1 and P2 . Interestingly, the convex
hull of P1 ∪ P2 is not necessarily closed as is illustrated in the following example.
Example 2.1
Consider the 2-dimensional polyhedra P1 and P2 defined by:



 



 
1
0
0 



1 −1
0 







−1
0
2 
 0


P1 = ~x ∈ Q2 
~
x
≤
P
=
~
x
∈
Q
~
x
≤
−1
1
0
2
 0
 1



1




0
−1
0
0 −1
−1

Observe that P1 = {h0, 1i} is a point whereas P2 = hx, yi ∈ Q2 | x = y ∧ 0 ≤ x is
a half-line. Note too that P1 and P2 are closed whereas the convex hull of P1 ∪ P2
excludes the points hx, yi ∈ Q2 | x > 0 ∧ y = x + 1 and hence is not closed (see
the diagram below).
y
y
y
✒
✒
✻
✻
✻
3
3
3
✒
✒
2
1
0

2
1

r
0

1

2

✲
3 x

P1 and P2

0

2
1

r
0

1

2

✲
3 x

convex hull of P1 ∪ P2

0

0

1

2

3

✲
x

closure of convex hull of P1 ∪ P2

Since the convex hull of P1 ∪ P2 is not necessarily closed, the convex hull cannot
always be represented by a system of non-strict linear inequalities; in order to
overcome this problem, the closure of the convex hull of P1 ∪ P2 is computed. The
starting point for our construction is the convex hull of P1 ∪ P2 that is given by:


~x = σ1 ~x1 + σ2 ~x2 ∧ σ1 + σ2 = 1 ∧ 0 ≤ σ1 ∧
PH = ~x ∈ Qn
~1
~ 2 ∧ 0 ≤ σ2
A1 ~x1 ≤ B
∧ A2 ~x2 ≤ B

Computing Convex Hulls with a Linear Solver

3

To avoid the non-linearity ~x = σ1 ~x1 +σ2 ~x2 the system can be reformulated (relaxed)
~ i to define:
by putting ~y1 = σ1 ~x1 and ~y2 = σ2 ~x2 so that ~x = ~y1 + ~y2 and Ai ~yi ≤ σi B


~x = ~y1 + ~y2
∧ σ1 + σ2 = 1
∧ 0 ≤ σ1 ∧
PCH = ~x ∈ Qn
~ 2 ∧ 0 ≤ σ2
~ 1 ∧ A2 ~y2 ≤ σ2 B
A1 ~y1 ≤ σ1 B
Observe that PH ⊆ PCH . Moreover, unlike PH , PCH is expressed in terms of a
system of linear inequalities. Note too that PCH is closed since the projection of a
system of non-strict linear inequalities is closed. In fact the following proposition
asserts that PCH coincides with the closure of the convex hull of P1 ∪ P2 .
Proposition 2.1
PCH is the closure of the convex hull of P1 and P2 .
The proof uses the concept of a recession cone. The recession cone of a polyhedron
P , denoted 0+P , is defined by: 0+ P = {~y ∈ Qn | ∀λ ≥ 0 . ∀~x ∈ P . ~x + λ~y ∈ P }.
The intuition is that 0+ P includes a vector ~y whenever P includes all the half-lines
in the direction of ~y that start in P .
Proof
~ i }. Theorem 19.6 of (Rockafellar 1970) states that
Suppose Pi = {~x ∈ Qn | Ai ~x ≤ B
the closure of the convex hull of P1 ∪ P2 is the set (0+P1 + P2 ) ∪ (P1 + 0+P2 ) ∪
(∪{σ1 P1 + σ2 P2 | σ1 + σ2 = 1 ∧ 0 < σ1 , σ2 }). Intuitively, 0+P1 + P2 is P2 extended
in the directions of half-lines contained within P1 . Let ~x ∈ Pi , then ~y ∈ 0+Pi if
~i for all λ ≥ 0 which holds if and only if Ai ~y ≤ ~0
and only if Ai (~x + λ~y ) ≤ B
(Rockafellar 1970)[pp 62]. Therefore 0+P1 + P2 = {~x ∈ Qn | ~x = ~y1 + ~y2 ∧ A1 ~y1 ≤
~ 2 } and similarly P1 + 0+P2 = {~x ∈ Qn | ~x = ~y1 + ~y2 ∧ A1 ~y1 ≤
~0 ∧ A2 ~y2 ≤ B
~
B1 ∧ A2 ~y2 ≤ ~0}. Furthermore, ∪{σ1 P1 + σ2 P2 | σ1 + σ2 = 1 ∧ 0 < σ1 , σ2 } =
~ 2 }.
~ 1 ∧ A2 ~y2 ≤ σ2 B
{~x ∈ Qn | σ1 +σ2 = 1 ∧ 0 < σ1 , σ2 ∧ ~x = ~y1 + ~y2 ∧ A1 ~y1 ≤ σ1 B
n
~ 2 } coincides with
~ 1 ∧ A2 ~y2 ≤ σ2 B
Observe that {~x ∈ Q | ~x = ~y1 + ~y2 ∧ A1 ~y1 ≤ σ1 B
the sets (i) 0+P1 +P2 , (ii) P1+0+P2 and (iii) ∪{σ1 P1 +σ2 P2 | σ1+σ2 = 1 ∧ 0 < σ1 , σ2 }
when (i) σ1 = 0 and σ2 = 1, (ii) σ1 = 1 and σ2 = 0 and (iii) σ1+σ2 = 1 and 0 < σ1 , σ2
respectively. Therefore PCH is the closure of the convex hull.
This result leads to an algorithm for computing the closure of the convex hull:
~ i by scaling the constant vectors B
~ i by σi , add
construct the systems Ai ~
yi ≤ σi B
the constraints ~x = ~y1 + ~y2 , σ1 + σ2 = 1 and 0 ≤ σi , then eliminate variables
other than ~x using projection to obtain PCH in terms of ~x. Hence the closure of the
convex hull can be computed without recourse to another representation. This is
illustrated below.
Example 2.2
~ i:
Returning to example 2.1, consider the systems Ai ~x ≤ B




x − y ≤ 0∧

x
≤
0
∧
−x
≤
0
∧
P1 = hx, yi ∈ Q2
P2 = hx, yi ∈ Q2 −x + y ≤ 0 ∧


y ≤ 1 ∧ −y ≤ −1
−x ≤ 0

4

Florence Benoy, Andy King and Fred Mesnard

Adding ~x = ~y1 + ~y2 , σ1 + σ2 = 1 and

x = x1 + x2




0
≤ σ1

2
PCH = hx, yi ∈ Q
x1 ≤ 0



y1 ≤ σ1


x2 − y2 ≤ 0

0 ≤ σi leads to the following system:
∧
∧
∧
∧
∧

y = y1 + y2
0 ≤ σ2
−x1 ≤ 0
−y1 ≤ −σ1
−x2 + y2 ≤ 0

∧ σ1 + σ2 = 1
∧
∧
∧
∧ −x2 ≤ 0

Eliminating the variables xi , yi and σi leads to the solution:


∧ 










PCH = {hx, yi ∈ Q2 | 0 ≤ x ∧ x ≤ y ∧ y ≤ x + 1}
Theorem 19.6 of (Rockafellar 1970), which is used in the proof, asserts that PCH
includes P1 + 0+P2 = P1 + P2 = {hx, yi ∈ Q2 | x ≥ 0 ∧ y = x + 1} and therefore
includes the points {hx, yi ∈ Q2 | x > 0 ∧ y = x + 1}, and hence ensures closure.
Note that calculating PCH without the inequalities 0 ≤ σ1 and 0 ≤ σ2 – the
relaxation advocated in (De Backer and Beringer 1993) for computing convex hull
– gives {hx, yi ∈ Q2 | 0 ≤ x} which is incorrect.

3 Implementation
This section shows how closure of the convex hull can be implemented elegantly
using a linear solver in particular the CLP(Q) library (Holzbaur 1995). The behaviour of a predicate is described with the aid of modes, that is, + indicates an
argument that should be instantiated to a non-variable term when the predicate is
called; - indicates an argument that should be uninstantiated; and ? indicates an
argument that may or may not be instantiated (Deransart et al. 1996).
3.1 Closed Polyhedra
Closed polyhedra will be represented by lists (conjunctions) of linear constraints
of the form c ::= e ≤ e | e = e | e ≥ e where expressions take the form e ::= x |
n | n ∗ x | −e | e + e | e − e and n is a rational number and x is a variable. A
convenient representation for a closed polyhedron is a (non-ground) list of constraints. This representation is interpreted with respect to a totally ordered (finite)
set of variables. The ordering governs the mapping of each variable to its specific
dimension. In practise, the ordering on variables is itself represented by the position of each variable within a list. Specifically, if C is a list of linear constraints
[c1 , . . . , cm ] and X is a list of variables [x1 , . . . , xn ], then the represented polyhedron
is PC,X = {hy1 , . . . , yn i ∈ Qn | (∧ni=1 xi = yi ) |=Q (∧m
j=1 cj )}. Note that although
the order of variables in X is significant, the order of the constraints in C is not.
Finally, let vars(o) denote the set of variables occurring in the syntactic object o.
Example 3.1
The polyhedron P1 from example 2.2 can be represented by the lists C1 = [x = 0,
y = 1] and X = [x, y], that is, P1 = PC1 ,X . Moreover, P2 = PC2 ,X where
C2 = [x = y, x ≥ 0] or alternatively C2 = [y + z ≥ x, x ≥ y + 2 ∗ z, y ≥ 0, z ≥ 0].

Computing Convex Hulls with a Linear Solver

5

Hence the dimension of PC,X is defined by the length of the list X rather than the
number of variables in C.

3.2 Projection
Projection is central to computing the convex hull. The desire, therefore, is to
construct a predicate project(+Xs,+Cxs,-ProjectCxs) that is true when for a
given list of dimensions Xs and a given list of constraints Cxs, ProjectCxs is the
projection of Cxs onto Xs. The specification of such a predicate is given below.
preconditions:
• Xs is a closed list with distinct variables as elements,
• Cxs is a closed list of linear constraints,
• Cxs is satisfiable.
postconditions:
•
•
•
•

Xs is a closed list with distinct variables as elements,
ProjectCxs is a closed list of linear constraints,
vars(ProjectCxs) ⊆ vars(Xs),
PCxs,Xs = PProjectCxs,Xs .

Such a predicate can be constructed by adding the given constraints to the store
and then invoking the projection facility provided in the CLP(Q) library, that is,
the predicate dump(+Target, -NewVars, -CodedAnswer) (Holzbaur 1995). Quoting from the manual: “[dump] reflects the constraints on the target variables into
a term, where Target and NewVars are lists of variables of equal length and
CodedAnswer is the term representation of the projection of constraints onto the
target variables where the target variables are replaced by the corresponding variables from NewVars”. This leads to the following implementation of project:
:- use_module(library(clpq)).
project(Xs, Cxs, ProjectCxs) :tell_cs(Cxs),
dump(Xs, Vs, ProjectCxs), Xs = Vs.
tell_cs([]).
tell_cs([C|Cs]) :- {C}, tell_cs(Cs).
Example 3.2
For example, the query project([X, Z], [X < Y, Y < Z], ProjectCs) will correctly bind Cs to [X-Z<0]. However, correctness of this predicate is compromised by
existing constraints in the store. For instance, the compound query {X = Z + 1},
project([X, Z], [X < Y, Y < Z], ProjectCs) will fail because constraints
posted within tell cs interact with those already in the store.

6

Florence Benoy, Andy King and Fred Mesnard

To insulate the constraints posted in tell cs, both the variables Xs and the constraints Cxs need to be renamed. Renaming is trivial with the builtin copy term
but care must be taken to ensure that Xs and Cxs are renamed consistently, that
is that variable sharing in Xs and Cxs is preserved in the copies. However, in SICStus Prolog copy term(Term, Cpy) copies any constraints in the store that involve variables in Term. For example, the query {X=Y}, copy term(X=Y+1, Cpy)
will bind Cpy to A= B+1 where A and B are fresh variables. It will also copy
the constraint X = Y by posting the new constraint A = B to the store. To nullify this effect, copy term is called within the scope of call residue. The call
call residue(copy term(X=Y+1, Cpy), Residue) residuates any new constraint
into Residue instead of posting it to the store, thereby copying the term without
copying any constraint. Whether residuation is required depends on the particular
Prolog system. This leads to the following (SICStus Prolog specific) revision:
project(Xs, Cxs, ProjectCxs) :call_residue(copy_term(Xs-Cxs, CpyXs-CpyCxs), _),
tell_cs(CpyCxs),
dump(CpyXs, Vs, ProjectCxs), Xs = Vs.
Example 3.3
Using this revision, the query {X = Z + 1}, project([X, Z], [X < Y, Y < Z],
ProjectCs) will succeed binding ProjectCs to [X-Z<0]. However, adding Z = 5 to
the list of constraints induces an error. The problem is that posting the constraints
binds Z to 5 so that dump is called with its first argument instantiated to a list that
contains a non-variable term.
A pre-processing predicate prepare dump is therefore introduced to ensure that
dump is called correctly. The following revision to project, in effect, extends the
facility provided by dump to capture constraints over both uninstantiated and instantiated variables:
project(Xs, Cxs, ProjectCxs) :call_residue(copy_term(Xs-Cxs, CpyXs-CpyCxs), _),
tell_cs(CpyCxs),
prepare_dump(CpyXs, Xs, Zs, DumpCxs, ProjectCxs),
dump(Zs, Vs, DumpCxs), Xs = Vs.
prepare_dump([], [], [], Cs, Cs).
prepare_dump([X|Xs], YsIn, ZsOut, CsIn, CsOut) :(ground(X) ->
YsIn = [Y|Ys],
ZsOut = [_|Zs],
CsOut = [Y=X|Cs]
;
YsIn = [_|Ys],
ZsOut = [X|Zs],

Computing Convex Hulls with a Linear Solver

7

CsOut = Cs
),
prepare_dump(Xs, Ys, Zs, CsIn, Cs).
The literal prepare dump(+Xs, +Ys, -Zs, ?CsIn, -CsOut) is true for a given list
Xs which contains either variables or numbers (or a mixture of the two) and a given
list Ys which contains only variables, if
• Zs is the list obtained by substituting the non-variable terms of Xs with fresh
variables and
• CsOut is an open ended list of equality constraints with CsIn at its end
that contains one equality constraint for each number in Xs. Each constraint
equates a numeric element of Xs with the element of Ys that is in the same
list position.
The call prepare dump([X1, 1, X3, 2], [A, B, C, D], Zs, CsIn, CsOut), for
instance, will bind Zs to [X1, A,X3, B] and CsOut to [B=1,D=2|CsIn]. The predicate ensures that dump is called with its first argument bound to a list of free
variables even when the list Xs includes numbers. In the CLP(Q) library, numbers
coincide with rationals which are represented as compound (ground) terms of the
form rat(n, d) where n and d are integers. The ground(X) test effectively checks
whether X is instantiated to a number; the test number(X) is inappropriate since it
would always fail.
Example 3.4
Consider again example 3.1. The second representation of P2 can be simplified by
using projection as follows:
| ?- Cs = [Y+Z>=X,X>=Y+2*Z,Y>=0,Z>=0], project([X,Y], Cs, ProjectCs).
ProjectCs = [Y>=0,X=Y] ? ;
no
The system Cs is expressed over 3 variables and therefore defines a 3 dimensional
space. Intuitively, the projection onto [X, Y] is the shadow cast by PCs,[X,Y,Z] onto
the 2 dimensional space over X and Y. The projection ProjectCs in fact defines
a half-line confined to the first quadrant since, by rearranging Cs, it follows that
PCs,[X,Y,Z] = {hx, y, zi ∈ Q3 | x = y ∧ 0 ≤ y ∧ z = 0}.
3.3 Convex Hull
The specification for the main predicate convex hull(+Xs, +Cxs, +Ys, +Cys,
-Zs, -Czs), and then its code, is given below.
preconditions:
• Xs is a closed list with distinct variables as elements and likewise for Ys,
• Xs and Ys have the same length,
• vars(Xs) ∩ vars(Ys) = ∅,
• Cxs and Cys are closed lists of linear constraints,
• Cxs and Cys are both satisfiable,
• vars(Cxs) ⊆ vars(Xs) and vars(Cys) ⊆ vars(Ys).

8

Florence Benoy, Andy King and Fred Mesnard

postconditions:
•
•
•
•
•

Xs, Ys and Zs are closed lists with distinct variables as elements,
Zs is the same length as both Xs and Ys,
Czs is a closed list of linear constraints,
vars(Czs) ⊆ vars(Zs) and (vars(Xs) ∪ vars(Ys)) ∩ vars(Zs) = ∅,
PCzs,Zs is the closure of the convex hull of PCxs,Xs ∪ PCys,Ys .
convex_hull(Xs, Cxs, Ys, Cys, Zs, Czs) :scale(Cxs, Sig1, [], C1s),
scale(Cys, Sig2, C1s, C2s),
add_vect(Xs, Ys, Zs, C2s, C3s),
project(Zs, [Sig1 >= 0, Sig2 >= 0, Sig1+Sig2 = 1|C3s], Czs).
scale([], _, Cs, Cs).
scale([C1|C1s], Sig, C2s, C3s) :C1 =.. [RelOp, A1, B1],
C2 =.. [RelOp, A2, B2],
mul_exp(A1, Sig, A2),
mul_exp(B1, Sig, B2),
scale(C1s, Sig, [C2|C2s], C3s).
mul_exp(E1, Sigma, E2) :- once(mulexp(E1, Sigma, E2)).
mulexp( X,
mulexp(N*X,
mulexp( -X,
mulexp(A+B,
mulexp(A-B,
mulexp( N,

_,
X)
_,
N*X)
Sig,
-Y)
Sig,
C+D)
Sig,
C-D)
Sig, N*Sig)

::::::-

var(X).
ground(N), var(X).
mulexp(X, Sig, Y).
mulexp(A, Sig, C), mulexp(B, Sig, D).
mulexp(A, Sig, C), mulexp(B, Sig, D).
ground(N).

add_vect([], [], [], Cs, Cs).
add_vect([U|Us], [V|Vs], [W|Ws], C1s, C2s) :add_vect(Us, Vs, Ws, [W = U+V|C1s], C2s).
The predicate mulexp(?E1, ?Sigma, -E2) scales the numeric constants that occur within E1 by the variable Sigma, providing they are not coefficients of variables,
to obtain the expression E2. Note that Sigma is a variable and the expression E1 may
be a variable, hence both E1 and Sigma have mode ? rather than +. Since a nonground representation is employed for expressions, the test var(X) is used to determine whether the expression is a variable. As before, the test ground(N) detects numeric constants – rational numbers – which are the only type of subexpressions that
are ground. Observe that mulexp can return more than one solution, for example,
mulexp(X, Sig, E2) generates E2 = X; X = -( A), E2 = -( A); X = -(-( A)),
E2 = -(-( A)) etc as solutions. Thus the pruning operator once is applied within
mul exp(?E1, ?Sigma, -E2) to prevent erroneous solutions.

Computing Convex Hulls with a Linear Solver

9

The predicate scale(+C1s, ?Sigma, ?C2s, -C3s) scales each constraint within
the list C1s by the variable Sigma. Each constraint consists of a binary operator and
two expressions, and scaling is applied to the numeric constants in each expression as
specified by mul exp. For example, scale([X+2 >= 1+Y, Y = Z], Sigma, Tail,
ScaledCs) binds ScaledCs to [Y = Z, X+2*Sigma >= 1*Sigma+Y | Tail]. Note
that scale finesses the problem of putting Cxs and Cys into the standard form
~ i before applying scaling. In standard form, X+2 >= 1+Y is Y-X =< 1 but
Ai ~yi ≤ B
scaling constants on both sides of the relational operator preserves equivalence in
that X+2*Sig >= 1*Sig+Y is equivalent to Y-X =< 1*Sig. The use of a difference
list avoids an unnecessary call to append in the body of convex hull.
The predicate add vect(+Us, +Vs, -Ws, ?C1s, -C2s) operates on the lists Us
= [U1 , . . ., Un ] and Vs = [V1 , . . ., Vn ] which correspond to the vectors ~y1 and
~y2 (as introduced in section 2). The argument Ws is instantiated to another list
of variables [W1 , . . ., Wn ], which corresponds with ~x. The predicate creates the
system of equalities [W1 = U1 +V1 , . . ., Wn = Un +Vn ] corresponding to the system
~x = ~y1 + ~
y2 . The scaled constraints output by the two calls to scale are passed to
add vect via its accumulator and thereby combined with the system of equalities.
For example, the call add vect([X1,Y1], [X2, Y2], Ws, Tail, Cs) returns the
bindings Cs = [ A=Y1+Y2, B=X1+X2|Tail] and Ws = [ B, A].
The predicate convex hull(Xs, Cxs, Ys, Cys, Zs, Czs) takes, as input, two
lists of constraints (Cxs and Cys) and their corresponding lists of variables (Xs
and Ys) and produces as output a single list of constraints Czs over the variables Zs that represents the closure of the convex hull of the two input polyhedra. If Xs and Ys are not variable disjoint, then the pre-requisite can be satisfied by appropriately renaming variables. Specifically, the variables Xs and constraints Cxs can be renamed with copy term(Xs-Cxs, CpyXs-CpyCxs) and the call
convex hull(Xs, Cxs, Ys, Cys, Zs, Czs) replaced with convex hull(CpyXs,
CpyCxs, Ys, Cys, Zs, Czs). Since the integrity of the constraint store is preserved by project and since project is the only source of interaction with the
store, then it follows that convex hull also does not side-effect any existing constraints. The following is an illustrative example.
Example 3.5
Running this code on the data of Example 2.2 gives:
| ?- convex_hull([X1,Y1],[X1=0,Y1=1],[X2,Y2],[X2>=0,Y2=X2],V,S).
S = [_A>=0,_A-_B>=-1,_A-_B=<0],
V = [_A,_B] ? ;
no

4 Discussion
This section discusses the method proposed in the paper, comparing it with related
techniques. The Chernikova method is exponential in the worst-case (Le Verge 1992)
and the Fourier-Motzkin method, like all projection techniques over linear inequalities (Chandru et al. 2000), is also exponential. The exponential behaviour of both

10

Florence Benoy, Andy King and Fred Mesnard
y
1

-1❅

❅
❅
P2 ❅
❅ ✲

❅ P1
❅
❅
-1
y
1

-1❅

y

✻

1

x

y

✻
❅
❅
❅
❅ ✲

❅
❅
❅
-1

✻
❍❍
❏
✯x
✟
❍❍♣
❏ ❍♣ ❍
♣ ♣1 ✟
♣♣
♣✟
♣
♣
❍
♣
-1❍
✟
♣
Q
♣ ♣❍
2❏
♣♣♣♣
♣❍
❍❍
♣♣ ❍
♣ ♣ ♣ ♣ ♣ ♣✟❏
♣✟
❍♣❍
✟
♣ ♣ ✟✟
❍❍
♣♣
❍❍❏
✟✟
♣♣ 1
Q
❍
✟✟
✟
❍
❏❍
♣
✟
♣
❍
1 ❍
✟ -1 ❍❍ ♣ ♣ ♣
❥
❍
z
♣
❍
-1
1

1

x

✻
1❍
✟
✟ ✡❏ ❍
❍ ♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣♣
✯x
✟
✟
❍ ♣ ✟
❍✟
✡ ❏ ❍✟
♣1
♣
♣♣
-1❍
✟
♣
❏
♣ ♣ ❍✡❍
✟ ♣♣
♣♣
❍✟
✡
♣ ♣ ✟✟
♣♣
❍ ♣❏
❍
♣♣ ❍
❏
✡ ✟♣ ♣ ♣
✟ ♣♣
❏❍
✡
♣
✟
♣
❍
♣
❥
❍
✟ -1 ❍❍ ♣ ♣ ♣ ♣ ✟✟✟1 ❍
♣
z
♣
❍♣ ✟
-1

Fig. 1. (i) P1 and P2 , (ii) Q1 and Q2 , (iii) conv(P1 ∪ P2 ), (iv) conv(Q1 ∪ Q2 )
methods stems from the same source: the possibly exponential relationship between
the number of vertices and the number of half-spaces that define a polyhedron. In
fact the problem of calculating the closure of the convex hull of two polyhedra is also
exponential even for bounded polyhedra (polytopes). This can be demonstrated by
considering the so-called cross polytope in n-dimensions which is the polyhedron
with the vertex set {h±1, 0, . . . , 0i, h0, ±1, . . . , 0i, . . . , h0, 0, . . . , ±1i}. The cross polytope can be defined by no less than 2n inequalities yet can arise as the convex hull
of two polyhedra both of which can be defined with O(n) inequalities. Specifically
consider the n-dimensional polyhedra
Pn
P1 = {hx1 , . . . , xn i ∈ Qn | ( i=1 −xi ≤ 1) ∧ (∧nj=1 xj ≤ 0)}
Pn
P2 = {hx1 , . . . , xn i ∈ Qn | ( i=1 xi ≤ 1) ∧ (∧nj=1 − xj ≤ 0)}

Because P1 and P2 are polytopes, they can be expressed in terms of their vertices:
P1 = conv({h0, 0, . . . , 0i, h−1, 0, . . . , 0i, h0, −1, . . . , 0i, . . . , h0, 0, . . . , −1i})
P2 = conv({h0, 0, . . . , 0i, h1, 0, . . . , 0i, h0, 1, . . . , 0i, . . . , h0, 0, . . . , 1i})

Since h0, 0, . . . , 0i is convexly spanned by h1, 0, . . . , 0i and h−1, 0, . . . , 0i, it follows
that cl(conv(P1 ∪ P2 )) = conv(P1 ∪ P2 ) = conv({h±1, 0, . . . , 0i, h0, ±1, . . . , 0i, . . .,
h0, 0, . . . , ±1i}) which is the n-dimensional cross polytope. The 2 and 3 dimensional
cases are denoted in Figure 1 by (i) P1 and P2 and (ii) Q1 and Q2 respectively for
which the cross polytopes are a solid square and an octahedron. Hence the problem
of calculating the closure of the convex hull is intrinsically exponential irrespective
of the algorithm employed.
Example 4.1

Computing Convex Hulls with a Linear Solver

11

The following query illustrates how the hull algorithm yields an exponential number
of inequalities for the 4 dimensional case.

| ?- Xs = [X1, X2, X3, X4], Ys = [Y1, Y2, Y3, Y4],
Cxs = [-1 =< X1+X2+X3+X4, X1 =< 0, X2 =< 0, X3 =< 0, X4 =< 0],
Cys = [ Y1+Y2+Y3+Y4 =< 1, 0 =< Y1, 0 =< Y2, 0 =< Y3, 0 =< Y4],
convex_hull(Xs, Cxs, Ys, Cys, Zs, Czs),
Zs = [A, B, C, D].
Czs = [A-B+C+D>=-1,
A-B-C+D>=-1,
A-B+C-D>=-1,
A-B-C-D>=-1,

A+B-C-D=<1,
A+B+C-D=<1,
A+B-C+D=<1,
A+B+C+D=<1,

A+B+C+D>=-1,
A+B-C+D>=-1,
A+B+C-D>=-1,
A+B-C-D>=-1,

A-B-C-D=<1,
A-B+C-D=<1,
A-B-C+D=<1,
A-B+C+D=<1] ? ;

no
However, it would be wrong to conclude from these examples that the frame and ray
representation is preferable – inequalities are unavoidable since they are required
for other polyhedral operations.
Despite the scaling problems that are inherent to any convex hull algorithm, in
practise the technique proposed in this paper has been widely applied in logic programming (Codish and Taboch 1999; Genaim and Codish 2001; King et al. 1997;
Mesnard and Neumerkel 2001; Sağlam and Gallagher 1997), mostly to satisfaction.
For example, in the context of inferring termination conditions for logic programs
this method is feasible since it accounts for 42% of this first pass of the analysis and
the first pass itself constitutes only 23% of the total analysis time (Mesnard and Neumerkel 2001).
Whether the approach presented in this paper is applicable depends on the application context. When only standard domain operations are required and performance is not an issue, this method has much to commend it. However, when the
application has to additionally reason, say, about integral points (Ancourt 1991;
Quinton et al. 1997) or parameterised polyhedra (Loechner and Wilde 1997) then
specialised polyhedral libraries are required. Further, if performance is important,
then recourse should be made to a polyhedral library, since a state-of-the-art implementation employing the Chernikova algorithm (Bagnara et al. 2002), will outperform the approach presented here.
We have presented a Prolog program for computing convex hulls using linear
solver machinery. As Holzbaur’s library is also available for CIAO Prolog, ECLiPSe,
XSB and Yap Prolog, the technique can be easily adapted to these systems. The
method is a reasonable compromise between conciseness, clarity and efficiency and
variants of this program have now been widely deployed.
Acknowledgements Thanks are due to Mats Carlsson, Bart Demoen, Pat Hill,
Joachim Schimpf and Raimund Seidel and, of course, the anonymous referees.

12

Florence Benoy, Andy King and Fred Mesnard
References

Ancourt, C. 1991. Gén’eration Automatique de Code de Transfert pour Multiprocesseurs
à Mémoires Locales. Ph.D. thesis, Université Paris 6.
Bagnara, R., Ricci, E., Zaffanella, E., and Hill, P. M. 2002. Possibly Not Closed
Convex Polyhedra and the Parma Polyhedra Library. In Static Analysis Symposium,
M. V. Hermenegildo and G. Puebla, Eds. Lecture Notes in Computer Science, vol. 2477.
Springer-Verlag, 213–229. See also http://www.cs.unipr.it/ppl/.
Benoy, F. and King, A. 1996. Inferring Argument Size Relationships with CLP(R). In
Logic-based Program Synthesis and Transformation (Selected Papers), J. P. Gallagher,
Ed. Lecture Notes in Computer Science, vol. 1207. Springer-Verlag, 204–223.
Chandru, V., Lassez, C., and Lassez, J.-L. 2000. Qualitative Theorem Proving in Linear Constraints. In International Symposium on Artificial Intelligence and Mathematics.
Long version to appear in the Annals of Mathematics and Artificial Intelligence.
Codish, M., Genaim, S., Sondergaard, H., and Stuckey, P. 2001. Higher-Precision
Groundness Analysis. In International Conference on Logic Programming, P. Codognet,
Ed. Lecture Notes in Computer Science, vol. 2237. Springer-Verlag, 135–149.
Codish, M. and Taboch, C. 1999. A Semantic Basis for the Termination Analysis of
Logic Programs. The Journal of Logic Programming 41, 1, 103–123.
Cousot, P. and Halbwachs, N. 1978. Automatic Discovery of Linear Restraints among
Variables of a Program. In Principles of Programming Languages. ACM Press, 84–97.
De Backer, B. and Beringer, H. 1993. A CLP language handling disjunctions of linear
constraints. In International Conference on Logic Programming, D. S. Warren, Ed. MIT
Press, 550–563.
Deransart, P., Ed-Dbali, A., and Cervoni, L. 1996. Prolog: The Standard. SpringerVerlag.
Genaim, S. and Codish, M. 2001. Inferring Termination Conditions for Logic Programs
using Backwards Analysis. In International Conference on Logic for Programming,
Artificial Intelligence and Reasoning, R. Nieuwenhuis and A. Voronkov, Eds. Lecture
Notes in Artificial Intelligence, vol. 2250. Springer-Verlag, 681–690.
Holzbaur, C. 1995. OFAI clp(Q,R) Manual. Tech. Rep. TR-95-09, Austrian Research
Institute for Artificial Intelligence (ÖFAI), Schottengasse 3, A-1010 Vienna, Austria.
Horspool, N. 1990. Mode Analysis Techniques for Discovery of Lists in Prolog. In
Object Management, D. Tsichritzis, Ed. Centre Universitaire d’Informatique, University
of Geneva, 305–312.
King, A., Shen, K., and Benoy, F. 1997. Lower-bound Time-Complexity Analysis of
Logic Programs. In International Symposium on Logic Programming, J. Maluszynski,
Ed. MIT Press, 261–276.
Le Verge, H. 1992. A note on Chernikova’s algorithm. Tech. Rep. 635, IRISA, Campus
Universitaire de Beaulieu, Rennes, France.
Loechner, V. and Wilde, D. K. 1997. Parameterized Polyhedra and their Vertices. International Journal of Parallel Programming 25, 6, 525–549. See also
http://icps.u-strasbg.fr/polylib/.
Mesnard, F. and Neumerkel, U. 2001. Applying Static Analysis Techniques for Inferring Termination Conditions of Logic Programs. In Static Analysis Symposium,
P. Cousot, Ed. Lecture Notes in Computer Science, vol. 2126. Springer-Verlag, 93–110.
Quinton, P., Rajopadhye, S. V., and Risset, T. 1997. On Manipulating Z-Polyhedra
using a Canonical Representation. Parallel Processing Letters 7, 2, 181–194.
Rockafellar, R. T. 1970. Convex Analysis. Princeton University Press.
Sağlam, H. and Gallagher, J. P. 1997. Constrained Regular Approximation of Logic

Computing Convex Hulls with a Linear Solver

13

Programs. In Logic Programming Synthesis and Transformation (Selected Papers), N. E.
Fuchs, Ed. Springer-Verlag, 282–299.
Vanhoof, W. and Bruynooghe, M. 2001. Binding-Time Annotations Without BindingTime Analysis. In International Conference on Logic for Programming, Artificial Intelligence and Reasoning, R. Nieuwenhuis and A. Voronkov, Eds. Lecture Notes in Artificial
Intelligence, vol. 2250. Springer-Verlag, 707–722.

