Research Rutcor Report
Preprocessing of Unconstrained Quadratic Binary Optimization
Endre Borosa Peter L. Hammera Gabriel Tavaresa
RRR 102006, April, 2006
RUTCOR Rutgers Center for Operations Research Rutgers University 640 Bartholomew Road Piscataway, New Jersey 08854-8003 Telephone: Telefax: Email: 732-445-3804 732-445-5472 [e-mail]
http://rutcor.rutgers.edu/rrr
Rutgers University, 640 Bartholomew Road, Piscataway, NJ 088548003, USA. E-Mail: {boros,hammer,gtavares}@rutcor.rutgers.edu
a RUTCOR,
Rutcor Research Report
RRR 102006, April, 2006
Preprocessing of Unconstrained Quadratic Binary Optimization
Endre Boros
Peter L. Hammer
Gabriel Tavares
Abstract. We propose several ecient preprocessing techniques for Unconstrained Quadratic Binary Optimization (QUBO), including the direct use of enhanced versions of known basic techniques (e.g., implications derived from rst and second order derivatives and from roofduality), and the integrated use (e.g., by probing and by consensus) of the basic techniques. The application of these techniques is implemented using a natural network ow model of QUBO. The use of the proposed preprocessing techniques provides: (i) a lower bound of the minimum of the objective function, (ii) the values of some of the variables in some or every optimum, (iii) binary relations (equations, inequalities, or non-equalities) between the values of certain pairs of variables in some or every optimum, and (iv) the decomposition (if possible) of the original problem into several smaller pairwise independent QUBO problems. Extensive computational experience showing the eciency of the proposed techniques is presented. Substantial problem simplications, improving considerably the existing techniques, are demonstrated on benchmark problems as well as on randomly generated ones. In particular, it is shown that 100% data reduction is achieved using the proposed preprocessing techniques for MAX-CUT graphs derived from VLSI design, MAX-Clique in c-fat graphs derived from fault diagnosis, and minimum vertex cover problems in random planar graphs of up to 500 000 vertices.
Acknowledgements: The third author was partially nancially supported by the Portuguese FCT and by the FSE in the context of the III Quadro Comunitrio de Apoio. a
RRR 102006
Page i
Contents
1 Introduction 2 Denitions and Notations 2.1 Posiforms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Persistency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Implication network . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3 Basic Preprocessing Tools 3.1 First order derivatives . . . 3.2 Second order derivatives and 3.3 Roofduality . . . . . . . . . 3.3.1 Strong persistency . 3.3.2 Weak persistency . . 3.3.3 Decomposition . . . 1 2 3 4 5 8 8 8 9 9 11 13 14 14 17 18 21 21 22 23 24 26 27 27 27 31 32 33 34 35 36 44
........ coderivatives ........ ........ ........ ........
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
. . . . . .
4 Combining Basic Tools 4.1 Enhancing roofduality by probing . . . . . . . . . . . . . . . . . . . . . . . 4.2 Consensus . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Preprocessing Algorithm and Implementation 6 Test Problems 6.1 Benchmarks with prescribed density . . . . . . . . 6.2 Maximum cliques of graphs . . . . . . . . . . . . 6.3 Minimum vertex cover problems of planar graphs 6.4 Graphs for MAXCUT . . . . . . . . . . . . . . . 6.5 MAX2SAT formulas . . . . . . . . . . . . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
. . . . .
7 Computational Experiments 7.1 Test environment . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.2 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8 Application: Optimal vertex covers of planar graphs 8.1 Deriving minimum vertex cover from QUBOs optimum 8.2 Preprocessing . . . . . . . . . . . . . . . . . . . . . . . 8.3 Optimization . . . . . . . . . . . . . . . . . . . . . . . 8.4 Minimum vertex covers of very large planar graphs . . 9 Conclusions A Characteristics of Benchmarks
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
. . . .
Page ii B Preprocessing Statistics
RRR 102006 48
RRR 102006
Page 1
1
Introduction
Let n be a positive integer, B = {0, 1}, and let R denote the set of reals. A function f : Bn R is called a pseudoBoolean function. Every pseudoBoolean function f is known (see [41]) to have a unique multilinear polynomial representation, the degree of which is called the degree of f . In this paper we shall consider the problem of Quadratic Unconstrained Binary Optimization (QUBO) min f (x) , (1) n
xB
concerning the minimization of a quadratic pseudoBoolean function f given by
n
f (x1 , , xn ) = c0 +
j=1
ci xi +
1 i<j n
cij xi xj ,
(2)
where c0 , ci for i = 1, , n and cij for 1 i < j n are given reals. PseudoBoolean functions and the optimization problem (1) were introduced in [41], and have since been shown to be a common model for a wide variety of discrete optimization problems arising in VLSI design ([9, 17, 19, 24, 48, 52, 78]), statistical mechanics ([9, 80, 71]), reliability theory and statistics ([67, 74, 75]), manufacturing ([2, 7, 8, 27, 53]), economics and nance ([42, 45, 54, 55, 57]), operations research and management science ([30, 37, 73, 84, 85]), and in numerous algorithmic problems of graph theory ([29, 38, 41, 69, 70, 72, 79]), etc. Several exact algorithms (see e.g., [10, 14, 25, 80, 36, 44, 49, 62, 68, 86]) as well as numerous heuristic procedures (see e.g., [3, 11, 28, 33, 34, 35, 50, 51, 59, 60, 63, 64, 65, 66]) were proposed in the literature for the solution of (1), which is generally an NP-hard optimization problem. In this paper we propose a family of preprocessing techniques, based on several computationally ecient transformations of problem (1), and aimed at simplifying it and possibly decomposing it into smaller problems of the same type. More precisely, the proposed preprocessing involves a series of transformations of the quadratic pseudoBoolean objective function f , including: (i ) The xation of some of the variables at values which must hold at every optimum, and the enforcement of certain binary relations (e.g., equations, inequalities, or nonequalities) between the values of certain pairs of variables, which must hold in every optimum; (ii ) The xation of some of the variables and the enforcement of some binary relations between certain pairs of variables, which must hold in at least one optimal solution of QUBO; (iii ) The possible decomposition of the problem into several smaller QUBO problems involving pairwise disjoint subsets of the original variables.
Page 2
RRR 102006
As a result, we obtain a constant K and quadratic pseudoBoolean functions fr : BSr R, r = 1, , c, where the sets of indices Sr {1, , n}, r = 1, , c are pairwise disjoint, such that c
xBn
min f (x) = K +
r=1
yBSr
min fr (y).
The proposed method has several key ingredients. The rst ingredient is the representation of a quadratic pseudoBoolean function by means of a network (see [18, 20, 82]). This representation allows an ecient derivation via network ow computations of a lower bound on the minimum of the function and allows the identication of variables whose values in the optimum can be easily predicted, as well as of other information which makes the simplication of the problem possible. The second ingredient is the methodology primarily based on the theory of roofduality, introduced in [40], combined with the use of rst and second order derivatives and co-derivatives. This methodology provides an eective lower bound for the minimum of a quadratic pseudoBoolean function, as well as information about the values of a subset of the variables in the optimum (so called linear persistencies), and about binary relations which must hold between certain pairs of variables in the optimum (so called quadratic persistencies). An additional component of the methodology is the decomposition of a large problem into several smaller ones, which can be derived at a low extra cost by combining the conclusions of roofduality with those oered by the network representation of the problem (see Section 3.3.3). The third ingredient of the proposed method consists of the integration of conclusions obtained from subproblems associated to the original problem. This integration is realized by combining the conclusions of probing (i.e. multiple application of roofduality to naturally associated subproblems), and of Boolean consensus (i.e. the exhaustive expansion of the detected linear and quadratic persistencies). We present an extensive computational evaluation of the proposed preprocessing method in Section 7, using various benchmark sets and randomly generated test problems of various types, involving thousands of variables, as described in Section 6. Our experience shows that for dense problems the proposed preprocessing technique is less eective as the size of the problems grows. It is demonstrated on numerous publicly available test problems that for relatively sparse problems, including in particular certain classes of structured problems, the proposed preprocessing method achieves substantial reductions in size at a very low computational cost. For instance, applying the method to QUBO problems corresponding to vertex cover problems in planar graphs involving up to 500 000 nodes, lead to the optimal xation of 100% of the variables, i.e. to the exact solution of the problem, in every single instance considered (see Section 8).
2
Denitions and Notations
Let V = {1, , n}, and for a subset S V, denote by 1S BV its characteristic vector, i.e. 1S = 1 for j S and 1S = 0 for j S. j j
RRR 102006
Page 3
We shall consider quadratic pseudoBoolean functions in n binary variables, as in (2). We shall use x = (x1 , , xn ) to denote a binary vector, or the vector of variables on which the pseudoBoolean function depend. The complement of a binary variable xi is denoted def by xi = 1 xi . Variables and their complements are called literals, and we denote by def L = {x1 , , xn , x1 , , xn } the set of literals. In the sequel, the letters u, v and w refer to literals, while bold face letters x, y, etc., denote vectors. Let us note that each literal u L itself can be considered as a simple pseudoBoolean function, since for a particular binary assignment y BV to the variables, we have u (y) = yj B if u = xj and u (y) = 1 yj B if u = xj .
2.1
Posiforms
PseudoBoolean functions can also be represented as posiforms, i.e. multilinear polynomial expressions depending on all the literals (x) =
T L
aT
uT
u,
(3)
0 whenever T = . If {u, u} T for some u L, then uT u is identically where aT zero over Bn and, therefore we shall assume that aT = 0. Similarly to the case of polynomial expressions, the size of the largest subset T L for which aT = 0 is called the degree of the posiform , and is denoted by deg (). The posiform is called quadratic (linear, cubic, etc.) if deg () 2 (1, 3, etc.) Note that a quadratic pseudoBoolean function may also have representations by a posiform of degree higher than two. For instance, the quadratic function f (x1 , x2 , x3 ) = 1 x1 x2 x3 + x1 x2 + x1 x3 + x2 x3 can also be represented by the cubic posiform = x1 x2 x3 + x1 x2 x3 . In this paper we shall be concerned only by quadratic posiforms of quadratic pseudo Boolean functions. Let us note that any quadratic pseudoBoolean function can be transformed to a quadratic posiform in time linear in its size, where the size of a pseudoBoolean function e.g., as given by (2), or of a posiform e.g., as given in (3), is the number of its nonzero terms. To do this, in each quadratic term cij xi xj in (2) having cij < 0, i < j we replace xj by 1 xj , and use the resulting transformation cij xi xj = cij xi (1 xj ) = cij xi + (cij )xi xj . If we do this for all quadratic terms of (2) with negative coecients we nd
n
f (x) = c0 +
j=1
cj xj +
1 i<j n cij >0
cij xi xj +
1 i<j n cij <0
(cij )xi xj ,
where ci = ci +
i<j n cij <0
cij
Page 4
RRR 102006
for i = 1, , n. Using now the same substitution xj = 1 xj , we rewrite those linear terms cj xj for which cj < 0 as cj xj = cj + (cj )xj , yielding at the end the quadratic posiform
n n
f (x) = c0 +
j=1 c >0 j
cj xj +
j=1 c <0 j
(cj )xj +
1 i<j n cij >0
cij xi xj +
1 i<j n cij <0
(cij )xi xj ,
where c0 = c0 +
n
cj .
j=1 c <0 j
It is also easy to see that from a quadratic posiform we can also derive the unique polynomial form of the represented quadratic pseudoBoolean function in linear time of its size. Thus, we can assume without any loss of generality that the objective function of problem (1) is given as a quadratic posiform. To simplify notation, we shall assume in the sequel that the quadratic pseudoBoolean function f is given (also) as a posiform (x) = a0 +
uL
au u +
u,vL
auv uv,
(4)
where as before, L denotes the set of literals, and au 0 and auv 0 for all u, v L. In the sequel, the letters f and g will denote pseudoBoolean functions as well as their unique multi-linear polynomial expressions, while the Greek letters and will denote posiforms.
2.2
Persistency
The concept of persistency was introduced in [40], referring to assignments which are valid at some or all optimal solutions of an optimization problem. More precisely, we say that for some index j V and binary value B strong (respectively, weak ) persistency holds for problem (1) if xj = holds for all (respectively, for some) minimizing points of f . We shall equivalently say that the assignment xj = is strongly (respectively, weakly) persistent. It will also be very useful to consider the somewhat more specic denition of autarkies. For a binary vector x BV and subset S V let us denote by x [S] BS the partial assignment (or subvector) dened as x [S] = (xj |j S). For a subset S V, a partial assignment y BS , and a binary vector x BV , let us denote by z = x[S y] the vector obtained from x by switching its S-components to y, i.e. zj = yj for j S, and zj = xj for j S. Given a subset S V we say that a partial assignment y BS is a strong (respectively, weak ) autarky for f if the inequality f (x[S y]) < f (x) (respectively, f (x[S y]) f (x)) holds for all binary vectors x BV , x [S] = y.
RRR 102006
Page 5
It is easy to see that if a partial assignment y BS is a strong (respectively, weak) autarky for f , then the assignments xj = yj for all j S are strong (respectively, weak) persistencies for (1). Note also that while persistency is a property of individual assignments to variables, autarky is a property of partial assignments. For instance, if y BV is a minimizing point of f , then it is also a weak autarky for f (strong if y is a unique minimum of f ), though obviously not all partial assignments are autarkies of an arbitrary quadratic pseudoBoolean function. In this paper we also extend the notion of persistency for binary relations, which we call quadratic persistencies. Given two literals u, v L, we say that the inequality u v is strongly (respectively, weakly) persistent for problem (1), if u (y) v (y) holds for all (respectively, for some) optimal solutions y BV of (1).
2.3
Implication network
Following [16, 18, 20, 82] we associate to a quadratic posiform given in (3) a capacitated network G with node set N and with nonnegative capacities uv 0 for all u, v N, u = v. We shall assume for the sake of simplicity that all the possible ordered pairs of nodes represent arcs of the network, some of which may have capacity 0. The complete set of arcs will be sometimes denoted by A. This assumption will simplify the statement of several of the properties which we shall need later, without limiting generality, since the complexity of computations to be performed in such networks will depend anyway only on the number of those arcs which have a positive capacity. In order to dene this network, let us rst introduce an additional variable x0 for which we shall assume the constant assignment x0 = 1. Using this new variable (and the assumption that it is always equal to 1), we can rewrite (3) as (x) a0 = uL au x0 u + u,vL auv uv, obtaining a homogeneous quadratic multilinear polynomial expression of a0 with nonnegative coecients. To simplify notation, we introduce N = L {x0 , x0 } as the set of nodes of G , and write the above homogeneous form as (x) a0 =
u,vN,u=v
auv uv,
(5)
where we assume that the indices of the coecients are sets of two literals (and not ordered pairs of literals). For the purpose of network denition, we also assume that (5) involves all pairs u, v N, u = v, in which we have auv 0 for all terms with the possibility that auv = 0 for some of them. The capacities of the arcs in G are dened by 1 uv = vu = auv 2 for all u, v N, u = v.
Page 6 Example 1. Consider the quadratic posiform given by
RRR 102006
= 8 + 2x1 + 2x3 + 4x4 + 2x5 + 2x1 x3 + 2x1 x5 + x1 x2 + 2x1 x4 + x2 x3 + x2 x4 + x2 x5 + 2x3 x4 + 2x3 x5 + 2x4 x5 . Then, its homogeneous representation is + 8 = 2x0 x1 + 2x0 x3 + 4x0 x4 + 2x0 x5 + 2x1 x3 + 2x1 x5 + x1 x2 + 2x1 x4 + x2 x3 + x2 x4 + x2 x5 + 2x3 x4 + 2x3 x5 + 2x4 x5 . and the associated network G is shown in Figure 1.
x1
0.5 1
x1
0.5
x2
1 0.5
1
1 1
x2
1
0.5 1
0.5
x0
1
x3
2
0.5 1
x3
2
1
x0
0.5
1
0.5
x4
1 1 1 1
x4
1 1
x5
x5
Figure 1: The network G corresponding to the posiform of Example 1. We indicate only those arcs which have positive capacities.
Let us note that a network G associated in this way to a quadratic posiform has a special symmetry, namely that we have uv = vu for all pairs of literals u, v N, u = v. (6)
We shall call such a pair (G , ) an implication network. To justify this term, observe that if auv uv is a term of then either it contributes auv to the value of , or we must have v (and equivalently v u). Similarly dened graphs, called implication graphs were u associated to quadratic disjunctive normal forms in [4].
RRR 102006
Page 7
We remark here that given an arbitrary capacitated network G with capacity function , we can obtain an implication network (G, ) from it by dening uv = vu = for all pairs of literals u, v N, u = v. Let us also note that given an arbitrary capacitated network (G, ) on node set N, we can associate to it a quadratic posiform by dening (G,) =
vL
uv + vu , 2
x0 ,v v +
vL
v,x0 v +
u,vL,u=v
uv uv,
(7)
and observe that if = (G,) , then G is the implication network (G, ), the symmetrized version of (G, ) as dened above. Thus, we have a one-to-one correspondence between implication networks and quadratic posiforms. As we shall see later, there is a much deeper connection between these objects. Given a network G on node set N, we assume in the sequel for the sake of simplicity that A = A (G) denotes the set of arcs of the network G, and that it consists of all pairs u, v N, u = v. To be able to state further connections between quadratic posiforms and implication networks, we need to recall rst that in network optimization theory, a mapping : A R+ is called a feasible ow from source x0 N to sink x0 N, in a capacitated network G on node set N, with capacity function RA , if it satises the following constraints: + 0 (u, v) uv (v, w)
(v,w)A
for all arcs (u, v) A for all nodes v N \ {x0 , x0 } . (8)
(u, v) =
(u,v)A
The value of a feasible ow is dened as v () =
(x0 ,v)A
(x0 , v) .
(9)
We shall denote by v (G, ) the maximum value of a feasible ow in the network (G, ), i.e. v (G, ) = max {v () | satises (8)} . Finally, the residual network is dened by the capacity function
uv = uv (u, v) + (v, u)
(10)
for all (u, v) A. Let us note that if is a feasible ow, then uv as dened in (7), is a posiform.
0 for all (u, v) A, and hence (G, ) ,
Page 8
RRR 102006
3
3.1
Basic Preprocessing Tools
First order derivatives
The ith rst order partial derivative (i V) of a quadratic pseudoBoolean function f is given by
i1 n
i (x1 , , xn ) = ci +
j=1
cij xj +
j=i+1
cij xj .
It is easy to show that in every local optimum, Proposition 1. (i) If i (x) (ii) If i (x) 0 for all x BV , then xi = 1 is a weak persistency; 0 for all x BV , then xi = 0 is a weak persistency.
If i (x) is strictly negative, respectively strictly positive, the above implications represent strong persistencies. These simple relations have been already noticed in Hammer and Rudeanu [41], and represent an essential component of even the most recent work on preprocessing (see [6]). Hammer, Hansen and Simeone [40] have shown that those strong persistencies which can be obtained from the analysis of partial derivatives, can also be obtained by roofduality (a tool to be described in the next subsection). Moreover, roofduality is a stronger preprocessing technique, as shown in [40], where an example is presented displaying persistencies found by roofduality, but not following from the analysis of the signs of partial derivatives. In view of these results, the preprocessing algorithm to be described in Section 5, which will exploit heavily roofduality, will not explicitly include an analysis of the signs of partial derivatives, since the conclusions derivable from such an analysis will be automatically included among those provided by roofduality.
3.2
Second order derivatives and coderivatives
A natural generalization of the concept of the rst order derivative allows us to establish some persistent binary relations to hold among the values taken in the optimum by certain pairs of variables. Hammer and Hansen [39] called the linear function i (x) j (x) + (xi xj ) cij the (i, j)th second order derivative of f and denoted it by ij . It can be seen that the linear function i (x) + j (x) + cij (1 (xi + xj )) has a similarly important role as ij ; it will be called (i, j)th second order co-derivative and will be denoted by ij . With these notations, ij (x) = f (x [{i, j} (1, 0)]) and ij (x) = f (x [{i, j} (1, 1)]) f (x [{i, j} (0, 1)]) f (x [{i, j} (0, 0)]) ,
RRR 102006
Page 9
i.e. evaluate the eect of simultaneously changing the values of xi and xj , while keeping the values of the other variables unchanged. The following statement identies a series of weakly, respectively strongly, persistent binary (order) relations of the form xi xj = 0 (i.e. xi xj ), xi xj = 0 (i.e. xi xj ), xi xj = 0 (i.e. xi xj ), and xi xj = 0 (i.e. xi xj ), between certain pairs (xi , xj ) of variables. Proposition 2. Let f (x1 , , xn ) be a quadratic pseudoBoolean function, and let xi and xj be two of its variables. Then, (i) If (ii) If
ij ij
(x) (x)
0 for every x BV then xi xj = 0 is a weak persistency; 0 for every x BV then xi xj = 0 is a weak persistency; 0 for every x BV then xi xj = 0 is a weak persistency; 0 for every x BV then xi xj = 0 is a weak persistency.
(iii) If ij (x) (iv) If ij (x)
If in any of the implications above the left hand side inequality is strict, then the corresponding persistency is strong.
3.3
Roofduality
Let us recall rst that the roofdual bound proposed in [40], can be determined eciently by maximum ow computations in the implication network G ([16, 18, 20]). Proposition 3. For any quadratic posiform given by (4) and any feasible ow in the corresponding implication network G the equality = a0 + v () + (G , ) holds, where the right hand side is a quadratic posiform. Therefore, a0 + v (G, )
xBV
(11)
min (x) .
(12)
Incidentally, for any network (G, ) we have v (G, ) v (G, ), thus we get the best lower bound in (12) exactly for the implication network form of a quadratic posiform. In fact the bound in (12) was shown in [16, 18] to be the same as the roofdual value of , introduced in [40]. 3.3.1 Strong persistency
Let us observe next that if is a feasible ow in an implication network (G, ) then in view of conditions (6) and of the equality v () = v () ,
Page 10 the symmetrized ow dened by (u, v) = (v, u) = (u, v) + (v, u) for all u, v N, u = v 2
RRR 102006
is also a feasible ow in (G, ). This implies that in any implication network, among the maximum ows there always exists a symmetric one, for which = holds. Let us consider therefore a symmetric maximum ow in the implication network (G , ), and let = (G , ) . As we observed above, the corresponding implication network G is exactly the residual network of G corresponding to the ow . Let us then dene S N as the set of nodes v N to which there exists a directed x0 v path in G , each arc of which has a positive residual capacity (we assume that x0 S). Furthermore, let T = {v | v S}. Since is a maximum ow in G , we cannot have x0 S, and consequently T S = follows by the symmetry (6) of the implication network G . The following fact is wellknown in network optimization: Proposition 4. The set S is unique, and is independent of the choice of the maximum ow . It is in fact the intersection of the source sides of all minimum cuts of G . Let us note also that {u, v} S cannot hold for any quadratic term auv uv of with positive coecient auv , since otherwise we would have a positive capacity arc (u, v) from u S to v T by the denition of arc capacities in the implication network associated to , leading to a positive capacity path from x0 to x0 , in contradiction with the maximality of . Thus, it follows that the assignment which sets all literals in S to 1 (there exists such an assignment, since T S = ) makes all terms of which involve literals from S or T vanish. Introduce J = {j V | {xj , xj } S = } and let y BJ be the partial assignment for which u (y) = 1 holds for all u S (and consequently, v (y) = 0 for all v T ). Proposition 5. The partial assignment y BJ is a strong autarky for (and hence for ). Consequently, the assignments xj = yj for all j J are strongly persistent for problem (1). Proof. See Theorem 11 in [18]. In fact the set of variables xj , j J consists exactly of those involved in the so called master roof as dened in [40]. This discussion shows that as a byproduct of a maximum ow computation, the above approach determines J in additional computing time, which is linear in the number of nonzero terms of , i.e. linear in the size of .
RRR 102006 3.3.2 Weak persistency
Page 11
Let us consider now the directed graph G obtained from G by keeping only those arcs which have a positive residual capacity. Since we can assume that the maximum ow is symmetric, we shall have (v, u) A G whenever (u, v) A G . (13)
Let us compute the strong components of this directed graph (necessitating linear time in the number of arcs, i.e. linear time in the size of [83]), and denote these strong components by Ki , i = 1, , c. It is easy to see that the symmetry (13) implies the following Proposition 6. For every strong component Ki of G we have either {v | v Ki } = Ki or {v | v Ki } = Ki for some i = i. Proof. Follows readily by (13). Let us label now as K1 , K1 , K2 , K2 , , K , K , those strong components which satisfy (15), in such a way that there is no directed path in G from Ki to Ki for i = 1, , , and there is no directed path in G from x0 to Ki for i = 1, , . Since is a maximum ow, we cannot have a directed path from x0 to x0 in G, and hence the symmetry conditions (13) imply the existence of such a labeling. Let Ji = {j V | {xj , xj } Ki = } and let yi BJi be the partial assignment for which u (yi ) = 1 for all u Ki , for i = 1, , . Proposition 7. The partial assignment yi is an autarky for , for i = 1, , . Moreover, it is a strong autarky if there is a directed path in G from x0 , or if there is a directed arc between Ki and Ki . Consequently, the assignments xj = yij for all j Ji and i = 1, , are all persistent assignments for problem (1) . Proof. The claim is implied by the fact that the terms of including variables in Ji , i = 1, , vanish. Example 2. Consider the quadratic pseudoBoolean function of Example 1, but represented now by a dierent posiform (= ), = 3 + x1 x2 + 2x1 x3 + 2x1 x4 + x1 x5 + x1 x5 + x2 x3 + x2 x5 + x2 x4 + 2x3 x4 + 2x3 x5 + 2x4 x5 . (15) (14)
Page 12
RRR 102006
In this case the maximum ow is = 0, and we have = . The residual implication network G = G has two nontrivial strong components besides {x0 } and {x0 }, K1 = {x1 , x2 , x3 , x4 , x5 } and K1 = {x1 , x2 , x3 , x4 , x5 }, both of them of type (15) (see Figure 2). Then, by Proposition 7 the minimum value of is (0, 0, 1, 0, 1) = 3, which is also its roofdual value. In fact, there is an arc from K1 to K1 , and hence the assignment y = (0, 0, 1, 0, 1) is strongly persistent, implying that it is the unique minimum of , in this particular case.
x1
x1
0.5 1 0.5 1 0.5
0.5 1
x2
0.5
1
x2
0.5
x0
x3
1
0.5
0.5 1
x3
x0
x4
0.5 1 1
0.5
0.5 0.5
x4
1 1
x5
x5
Figure 2: The implication network G of Example 2. We indicated only those arcs which have positive capacities.
Let us note that if there is a directed arc (u, v) between Ki and Ki for some i = 1, , , then symmetry (13) implies that an arc (v, u) must also exist between Ki and Ki . It is this property that makes the persistency result of Proposition 7 to be strong, in this particular situation. In general, deciding if a partial assignment yi is a strong autarky for , for i = 1, , is a NPhard decision problem. This result can easily be established, since the outcome of this decision depends on the optimization of NPhard sub-problems, each one associated to a component of type (14). Example 3. Consider the quadratic posiform given by = 2x1 x2 + 4x1 x2 + 2x2 x3 + 2x2 x3 + 2x1 x3 + 4x1 x3 + 2x4 x5 + 2x4 x5 + 2x5 x1 ,
RRR 102006
Page 13
and the associated network G shown in Figure 3. The strong components of G are: K1 K1 K2 K2 K3 = = = = = {x0 } , {x0 } , {x4 , x5 } , {x4 , x5 } and {x1 , x2 , x3 , x1 , x2 , x3 } .
Let us note that there is no directed path from K2 to K2 , there is no directed path from x0 to K2 , and that K2 and K2 satisfy condition (15). Let J2 = {4, 5} and let y2 BJ2 be the partial assignment for which x4 = x5 = 1. By Proposition 7, x4 = x5 = 0 must hold in a minimizer of . Let us now show that these two assignments are not strongly presistent. In the rst place one can verify that x1 = 1 must hold in all minimizers of . Therefore, the only term connecting K3 to the other components vanishes. So, any solution satisfying the equation x4 x5 + x4 x5 = 0, including x4 = x5 = 1, must hold in a minimizer of .
1 1
x2
2 1
x3
1 2
x0
x4
1 1
x5
1
x1
1 2 2 1
x1
1
x5
1 1
x4
x0
x3
1 1
x2
Figure 3: The network G corresponding to the posiform of Example 3. We indicate only those arcs which have positive capacities.
3.3.3
Decomposition
Let us assume now that we have already xed all weakly and strongly persistent variables (e.g., by Proposition 7), and that the strong components of the residual posiform , Ki , i = 1, , c, are all of type (14). Clearly, in this case the symmetry conditions (13) imply that there are no arcs between dierent strong components, i.e. does not involve a quadratic term auv uv, auv > 0 for which u Ki , v Kj and i = j. Thus, denoting by i the posiform corresponding to the induced subnetwork Ki , i = 1, , c, we have
c
=
i=1
i .
Page 14
RRR 102006
Furthermore, due to property (14), these posiforms involve disjoint sets of variables. Hence, we have Proposition 8.
xBV c
min (x) =
i=1
xBJi
min i (x) .
(16)
A similar decomposition of quadratic posiforms which does not involve linear terms was already proposed in [13]. Let us note that after computing the maximum ow in the implication network G , both the persistent assignments, as well as the above decomposition can be determined in linear time of the size of . Example 4. Consider to be an homogeneous quadratic posiform whose nonzero terms are: x1 x2 , x1 x2 , x2 x3 , x2 x3 , x1 x3 , x1 x3 , x4 x5 , x4 x5 , x5 x6 , x5 x6 , x4 x6 , x4 x6 , x7 x8 , x7 x8 , x1 x7 , x4 x8 , x9 x10 , x9 x10 , x7 x9 , x11 , x11 x12 , x11 x12 , x10 x11 . The associated network G , shown in Figure 4, has the following strong components: K1 K2 K3 K4 K5 K6 = = = = = = {x0 } , {x11 , x12 } , {x9 , x10 } , {x7 , x8 } , {x1 , x2 , x3 , x1 , x2 , x3 } and {x4 , x5 , x6 , x4 , x5 , x6 } . K1 K2 K3 K4 = = = = {x0 } , {x11 , x12 } , {x9 , x12 } , {x7 , x8 } ,
Let us rst note that there is no directed path from x0 to x0 . Thus, by strong persistency (Proposition 5) x11 = x12 = 1 must hold for all minimizers of . Also, regardless of the values of the coecients in the nontrivial terms of , by weak persistency (Proposition 7) the partial assignment x7 = x8 = x9 = x10 = x11 = x12 = 1 must hold in a minimizer of . This partial assignment automatically cancels all those terms of which involve at least a variable from the set {x7 , x8 , x9 , x10 , x11 , x12 }. After eliminating these terms, the original problem is decomposed into two subproblems, involving disjoint sets of variables, coming respectively from K5 and K6 . Obviously, these two subproblems can be optimized separately, and the sum of their optimal values will coincide with the minimum of .
4
4.1
Combining Basic Tools
Enhancing roofduality by probing
It will be seen in this section that the results of roofduality can be substantially strengthened by analyzing the roofduals of several quadratic pseudoBoolean functions associated to the original one. The analysis of the roofduals of the 2n quadratic pseudoBoolean functions
RRR 102006
x2 x3
Page 15
x1
x1
x11 x0 x9
x12
x7
x3
x2
x7
x12
x11 x0
x10
x8
x6
x5
x8
x10
x9
x4
x4
x5
x6
Figure 4: The network G corresponding to the posiform of Example 4. We disregarded the values of the capacities, and indicated only those arcs which have positive capacities. The dashed arcs represent arcs connecting the strong components of G .
obtained by xing either to zero or to one the value of one of the n variables of a quadratic pseudoBoolean function f (x1 , , xn ), will be called probing. Among the specic results of probing we mention the possibility of identifying an improved lower bound of the minimum of a quadratic pseudoBoolean function, and of enlarging at the same time the set of variables and binary relations for which persistency conclusions apply. In view of the fact that nding the roofdual of a quadratic pseudoBoolean function can be achieved by simply solving a maximum ow problem in a network, the calculation of the roofduals of 2n quadratic pseudoBoolean functions associated to the initially considered one is a feasible, easytocarryout operation. Let us provide now some technical ideas on how to eciently calculate the 2n roofduals required by probing. Let us assume rst that using a simple heuristic (e.g. as proposed in [21, 22]) we have
Page 16 found some upper bound U on the minimum of (1), say
xBV
RRR 102006
min (x)
U.
(17)
Let us now consider the most typical branching procedure, in which we split the problem into two somewhat smaller ones by xing variable xj at the two possible binary values. Since is a posiform, all of its terms with the possible exception of the constant a0 contribute nonnegative quantities to the objective. Therefore, if M > U a0 , then
xBV
min (x) = min
xBV
min (x) + Mxj , min (x) + Mxj ,
xBV
(18)
where a0 is the constant in , as given in (4). The two subproblems in (18) have simple network representations. In order to calculate the roofduals of minxBV (x) + Mxj and of minxBV (x) + Mxj , and to derive persistency relations from these, we should add to the original network an arc, (x0 , xj ) in the rst case and (x0 , xj ) in the second case, and to assign to these the large capacity M. Clearly, computationally it is simpler to increase the capacity of two arcs than to substitute xj = 0, respectively xj = 1, implying the deletion of nodes xj and xj , and of all arcs incident to these nodes in the network. In addition, keeping the same network and updating the capacities of a few arcs at each branch evaluation, allows us to carry out computations without increasing the amount of computer memory needed to keep the network data necessary to nd a maximum ow for each subproblem. From an implementational point of view, this approach has the added advantage of allowing the easy restoration of the network corresponding to the original problem, by simply nding an additional maximum ow an option which turned out to be on the average to be much better than creating a copy to be reused after branching. It should be remarked that without these simplifying steps, the large scale QUBOs (including for instance those coming from nding optimal vertex covers of planar graphs with half a million vertices; see Section 8) could not have been handled. We can similarly administer more complicated branching policies, as well. For instance, if u, v L are two literals, u = v, then branching on the binary relation u v can be written as
xBV
min (x) = min
xBV
min (x) + Mu + Mv, min (x) + Muv
xBV
(19)
for some M > U a0 , resulting in the modication of 4 arc capacities in the rst branch corresponding to u = 1 and v = 0, and of two arc capacities on the other branch corresponding to u v. u We can also apply the above for handling persistencies. For instance, if we learn that v is a persistent binary relation, then we can rewrite (19) as
xBV
min (x) = min (x) + Muv.
xBV
(20)
RRR 102006
Page 17
Let us note that in all of the above cases, we had to increase the capacity of some of the arcs. Thus, as our procedure advances, and we learn more and more persistencies, at the same time the maximum ow value is also increasing. Hence, according to (12), as an added value we get better and better lower bounds on the minimum of (1). To describe probing and its benets, let us consider an arbitrary quadratic posiform , as given in (4). For a literal u L let us consider the posiform u = + Mu, where M > U a0 for an upper bound U satisfying (17). Let us further denote by Su L (Wu L) the set of strongly (weakly) persistent literals for u , as dened in Section 3.3.1 (3.3.2), and let Lu denote the roofdual bound for u . We can derive several consequences from the sets Su , Wu and lower bounds Lu when generating these for all literals u L. Proposition 9. Let U be an upper bound of minxBV (x), and let u L and j V. Then The value L = maxuL min {Lu , Lu } is a lower bound on the minimum of . If Lu > U then u = 0 is a strongly persistent assignment for . If v Sxj Sxj (v Wxj Wxj ) then v = 1 is a strongly (weakly) persistent assignment for . If v Sxj and v Sxj (v Wxj and v Wxj ) then xj = v is a strongly (weakly) persistent relation for . For all v Sxj and w Sxj (v Wxj and w Wxj ) the quadratic relations xj xj w and w v are all strongly (weakly) persistent for . v,
All these follow from the above denitions, by which the assignment v = 1 is strongly (weakly) persistent for u , for all v Su (v Wu ). Let us note that by adding these new persistencies to , as in (18) and (20), we may increase both the roofdual value as well as the set of strongly and weakly persistent literals of . Furthermore, the addition of these to may also change the sets Sv or Wv for some other literals v L, v = u. Let us remark that the lower bound derived from probing was also considered in [16, 20], and that analogous techniques were explored in the broader context of binary optimization in [5, 77].
4.2
Consensus
It has been remarked above that order relations between literals can be derived both from the signs of the second order derivatives and coderivatives, as well as during the process of probing. The interactions of the various binary relations, and the conclusions obtained
Page 18
RRR 102006
by combining them can be very easily derived by the introduction of a Boolean quadratic equation = 0, where the terms of are quadratic elementary conjunctions (representing the detected binary relations between literals) which must take the value zero in every minimum of (1). Moreover, the application of the consensus method to allows the polynomial detection of all of its prime implicants (see [26]). Taking into account that the prime implicants of this function represent either variables with xed values in the optimum, or order relations which must hold between pairs of literals in every optimum, it is clear that the detection of all prime implicants of provides an enlarged set of strong persistencies. Finally, we should remark that the conclusions that can be obtained from the knowledge of the prime implicants of can also be obtained directly (using Proposition 9) during our implementation of probing, by appropriately transforming the original functions via term additions (as explained in the previous section) corresponding to the persistent binary relations found by the preprocessing tools considered.
5
Preprocessing Algorithm and Implementation
The proposed preprocessing algorithm is presented in this section, in which the tools described in the previous sections are used iteratively and recursively. The structure adopted is based on the network ow model of Section 2.3. Our goal is to nd better and better lower bounds, weak (and strong) linear persistencies and weak (and strong) quadratic persistencies for , as well as to decompose the problem into several smaller problems, as explained in Section 3.3.3. The PrePro algorithm is described in Figure 5. The input of the algorithm is a quadratic posiform representing a quadratic pseudoBoolean function f , as in (4). The output returned by PrePro is a decomposed representation of the minimum of f , as in (16), where the subproblems on the right hand side of (16) involve pairwise disjoint sets of variables, together with a lower and an upper bound to the minimum of each subproblem, and with a partial (optimal) assignment to the derived persistent variables. Four main components are part of PrePro: Network This routine has a posiform as input. It rst nds a maximum ow in the network G as explained in Proposition 11. The maximum ow implementation that we have adopted is based on the shortest augmenting path algorithm, yielding a worst case time of O (n3 ), and is specially designed to obtain a residual network satisfying the symmetry conditions (6). When a minimum cut is found, a set of strong persistencies is obtained directly from a nonempty source side of a minimum cut (see Proposition 4) and saved accordingly in (V0 , V1 ). The nodes of the residual network and corresponding arcs, which are associated to the set of strong persistencies is removed from the network. Using a linear time algorithm for identifying the strongly connected components of the residual network, a subset of weak persistencies is identied in every component of type (15) (see also Proposition 7), and saved accordingly in (V0 , V1 ). The nodes of the residual network and corresponding arcs, which are associated to the
RRR 102006
PrePro(, (V0 , V1 , E) , ) Input: Step 1: Step 2: Step 3: Step 4: Step 5: Step 6: Step 7:
Page 19
A quadratic posiform representing a quadratic pseudoBoolean function f , as in (4). Call Network(, , (V0 , V1 , E)). If = then and STOP. Otherwise, for all (, L, U ) do call Heuristic(, U ). For all (, L, U ) do call Coordination(, U, E). P {(, L, U ) | is purely quadratic }. For all (, L, U ) P do call Probing(, L, U, (V0 , V1 , E)). P {(, L, U ) P | is purely quadratic}. For all (, L, U ) ( \ P) do call PrePro (, (V0 , V1 , E) , S). P S.
Output: It returns in a set of triplets, whose sum of all posiforms has the same minimum value as the minimum value of f . Additionally, a set of persistencies for the initial posiforms is returned as follows: V0 A set of variables with value 0, V1 A set of variables with value 1, and E A set of equality relations between pairs of literals, which must hold in a minimizer of f .
Figure 5: PrePro algorithm.
set of weak persistencies are removed from the network. What is left after applying Network is a disjoint set of strong components, each corresponding to a subproblem (included in ) that can be optimized separately from the other subproblems. Heuristic For each subproblem identied in Network, an upper bound U is found, and later used in the Coordination and Probing procedures. Any heuristic method valid for QUBO can be used. All of our experimental results include a fast one pass heuristic based on a greedy approach, which ranks the set of nonxed variables according to an approximated probability value of a variable to have a persistent value (some onepass variants are proposed in [33]). Coordination This procedure has as input a posiform , and the upper bound U found by Heuristic. Let us note that the minimum of any posiform called within PrePro is strictly larger than its constant value a0 . This routine identies binary persistent relations originated from the analysis of the second order derivative and co-derivative functions as explained in Proposition 2. The basic idea is to compute over all possible pair of variables i < j, the minimum and the maximum of the linear pseudoBoolean functions ij and ij . A key element to save computing time in this operation is to stop it as soon as one learns that the minimum is strictly negative and the maximum is strictly positive. If a quadratic persistency u v is found,
Page 20
RRR 102006
then is updated by adding a term auv uv with a large enough coecient (we use auv = 2(U a0 ) + 1). Since the implication network structure was adopted in all tools, this last operation can be eciently implemented by updating the coecient of the arcs (u, v) and (v, u). Our data structure is also able to quickly identifying if the reverse relation v u is a quadratic persistency. In the armative case, E is updated to include the equality persistency u = v, and is transformed by replacing v (v) by u (u). This routine stops as soon as a linear term or a equality persistency is found. Probing This procedure has as input a purely quadratic posiform (i.e. a posiform that does not have linear terms) and the upper bound U found by Heuristic. The implication network structure plays a crucial role in this routine. Independently, each variable xj is xed to one of the possible binary values, and the resulting function is analyzed in terms of roofduality and strong and weak persistencies. This operation can be accommodated easily in the network as explained in Section 4.1. For a given assignment to xj , a maximum ow algorithm is applied to the transformed network. All the strong and weak persistencies derived from the residual network are updated both in the network and in (V0 , V1 , E), as explained in Proposition 9. To analyze the complement assignment of xj , a maximum ow algorithm is applied to the residual network. All the strong and weak persistencies derived from the residual network are again updated as before. A third maximum ow algorithm is applied to obtain a network which represents the same function as . The use of maximum ow algorithms to recuperate the original function being optimized is possible due to Proposition 3. A clear advantage of this approach is that the amount of memory needed for the data structures remains about the same through every step of the procedure. Let us note that probing through the implication network is able to capture persistencies of transitive relations (see Section 4.2). For instance, suppose that u v and v w are quadratic persistencies, then if at some point w u is also found to be persistent, then the network and the set E are immediately updated with the equality relation u = w. The routine stops as soon as a linear term or a linear/equality persistency is found. Step 4 of PrePro selects the subproblems for which probing is applied. Step 6 of PrePro selects the subproblems for which PrePro is recursively applied. Obviously, the rules that govern these choices may vary. In our implementation, the rule adopted is to apply Network to a subproblem whenever a new linear persistency or a new linear term was found by Coordination or probing. All the tools considered in PrePro are polynomial time algorithms: Network O (n3 ) ; Heuristic O (n2 ) ; Coordination O (n3 log (n)) ; Probing O (n4 ) . As a consequence of the previous complexity times, each run from Step 1 to Step 5 of PrePro takes at most O (n4 ).
RRR 102006
Page 21
6
Test Problems
Most of the problem classes on which we have tested the proposed preprocessing procedures are benchmarks used in several other studies related to QUBO. These datasets include 85 problems with prescribed density, 13 graphs for MAX-Clique, 38 graphs for MAXCUT, and 34 MAX2SAT formulas. Beside the above classes we have also carried out computational experiments on 436 randomly generated planar graphs for vertex cover optimization problems.
6.1
Benchmarks with prescribed density
The class of benchmarks with prescribed density consists of quadratic multilinear polynomials, which are randomly generated with a predetermined expected density1 d, a discrete uniform distribution of values of the linear coecients in an interval [c , c+ ], and a discrete uniform distribution of values of the the quadratic coecients in an interval [2q , 2q + ]. The constant term of the function is zero. Since a quadratic terms probability to have a nonzero coecient is d, the expected number of quadratic terms which include a specic variable and have a nonzero coecient is (n 1) d. This group of problems includes the test problems of Glover, Kochenberger and Alidaee [34], and the problems of Beasley [11] with at most 500 variables. The basic generation parameters of the subfamilies containing these problems can be seen in Table 1, while the individual characteristics of the problems appear in Tables A.1 and A.2 of the Appendix. Obviously, Glovers and Beasleys maximization problems have been rst converted to minimization problems in order to make the direct application of our proposed procedures possible.
Table 1: Benchmarks with prescribed densities for QUBO.
Family GKA Sub-Family A B C D E F1 ORL-50 ORL-100 ORL-250 ORL-500 Number of Problems 8 10 7 10 5 5 10 10 10 10 Variables (n) 30 to 100 20 to 125 40 to 100 100 200 500 50 100 250 500 Density (d %) 6.5 to 50 100 10 to 80 6.5 to 50 10 to 50 10 to 100 10 10 10 10 Linear (c ) 100 0 50 50 50 75 100 100 100 100 Coef. (c+ ) 100 63 50 50 50 75 100 100 100 100 Quadr. (q ) 100 100 100 75 100 50 100 100 100 100 Coef. (q + ) 100 0 100 75 100 50 100 100 100 100
Beasley
The density d of a quadratic pseudoBoolean function represented by the polynomial expression (2) is dened as the number of nonzero coecients cij (1 i < j n) divided by n . 2
1
Page 22
RRR 102006
6.2
Maximum cliques of graphs
Let G = (V, E) be an undirected graph, with vertex set V and edge set E. The complement G of G is the graph G = V, V \E , i.e. the graph having the same vertex set V and having 2 as edge set the complement of the edge set E of G. A clique of the graph G = (V, E) is a set of pairwise adjacent vertices. Naturally, a maximum clique (or MAXClique in short) is a clique of maximum cardinality. The size of a maximum clique is commonly called the clique number of G. We shall denote it as (G). A subset S V is called independent (or stable), if no edge of G has both endpoints in S. A maximum independent set is a largest cardinality independent set; the cardinality of a maximum independent set will be denoted by (G) and called the stability number of G. An independent set in a graph G is a clique in the complement graph G. Thus, (G) = G . Moreover, since the complement G of G has a polynomial size representation of the input size of G, and since it can be obtained in a polynomial time function of the size of G, then there is a strong equivalence between the maximum clique and the maximum independent set problems. If a solution to one of the problems is available, a solution to the other problem can be obtained immediately. It is not dicult to show (see e.g. [18]) that the cardinality of a maximum independent set of G = (V, E) is given by where (i,j) ( 0) are arbitrary nonnegative reals for all (i, j) E. Furthermore, if x = 1S is a maximizing binary vector of (21), then a maximum cardinality set S ( S) of G can be obtained in O (n) time. In fact, if the coecients e , e E, are restricted to be positive, then S = S is a maximum cardinality independent set of G. We note here that in our experiments it was assumed that (i,j) = 0 for all edges (i, j). In order to facilitate comparisons among dierent exact and heuristic methods related to clique problems, a set of benchmark graphs has been constructed in conjunction with the 1993 DIMACS Challenge on maximum cliques, coloring and satisabilty ([47]). This study only reports preprocessing results for two families of this dataset, since for the other graphs PrePro did not nd any persistencies. Namely, we consider the benchmarks containing cfat graphs or Hamming graphs: A major step in the algorithm of the fault diagnosis problem proposed by Berman and Pelc [12] is to nd the maximum clique of a special class of graphs, called cfat rings. In order to dene a cfat graph G = (V, E), let us consider an arbitrary nite set of |V | vertices V . Let c be a real parameter, k = c log|V | , and let us consider a partition W0 , , Wk1 of V , such that c log |V | |Wi | c log |V | + 1 for all i = 0, , k 1. The edge set E is dened as the set of those edges (u, v) which link distinct pairs of (G) = max
xBV
xi
1+
(i,j)
iV
(i,j)E
xi xj ,
(21)
RRR 102006
Page 23
vertices u Wi and v Wj , such that |i j| {0, 1, k 1}. The DIMACS cfat rings were created by Panos Pardalos using the cfat rings generator of Hasselberg, Pardalos and Vairaktarakis [43]. The Hamming graphs arise from coding theory problems ([81]). The Hamming distance between the binary vectors u = (u1 , , un ) and v = (v1 , , vn ) is the number of indices i = 1, , n where ui = vi . The Hamming graph H (n, d) of size n and distance d is the graph whose vertex set is the set of all binary nvectors, and whose edges link any two nvectors at distance d or larger. Clearly, the graph H (n, d) has 2n vertices, 2n1 n n edges, and the degree of each vertex is n n . A binary code i=d i i=d i consisting of a set of binary vectors, any two of which have Hamming distance greater or equal to d, can correct d1 errors. Thus, a coding theorist (see [56]) would like 2 to nd the maximum number of binary vectors of size n with Hamming distance d, i.e. the maximum clique of H (n, d). The DIMACS Hamming graphs were created by Panos Pardalos (for details see [43]).
6.3
Minimum vertex cover problems of planar graphs
The complement V \S of an independent set S of G is called a vertex cover of the graph. Let us denote the size of the smallest vertex cover of G as (G) = |V | (G). Then using (21) it can be shown that (G) = min
xBV
xi +
1+
(i,j)
iV
(i,j)E
xi xj .
(22)
The knowledge of any one of the three numbers G , (G) or (G) implies that the other two values can be immediately determined. In general nding any of these numbers is a NPhard optimization problem ([31]). It is important to note that even for planar graphs it is known that solving the minimum vertex cover problem is NPhard ([32]). Motivated by a recent work by Alber, Dorn and Niedermeier [1] we have analyzed the performance of PrePro in this class of graphs. The only dierence between the two approaches is that the method of Alber, Dorn and Niedermeier considers the inuence of a clever, VERTEX COVERspecic data reduction (see [1], page 220), whereas the results obtained by PrePro are entirely due to its ways of simplifying QUBOs, since we have not introduced any specic adaptation of this method for the case of vertex cover problems. In their experiments, [1] used a set of planar graphs randomly generated by the LEDA software package ([58]). Using the same LEDA generator we tried to replicate the experiment reported in [1], although it should be noted that not having access to the seeds used in [1], the graphs generated by us are not exactly identical to the ones used by Alber, Dorn and Niedermeier [1]. In order to distinguish between the two planar vertex cover benchmarks, we shall call those of [1] ADN benchmark graphs and the new ones PVC LEDA benchmark.
Page 24
RRR 102006
The total number of planar graphs that we have generated with LEDA is 400, partioned into 4 sets of 100 graphs, each subset having a specic number of vertices: 1000, 2000, 3000 and 4000. The planar density of each graph G (V, E) was randomly determined, i.e. |E| discrete uniform (|V | 1, 3 |V | 6). Some comparative statistical numbers about these two benchmarks are displayed in Table 2.
Table 2: Comparative statistical numbers about the LEDA benchmarks.
Benchmark PVC LEDA Vertices 1000 2000 3000 4000 1000 2000 3000 4000 Number of Graphs 100 100 100 100 100 100 100 100 Average Number of Edges 2037.9 4068.6 6204.3 8207.1 1978.9 3960.8 6070.6 8264.5 Average Maximum Degree 73.0 106.7 132.4 149.2 73.3 104.9 129.6 146.6 Average Degree 4.08 4.07 4.14 4.10 3.96 3.96 4.05 4.13 Average Minimum Vertex Cover 460.6 921.1 1391.2 1848.2 453.9 917.3 1373.8 1856.8
ADN ([1])
In addition to the PVC LEDA planar graphs we have also generated a dataset containing larger graphs with up to 500 000 vertices. These graphs were generated in order to analyze the scalability of the routine PrePro. Because of size limitations associated to our trial license on LEDA, we used for this experiment Rinaldis ([76]) generator called RUDY. With the RUDY program, we generated a total of 36 graphs whose sizes are of 50 000, 100 000, 250 000 and 500 000 vertices; for each of these graph sizes, we generated nine graphs: three instances with density of 10%, three with density of 50% and three with density of 90%. This set of benchmark graphs is called PVC RUDY.
6.4
Graphs for MAXCUT
Let G = (V, E) be an undirected graph, with vertex set V and edges set E. To each edge e E of the graph, we may assign a real weight we R. A cut of the graph G = (V, E) is dened by a partition of the vertex set V into two subsets S and S, and consists of the set of edges with exactly one endpoint in S and another in S. We shall denote the cut dened by S as S, S . The maximum cut (or MAXCUT in short) problem is to nd a cut S, S with the largest cardinality. If x = 1S is the characteristic vector representing S, then it can be shown that max S, S
SV
The weighted MAXCUT problem in a graph G = (V, E) with weights wE is to nd a cut S, S for which the sum of the weights of the corresponding edges is maximum. If the
= max
xBV
(i,j)E
(xi xj + xi xj ) .
RRR 102006
Page 25
total weight of a cut is denoted by W S, S , then the weighted maximum cut can be found by solving the problem maxW S, S = max
SV xBV
(i,j)E
wij (xi xj + xi xj ) .
(23)
Let us remark that any optimal solution x = (x1 , , xn ) of the previous problem has a complementary solution (x1 , , xn ). Thus, before calling PrePro for a MAXCUT problem, we select a variable and assign to it a 01 value. Let us note that these problems are maximization problems. Therefore, we transformed the MAXCUT problems (23) into quadratic posiform minimization problems as explained previously for the benchmarks in the xed degree dataset. We tested the PrePro algorithm on the following graphs: Torus graphs: These graphs are 3Dtoroidal graphs, originated from the Ising model of spin glasses in physics. They were taken from the DIMACS library of mixed semidenite-quadratic-linear programs2 . Via graphs ([46]): These ten graphs are derived from layer assignment problems in the design process for VLSI chips. The characteristics of these graphs, and the largest cut information can be seen in Table A.5. Sparse random graphs ([46]): These eight graphs constitute the family R of Homer and Peinado [46], each having an edge probability of 10/n, where the number of vertices n ranges from 1000 to 8000. Gn.p graphs ([51]): These graphs are randomly generated with n vertices, where an edge is placed between any two vertices with probability d, independent of all other edges. The probability d is chosen so that the expected vertex degree, i.e. d (n 1), is p. Un.p graphs ([51]): These graphs are random geometric graphs on n vertices that lie in the unit square, and whose coordinates are chosen uniformly from the unit interval. There is an edge between any two vertices if their Euclidean distance is t or less, where p = nt2 is the expected vertex degree. The general characteristics of the previous graphs, and the largest cut information are shown in Tables A.4 to A.6 of the Appendix.
The DIMACS library of mixed semidenite-quadratic-linear programs: http://dimacs.rutgers.edu/Challenges/Seventh/Instances/.
2
Page 26
RRR 102006
6.5
MAX2SAT formulas
The maximum satisability problem (or MAXSAT in short) is a popular subject in applied mathematics and computer science, and it has a natural pseudoBoolean formulation. The input of a MAXSAT instance (usually called a formula) consists of a family C of subsets C L of literals, called clauses. A binary assignment x BV satises a clause C, if at least one literal in C takes value 1 (true) for this assignment. The maximum satisability problem consists in nding a binary assignment satisfying the maximum number of clauses in C. It is easy to see that a clause C is satised by x BV if and only if uC u = 0. Thus, MAXSAT is equivalent to the problem max
CC
xBV
aC
1
uC
u,
where aC = 1 for all C C. In the weighted maximum satisability problem there is also a nonnegative weight aC R+ associated with each clause C C, and the objective is to maximize the total weight of the satised clauses. If the clauses C have at most k literals, then the (weighted) MAXSAT problem is called the MAXkSAT (weighted) problem. In particular, the weighted MAX2SAT problem can be formulated as the optimization of a special quadratic negaform 3 : () = max = max
xB V xB V
a{u} (1 u) +
{u}C
{u,v}C
a{u,v} (1 u v) .
Let us denote the sum of the nontrivial coecients of any given posiform as A () = T = aT . If the previous negaform is denoted by , then we consider the minimum of the quadratic posiform = A () , i.e. () = min = min
xB V xB V
def
a{u} u +
{u}C {u,v}C
a{u,v} u v.
(24)
Since for any assignment of the formula C, () is the maximum number of true clauses and () is the minimum number of false clauses, then it is simple to notice that () + () = A () = A (). In this work, we tested algorithm PrePro in a set of random weighted and nonweighted MAX2SAT formulas proposed by Borchers and Furman [15]. These instances are publicly available on the Internet4 .
In a similar way to a posiform, a negaform of a pseudoBoolean function f (x1 , , xn ) is dened as a polynomial g (x1 , x1 , , xn , xn ), taking the same values as f in every binary nvector, and having the property that all its coecients (with the possible exception of the free term) are nonpositive 4 MAXSAT, A Davis-Putnam like code for MAX-SAT Problems (12/18/04): http://www.nmt.edu/borchers/maxsat.html
3
RRR 102006
Page 27
The list of problems contains 17 standard formulas and 17 formulas with weights associated to the clauses, ranging from one to ten. The parameters of the non-weighted formulas can be seen in Table A.7, and those of the weighted formulas can be seen in Table A.8. We solve the (weighted) MAX2SAT problem by associating to it a quadratic posiform (see (24)), for which the minimum value () is the minimum weighted set of unsatised clauses.
7
7.1
Computational Experiments
Test environment
The algorithm PrePro was implemented in C++, compiled using the Microsoft Windows 32-bit C/C++ Optimizing Compiler (version 12) for 80x86, and linked with the Microsoft Incremental Linker (version 6). Except for the experiments on large scale planar graphs (previously dened as the RUDY benchmark), all the remaining computations were carried out on a computer with a 2.80 GHz Pentium 4, 512 MB of RAM, and hyper-threading technology. The experiments on the large planar graphs were carried out on a computer system with 3.5 GB of RAM, and a 3.06 GHz Xeon. Both computers have installed the Windows XP Professional (version 2002) operating system.
7.2
Results
Given a posiform , its roofdual bound, a set of strong/weak persistencies and a decomposition of type (16) of it, can easily be derived from the residual network resulting after aplying a maximum ow algorithm to the implication network associated to . We consider this preprocessing step (entirely based on the roofduality theory) as a standard preprocessing tool in all experiments that we have carried out. We tested four possible preprocessing strategies: A Only the standard tool is considered; B Standard tool and coordination are considered; C Standard tool and probing are considered; and D All preprocessing tools are considered. Since strategy D usually provides the best data reduction of the problems, we have included in Appendix B several statistical results about the preprocessing performance of this strategy in all benchmarks. At rst glance, we have tried to understand how any of the preprocessing techniques would impact in the reduction of the problems size. Table 3 provides preprocessing results for the test beds, whose values are averages for groups of problems belonging to the same
Page 28
RRR 102006
family. Strategy D provides 100% data reduction for the following benchmarks: MAXClique problems in all HAM2 graphs and all cFAT graphs; minimum vertex cover problems in all PVC LEDA planar graphs; MAXCUT problems in all VIA graphs. These results clearly indicate that one can expect getting an oustanding data reduction level in these special well structured problems. It should be remarked that the border separating successful from unsuccessful preprocessing cases is very thin. For instance, all the Hamming graphs in HAM2 were optimally solved with the standard preprocessing tool. However, in the closely related family HAM4 there was no reduction found for any of the graphs, even when all the preprocessing tools were considered. We also remark the fact that strategy C provided optimallity for all MAXClique problems and all MAXCUT problem in the VIA graphs, using a substantial smaller computing time than the one corresponding to strategy D. Strategy C with an average value of 99.9%, provided also a very good data reduction on the minimum vertex cover problems in the PVC LEDA planar graphs. Table 4 suggests the particular preprocessing techniques which can be recommended for each of the problem families considered, in order to achieve (on the average) as high a data reduction as possible within a limited amount of time. In view of the relatively uniform behavior of problems within the same family, the recommendation of a common strategy for problems within a family seems reasonable. Here are some remarks and some recommendations for the examined groups: Problems with prescribed density Coordination does not have a practical inuence in the preprocessing results. Probing should be used in the cases where density is low. In general, the probing tool should be used in this class, if the condition nd 20 is satised. Let us also note that family B consists of very dense submodular maximization problems, and for which the preprocessing outcome changed considerably, in comparison with the other problems. Six of the 10 problems in the B group were solved optimally, and for the unsolved cases, a large number of quadratic persistencies was found. Minimum vertex cover of planar graphs Probing when used with the coordination method provides slighty better preprocessing data reduction, without degrading the computing times returned by probing only. MAXCUT in torus graphs The standard preprocessing tool should be used for the graphs with 1 weights in the edges (see also Table B.5). Probing should be used in the other weighted graphs. MAXCUT in VIA graphs All the problems in the VIA.CY family are solved optimally by the basic preprocessing tool (see also Table B.5). Every instance in this group of problems is solved in less than 0.2 seconds. The VIA.CN problems are all optimally solved with probing, taking an average computing time of 13.2 seconds. In two of the VIA.CN problems, the analysis of the starting implication network found 2 components which were preprocessed separately under the result of Proposition 16.
RRR 102006
Table 3: Average QUBO simplications and decomposition after preprocessing.
Type of Problem Family Name A B C D E F1 ORL-50 ORL-100 ORL-250 ORL-500 C-FAT-200 C-FAT-500 HAM-2 HAM-4 LEDA-1000 LEDA-2000 LEDA-3000 LEDA-4000 Torus R VIA.CN VIA.CY G500 G1000 U500 U1000 BF-50 BF-100 BF-150 BFW-50 BFW-100 BFW-150 Roof-Duality Total Variab. Time Reduc. 0.0 52.3% 0.0 0.9% 0.0 18.1% 0.0 0.6% 0.0 0.0% 1.1 0.0% 0.0 94.2% 0.0 4.9% 0.0 0.0% 0.1 0.0% 0.0 0.0% 0.1 0.0% 0.0 100.0% 0.1 0.0% 0.0 75.3% 0.1 74.5% 0.2 75.0% 0.4 75.1% 0.1 0.2% 0.4 0.0% 0.1 4.0% 0.1 100.0% 0.0 4.2% 0.1 2.7% 0.0 1.0% 0.1 0.6% 0.0 9.3% 0.0 12.6% 0.0 17.1% 0.0 6.2% 0.0 14.2% 0.0 17.1% Preprocessing Tools Used: Roof-Duality and Coordination Roof-Duality and Probing Total Relat. Quad. Variab. Total Relat. Quad. Variab. Time Gap Rel. Reduc. Time Gap Rel. Reduc. 0.0 8.2% 1 52.9% 0.0 5.1% 1 64.0% 0.1 85.3% 1762 0.9% 1.2 34.9% 785 72.2% 0.0 22.4% 0 18.5% 0.1 19.8% 1 30.2% 0.0 56.8% 0 0.6% 1.2 55.4% 3 1.6% 0.1 57.9% 0 0.0% 5.2 57.2% 6 0.0% 1.6 78.7% 0 0.0% 149.8 78.6% 0 0.0% 0.0 0.2% 2 94.4% 0.0 0.0% 0 100.0% 0.0 13.3% 0 5.2% 0.3 8.8% 26 36.3% 0.1 44.1% 0 0.0% 3.3 43.5% 1 0.0% 0.4 60.6% 0 0.0% 28.4 60.3% 0 0.0% 1.3 68.7% 29 0.0% 4.8 0.0% 0 100.0% 21.9 77.0% 56 0.0% 80.8 0.0% 0 100.0% 0.0 0.0% 0 100.0% 0.0 0.0% 0 100.0% 12.8 89.6% 0 0.0% 79.0 88.3% 0 0.0% 0.1 0.0% 0 99.9% 0.1 0.0% 0 99.9% 0.2 0.0% 0 99.9% 0.2 0.0% 0 99.8% 0.3 0.0% 0 99.9% 0.3 0.0% 0 99.8% 0.5 0.0% 0 99.8% 0.6 0.0% 0 99.9% 6.4 39.2% 0 0.2% 220.0 38.3% 2 3.1% 64.7 28.9% 0 0.0% 1818.1 28.7% 1 0.2% 1.1 4.6% 3 4.0% 13.2 0.0% 0 100.0% 0.1 0.0% 0 100.0% 0.1 0.0% 0 100.0% 0.5 23.1% 3 4.2% 6.2 22.3% 5 19.0% 2.2 23.5% 5 2.7% 36.1 23.0% 7 18.8% 2.6 35.2% 31 2.6% 9.1 36.0% 0 5.7% 14.7 35.2% 57 2.5% 39.0 35.9% 0 7.6% 0.0 240.0% 1 11.3% 0.1 121.2% 8 25.8% 0.0 691.1% 1 13.8% 0.1 271.9% 24 27.6% 0.0 908.1% 0 17.1% 0.3 258.3% 55 42.2% 0.0 391.4% 1 6.2% 0.1 133.2% 8 24.2% 0.0 1731.4% 5 14.2% 0.2 335.5% 39 26.2% 0.0 3214.3% 4 17.1% 0.7 263.1% 68 45.6% Total Time 0.0 2.9 0.1 1.2 5.3 150.1 0.0 0.4 3.4 28.5 14.1 327.4 0.0 91.5 0.1 0.2 0.3 0.5 853.3 3187.8 100.5 0.1 11.8 100.8 14.3 76.2 0.1 0.2 0.5 0.1 0.3 0.9 ALL Tools Relat. Quad. Gap Rel. 5.1% 1 41.9% 801 19.8% 1 55.4% 3 57.2% 6 78.6% 0 0.0% 0 8.8% 26 43.5% 1 60.3% 0 0.0% 0 0.0% 0 0.0% 0 88.3% 0 0.0% 0 0.0% 0 0.0% 0 0.0% 0 38.3% 2 28.7% 1 0.0% 0 0.0% 0 22.3% 5 23.0% 7 33.8% 32 34.2% 83 121.2% 8 271.9% 25 258.3% 55 133.2% 8 335.5% 39 265.2% 67 Variab. Reduc. 64.0% 68.2% 30.2% 1.6% 0.0% 0.0% 100.0% 36.3% 0.0% 0.0% 100.0% 100.0% 100.0% 0.0% 100.0% 100.0% 100.0% 100.0% 3.1% 0.2% 100.0% 100.0% 19.0% 18.8% 11.5% 9.9% 25.8% 27.6% 42.2% 24.2% 26.2% 45.8%
Fixed Degree
MAX Clique
MIN Vert. Cov.
MAX Cut
MAX 2Sat
Page 29
Page 30
RRR 102006
Table 4: Preprocessing strategies recommended for the benchmarks.
Type of Problem Family Name A B C D E F1 ORL-50 ORL-100 ORL-250 ORL-500 C-FAT-200 C-FAT-500 HAM-2 HAM-4 LEDA-1000 LEDA-2000 LEDA-3000 LEDA-4000 Torus R VIA.CN VIA.CY G500 G1000 U500 U1000 BF-50 BF-100 BF-150 BFW-50 BFW-100 BFW-150 Number of Instances 8 10 7 10 5 5 10 10 10 10 3 4 3 3 100 100 100 100 4 8 5 5 4 4 4 4 9 5 3 9 5 3 Density 18.78% 98.83% 36.01% 54.31% 29.62% 51.56% 9.75% 9.75% 9.94% 9.90% 77.82% 83.28% 4.55% 39.42% 0.41% 0.20% 0.14% 0.10% 0.68% 0.34% 0.32% 0.37% 1.87% 0.95% 3.40% 1.72% 19.98% 7.53% 3.91% 21.18% 7.77% 3.93% Best Strategy C C C A A A C C A A C C A A D D D D C C C A C C D D C C C C C C
Fixed Degree
MAX Clique
MIN Vertex Cover
MAX Cut
MAX 2SAT
RRR 102006
Page 31
MAXCUT in Gn.p graphs The preprocessing eciency decreases with density for the graphs with the same number of vertices (see also Table B.5). Probing helped increasing data reduction for graphs with densities up to 5%, and attained as expected better performance for graphs with 500 vertices, than for those with 1000 vertices. MAXCUT in Un.p graphs The preprocessing eciency decreases with density for the graphs with the same number of vertices (see also Table B.5). In this category, the standard preprocessing tool found some non trivial decomposition, and both coordination and probing helped improving the average data reduction rates from 23% to about 11%. MAX2SAT The preprocessing eciency decreases with the number of clauses when the number of variables is xed (see also Table B.6). Both in the nonweighted and weighted formulas, the probing technique provided better reduction indicators. In conclusion it can be seen that the choice best strategy is highly problem family dependent. It should also be remarked that only three out of the four examined strategies turn out to provide the best performance for some of the considered group of problems; strategy B (consisting of the application of the standard tool and coordination, but not of probing) did not give best results in any of the examined cases. Table 4 indicates the best recommended strategies for each of the examined families of problems.
8
Application: Optimal vertex covers of planar graphs
In view of the outstanding results obtained by applying PrePro to the minimum vertex cover problem in random planar graphs, we have tried to rene this method to the point where it would not only preprocess the problem but actually produce an optimal solution of it. As it will be seen in this section, the resulting method allowed the ecient detection of minimum vertex covers in planar graphs of impressive dimensions, including some having 500 000 vertices. Although the vertex cover problem is known to be NPhard in the class of planar graphs ([32]), our computational experiments with a large collection of benchmark planar graphs indicate that, in all likelihood, nding vertex covers in planar graphs may be frequently tractable. This conclusion provides the motivation for the work reported in this section. Before presenting the results of our computational experiments we would like to emphasize that PrePro is not capable of solving the QUBO problems associated to every planar graph, and that it may encounter problems even in the case of very small graphs. For example there are no persistencies in the QUBO associated to the toy box graph of Figure 6.
Page 32
RRR 102006
(a) 3dimensional drawing.
(b) Planar drawing.
Figure 6: Planar graph for which PrePro does not nd any persistent result.
8.1
Deriving minimum vertex cover from QUBOs optimum
We have seen in Section 6.3 that nding a minimum vertex cover of a graph G = (V, E) is a QUBO, and we have also noticed in Section 7.2 that out of the 400 QUBOs coming from vertex cover problems in randomly generated planar graphs, every single QUBO was solved to optimality by PrePro. The only matter which needs special attention is that due to the fact that in (22) we have xed to zero the values of (i,j) for every edge (i, j) it may happen that the optimal 01 solution of a QUBO denes a vertex set which does not cover every edge. Let us prove now that there is a simple polynomial time transformation, which associates to the optimal solution of QUBO an optimal vertex cover of G. Proposition 10. Let G = (V, E) be a graph, and let us associate to it the quadratic pseudo Boolean function f (x1 , , xn ) = iV xi + (i,j)E xi xj and the QUBO (22). Let further f (x1 , , x ) be a minimum of f , and let S be the set of vertices j V for which x = 1. j n Then, the size of a minimum vertex cover is f (x , , x ), and the set S can be enlarged 1 n to a minimum vertex cover S S in O (|E|) time. Proof. If S is a vertex cover, it is clearly minimal, and since x x = 0 for every (i, j) E, ij the problem is solved. Let us assume that S is not a vertex cover, and let h and k be two adjacent vertices with x = x = 0. Let us consider now the vector (x , , x ) dened by n 1 k h x = l
xl if l = k . 1 if l = k
If Nk is the neighborhood of k in the set V \ S , then clearly
f (x , , x ) = f (x , , x ) + 1 1 n 1 n
tNk
x , t
and since h Nk , the set Nk is not empty, and therefore
f (x , , x ) 1 n
f (x , , x ) . 1 n
RRR 102006 On the other hand, from the minimality of f in (x , , x ) it follows that 1 n f (x , , x ) = f (x , , x ) . 1 n 1 n
Page 33
By repeating this transformation several times (at most |V \ S | times), the set S will be eventually transformed to a vertex cover S of the graph G, which has to have a minimum size. While the above proposition holds in any graph, it is particularly useful in classes of graphs for which the associated QUBO (22) can be solved to optimality. As a consequence of the discussion above, we have supplemented PrePro with the simple procedure outlined in the proof of Proposition 10 to derive an optimal vertex cover from the optimal solution of the corresponding QUBO problem (22). This amended version of the proposed algorithm will be called PrePro+ . In this section the preprocessing strategies A, B and C of the previous section will not be considered, i.e. all the experiments below were carried out by using strategy D, which involves all the preprocessing tools of Section 7.2. The results obtained by PrePro+ for moderately sized graphs (i.e. having up to a few thousand vertices), have been compared with those of the recent paper of Alber, Dorn and Niedermeier (or ADN in short) reported in [1], which is essentially based on the data reduction results of Nemhauser and Trotter [61].
8.2
Preprocessing
Table 5 provides averages of results obtained by preprocessing minimum vertex cover problems on 400 random planar graphs generated by the LEDA package (see Section 6.3). Four groups of 100 graphs each have been considered, each graph in these sets containing respectively 1000, 2000, 3000 and 4000 vertices.
Table 5: Comparative preprocessing results for minimum vertex cover problems in planar graphs. Number of Time Variables Fixed Size of Residual Graphs Vertices (sec) (%) Problem per Family per Graph ADN ([1]) PrePro ADN ([1]) PrePro ADN ([1]) PrePro 100 1 000 4.06 0.05 68.4 100 315.8 0 100 2 000 12.24 0.16 67.4 100 652.9 0 100 3 000 30.90 0.27 65.5 100 1036.1 0 100 4 000 60.45 0.53 62.7 100 1492.9 0
Remarkably, PrePro achieved 100% data reduction in all PVC LEDA graphs, whereas the ADN method obtained between 63% and 68% average data reduction on their LEDA benchmarks, which have similar characteristics to the PVC LEDA graphs (see Table 2). It can also be seen that the best performance of the ADN method (68.4% reduction of vertex set) occurs on the group of relatively smaller graphs, while the performance of PrePro
Page 34
RRR 102006
(100% reduction of vertex set) does not seem to be inuenced by the size of the graphs considered.
8.3
Optimization
While the results reported in Table 5 refer to the preprocessing by PrePro of the minimum vertex cover problem, we shall discuss below the results of applying PrePro+ for actually nding optimal solutions to this problem. It is important to remark that PrePro+ assumes the knowledge of the exact optimum of the associated QUBO. If this optimum is not available PrePro+ is further enhanced to an algorithm PrePro , by adding a branchandbound component to it, in order to handle minimum vertex cover problems even in this case. However, the use of PrePro turned out not to be necessary in any of the 400 test problems randomly generated with the LEDA software package, which were all solved to optimality without the branchandbound component having been called.
Table 6: Average computing times of optimal vertex covers for graphs belonging to the LEDA benchmarks. Algorithm ADN ([1]) PrePro in the PVC LEDA Benchmark 750 MHz 500 MHz Pentium III 2.8 GHz Pentium 4 Computer Linux Windows 98 Windows XP System 720 MB RAM 96 MB RAM 512 MB RAM
(speed) (slower) (faster)
1 000 vertices 2 000 vertices 3 000 vertices 4 000 vertices Average Speedup
5.75 19.93 51.54 109.84
sec sec sec sec
0.24 0.64 1.07 1.71
sec sec sec sec
0.06 0.18 0.31 0.56
sec sec sec sec
51 times
169 times
Table 6 provides comparative results for nding optimal vertex covers for graphs belonging to the LEDA benchmarks. It includes computing times for the exact algorithm of Alber, Dorn and Niedermeier [1] and solution times for PrePro (which coincide with PrePro+ for all the 400 cases). We would like to recall the fact that not having had access to the test problems of [1] we have randomly generated our problems, but made sure (as explained in Section 6.3) that the parameters used in the random graph generation process were chosen so as to match exactly does of [1]. In order to be able to dierentiate between the acceleration due to computer systems and those due to algorithms, all the experiments reported in Table 6 have been carried out twice, rst on a somewhat slower computer system (500 MHz Pentium III, 98 MB RAM, Windows 98) than the one used by [1] (715 MHz, 720 MB RAM, Linux), and second on a faster system (2.8 GHz Pentium 4, 512 MB RAM, Windows XP).
RRR 102006
Page 35
The basic conclusion of this set of experiments is that using the slower computer system, PrePro is about 50 times faster than that of [1], on average.
8.4
Minimum vertex covers of very large planar graphs
Based on the high eciency of PrePro when applied to the optimization of vertex cover problems in planar graphs, we have investigated the possibility of using it on substantially larger planar graphs. The relevant experiments were carried out on the set of 36 benchmark problems contained in the RUDY list (described in Section 6.3), which contains graphs whose vertex sets contain 50 000, 100 000, 250 000 and 500 000 vertices, and have planar densities of 10%, 50% and 90%. For each particular number of vertices and each density the list contains three graphs.
Table 7: Average computing times over 3 experiments of optimal vertex covers for graphs belonging to the PVC RUDY benchmark. Planar Density Vertices 10% 50% 90% 50 000 1.2 min 3.7 min 1.8 min 100 000 4.8 min 16.2 min 7.4 min 250 000 30.4 min 107.7 min 48.2 min 500 000 124.7 min 422.4 min 195.3 min
Table 7 presents the average computing times needed by PrePro for nding minimum vertex cover sets in all the graphs contained in the RUDY list. Each of the computing times reported in the table represents the average needed for solving the three problems with a xed number of vertices and a xed planar density contained in the RUDY list. The average computing times range from 2.2 minutes for the graphs with 50 000 vertices up to 4.1 hours for the graphs with 500 000 vertices. Clearly, the computing times vary with the size of the vertex set. A similarly foreseable phenomenon happens with the dependency of computing times and densities. Indeed, the average computing time for the low density graphs is 40 minutes, for medium density graphs this increases to 2.3 hours, and for high density graphs it drops to 1 hour. More detailed information about the performance of PrePro can be read from the statistics presented in Table B.4 in the Appendix, where specic data are given for each of the 36 problems of the RUDY list. First, it can be seen that almost all of the computing time (78.7%) is spent on calculating the roof duals; moreover, most of this time (99.9%) is spent on calculating the very rst roof dual. The large investment of computing time in the calculation of roof duals brings returns in the form of graph size reductions (which are due to strong and weak persistency) and in the form of decompositions. The detailed analysis of the problem size reductions occurring in PrePro shows that
Page 36
RRR 102006
roofduality accounts for 99.8% of these reductions for planar graphs of density 10%, 93.2% for the 50% dense graphs, and 51.8% for the 90% dense graphs. It is interesting to note the extremely small size of the average components of the graphs left after applying decomposition and strong and weak persistency. Indeed, the average size of these components for graphs of 10%, 50% and 90% density is of 3.1, 4.4 and 14.4 vertices, respectively. Beside roofduality, important simplications of the remaining QUBOs were obtained by the coordination method and by probing. It can be seen in column (ne ) of Table B.4 of the Appendix that the number of equality relations between pairs of variables or their complements, discovered by the coordination method is an increasing monotone function of planar density. Also, column (nf ) of Table B.4 shows that the number of variables whose values are xed by probing reaches maximum values for the medium density graphs. In conclusion it can be seen that there is substantial complementarity in the eect of applying the basic preprocessing techniques considered in this paper. Indeed, 10% dense planar graphs derive almost the entire solution from the application of roofduality; 50% dense planar graphs derive a considerable reduction through probing; and 90% dense planar graphs derive a considerable reduction through the coordination method. However, the most important conclusion is that PrePro+ found optimal vertex covers for all the 36 benchmarks in the RUDY list.
9
Conclusions
This study is devoted to the systematic simplication of QUBOs. The proposed method uses enhanced versions of several basic techniques (e.g., extraction of conclusions from the analysis of rst and second order derivatives [39], and from roofduality [40]) and several integrative techniques (e.g., probing, consensus) for combining the conclusions provided by the basic techniques. The application of these techniques is implemented using the network ow model of [18, 20]. The use of the proposed preprocessing techniques provides: (i) A lower bound of the minimum of the objective function; (ii) The values of some of the variables in some or every optimum; (iii) Binary relations (equations, inequalities, or non-equalities) between the values of certain pairs of variables in some or every optimum;
RRR 102006
Page 37
(iv) The decomposition (if possible) of the original problem into several smaller pairwise independent minimization problems. The eciency of the proposed techniques is demonstrated through numerous computational experiments carried both on benchmark problems and on randomly generated ones. The simplications obtained with the proposed methods exceed substantially those reported in the literature. An interesting example is the minimum vertex cover problem for which [1] reports a preprocessing stage reduction of dimensionality by 62.7%68.4%, while the method proposed here achieves 100% reduction (i.e. exact solution) in each of the test problems. Moreover, while the computing times reported in [1] for nding optimal vertex covers for graphs from 1 000 to 4 000 vertices range from 5.75 to 109.84 seconds, those required by the proposed method range from 0.24 to 1.71 seconds using a somewhat slower computer, or from 0.06 to 0.56 seconds using a somewhat faster one. The experiments show that the methods can be applied successfully to problems of unusually large size, for instance: MAXClique on graphs derived from fault diagnosis having up to 500 vertices; MAXCUT problems on graphs derived from VLSI design having thousands of vertices; Minimum vertex cover problems on randomly generated planar graphs (an NP-hard problem [32]) having up to 500 000 vertices. It should be added that the proposed preprocessing technique have not only simplied the above problems but have in fact produced their exact optimum solutions. As far as we know there are no reports in the literature about methods capable of providing optimal solutions to vertex cover problems in planar graphs with the investigated sizes.
References
[1] Alber, J., F. Dorn and R. Niedermeier. Experimental evaluation of a tree decomposition based algorithm for vertex cover on planar graphs. Discrete Applied Mathematics 145, (2005), pp. 219231. [2] Alidaee, B., G.A. Kochenberger and A. Ahmadian. 0-1 quadratic programming approach for the optimal solution of two scheduling problems. International Journal of Systems Science 25, (1994), pp. 401408. [3] Amini, M.M., B. Alidaee and G.A. Kochemberger. A scatter search approach to unconstrained quadratic binary programs. In New ideas in optimisation (D. Corne, M. Dorigo and F. Glover, eds.), pp. 317329 (McGraw-Hill, London, 1999). [4] Aspvall, B., M.F. Plass and R.E. Tarjan. A linear time algorithm for testing the truth of certain quantied boolean formulas. Information Processing Letters 8, (1979), pp. 121123.
Page 38
RRR 102006
[5] Atamt rk, A., G.L. Nemhauser and M.W.P. Savelsbergh. Conict graphs in solving u integer programming problems. European Journal of Operational Research 121, (2000), pp. 4055. [6] Axehill, D. and A. Hansson. A preprocessing algorithm for MIQP solvers with applications to MPC. Tech. rep., Department of Electrical Engineering, Linkping University o (2004). [7] Badics, T. Approximation of some nonlinear binary optimization problems. Ph.D. thesis, RUTCOR, Rutgers University (1996). [8] Badics, T. and E. Boros. Minimization of half-products. Mathematics of Operations Research 23, (1998), pp. 649660. [9] Barahona, F., M. Grtschel, M. J nger and G. Reinelt. An application of combinatorial o u optimization to statistical physics and circuit layout design. Operations Research 36, (1988), pp. 493513. [10] Barahona, F., M. J nger and G. Reinelt. Experiments in quadratic 0-1 programming. u Mathematical Programming 44, (1989), pp. 127137. [11] Beasley, J.E. Heuristic algorithms for the unconstrained binary quadratic programming problem. Tech. rep., Management School, Imperial College, London, UK (1998). [12] Berman, P. and A. Pelc. Distributed fault diagnosis for multiprocessor systems. In Proceedings of the 20th annual international symposium on fault-tolerant computing (Newcastle, UK, 1990), pp. 340346. [13] Billionnet, A. and B. Jaumard. A decomposition method for minimizing quadratic pseudo-boolean functions. Operations Research Letters 8(3), (1989), pp. 161163. [14] Billionnet, A. and A. Sutter. Minimization of a quadratic pseudo-boolean function. European Journal of Operational Research 78, (1994), pp. 106115. [15] Borchers, B. and J. Furman. A two-phase exact algorithm for max-sat and weighted max-sat problems. Journal of Combinatorial Optimization 2(4), (1998), pp. 299306. [16] Boros, E. and P.L. Hammer. A max-ow approach to improved roof-duality in quadratic 0-1 minimization. Research Report RRR 15-1989, RUTCOR, Rutgers University (1989). [17] Boros, E. and P.L. Hammer. The max-cut problem and quadratic 0-1 optimization, polyhedral aspects, relaxations and bounds. Annals of Operations Research 33, (1991), pp. 151180. [18] Boros, E. and P.L. Hammer. Pseudo-boolean optimization. Discrete Applied Mathematics 123, (2002), pp. 155225.
RRR 102006
Page 39
[19] Boros, E., P.L. Hammer, M. Minoux and D. Rader. Optimal cell ipping to minimize channel density in vlsi design and pseudo-boolean optimization. Discrete Applied Mathematics 90, (1999), pp. 6988. [20] Boros, E., P.L. Hammer, R. Sun and G. Tavares. A max-ow approach to improved lower bounds for quadratic 0-1 minimization. Research Report RRR 7-2006, RUTCOR, Rutgers University (2006). [21] Boros, E., P.L. Hammer and X. Sun. The DDT method for quadratic 0-1 optimization. Research Report RRR 39-1989, RUTCOR, Rutgers University (1989). [22] Boros, E., P.L. Hammer and G. Tavares. Probabilistic One-Pass Heuristics for Quadratic Unconstrained Binary Optimization (QUBO). Technical Report 1-2006, RUTCOR, Rutgers University (2006). [23] Burer, S., R.D.C. Monteiro and Y. Zhang. Rank-two relaxation heuristics for max-cut and other binary quadratic programs. SIAM Journal on Optimization 12(2), (2001), pp. 503521. [24] Bushnell, M.L. and I.P. Shaik. Robust delay fault built-in self-testing method and apparatus. United States Patent # 5,422,891 . [25] Carter, M.W. The indenite zero-one quadratic problem. Discrete Applied Mathematics 7, (1984), pp. 2344. [26] Crama, Y. and P.L. Hammer. Boolean functions: Theory, algorithms and applications (Cambridge University Press, 2006, forthcoming). [27] Crama, Y. and J.B. Mazzola. Valid inequalities and facets for a hypergraph model of the nonlinear knapsack and fms part-selection problems. Annals of Operations Research 58, (1995), pp. 99128. [28] Festa, P., P.M. Pardalos, M.G.C. Resende and C.C. Ribeiro. Randomized heuristics for the max-cut problem. Optimization Methods and Software 7, (2002), pp. 10331058. [29] Fraenkel, A.S. and P.L. Hammer. Pseudo-boolean functions and their graphs. Annals of Discrete Mathematics 20, (1984), pp. 137146. [30] Gallo, G., P.L. Hammer and B. Simeone. Quadratic knapsack problems. Mathematical Programming 12, (1980), pp. 132149. [31] Garey, M.R. and D.S. Johnson. Computers and intractability: An introduction to the theory of NP-completeness (W. H. Freeman, San Francisco, 1979). [32] Garey, M.R., D.S. Johnson and L. Stockmeyer. Some simplied np-complete graph problems. Theoretical Computer Science 1, (1976), pp. 237267.
Page 40
RRR 102006
[33] Glover, F., B. Alidaee, C. Rego and G. Kochenberger. One-pass heuristics for large-scale unconstrained binary quadratic problems. European Journal of Operational Research 137, (2002), pp. 272287. [34] Glover, F., G. Kochenberger and B. Alidaee. Adaptative memory tabu search for binary quadratic programs. Management Science 44(3), (1998), pp. 336345. [35] Glover, F., G.A. Kochenberger, B. Alidaee and M. Amini. Tabu search with critical event memory: An enhanced application for binary quadratic programs. In Metaheuristics - Advances and trends in local search paradigms for optimization (S. Voss, S. Martello, I. Osman and C. Roucairol, eds.), pp. 83109 (Kluwer Academic Publishers, 1998). [36] Gulati, S.K., S.K. Gupta and A.K. Mittal. Unconstrained quadratic bivalent programming problem. European Journal of Operational Research 15, (1980), pp. 121125. [37] Hammer, P.L. Plant location - a pseudo-boolean approach. Israel Journal of Technology 6, (1968), pp. 330332. [38] Hammer, P.L. Pseudo-boolean remarks on balanced graphs. International Series of Numerical Mathematics 36, (1977), pp. 6978. [39] Hammer, P.L. and P. Hansen. Logical relations in quadratic 0-1 programming. Romanian Journal of Pure and Applied Mathematics 26(3), (1981), pp. 421429. [40] Hammer, P.L., P. Hansen and B. Simeone. Roof duality, complementation and persistency in quadratic 0-1 optimization. Mathematical Programming 28, (1984), pp. 121155. [41] Hammer, P.L. and S. Rudeanu. Boolean methods in operations research and related areas (Springer-Verlag, Berlin, Heidelberg, New York, 1968). [42] Hammer, P.L. and E. Shlier. Applications of pseudo-boolean methods to economic problems. Theory and Decision 1, (1971), pp. 296308. [43] Hasselberg, J., P.M. Pardalos and G. Vairaktarakis. Test case generators and computational results for the maximum clique problem. Journal of Global Optimization 3(4), (1993), pp. 463482. [44] Helmberg, C. and F. Rendl. Solving quadratic (0,1)-problems by semidenite programs and cutting planes. Mathematical Programming 82, (1998), pp. 291315. [45] Hillier, F.S. The evaluation of risky interrelated investments (North-Holland, Amsterdam, 1969).
RRR 102006
Page 41
[46] Homer, S. and M. Peinado. Design and performance of parallel and distributed approximation algorithms for maxcut. Journal of Parallel and Distributed Computing 46, (1997), pp. 4861. [47] Johnson, D.S. and M.A. Trick. Cliques, coloring, and satisability. In 2nd DIMACS Implementation Challenge in New Brunswick, NJ, October 11-13, 1993 (D.S. in Discrete Mathematics and Theoretical Computer Science (26), eds.) (American Mathematical Society, Providence, RI, 1996). [48] J nger, M., A. Martin, G. Reinelt and R. Weismantel. Quadratic 0-1 optimization u and a decomposition approach for the placement of electronic circuits. Mathematical Programming 63, (1994), pp. 257279. [49] Kalantari, B. and A. Bagchi. An algorithm for quadratic zero-one programs. Naval Research Logistics 37, (1990), pp. 527538. [50] Katayama, K. and H. Narihisa. Performance of simulated annealing-based heuristic for the unconstrained binary quadratic programming problem. European Journal of Operational Research 134, (2001), pp. 103119. [51] Kim, S.H., Y.H. Kim and B.R. Moon. A hybrid genetic algorithm for the max cut problem. In Genetic and Evolutionary Computation Conference (2001), pp. 416423. [52] Krarup, J. and P.M. Pruzan. Computer-aided layout design. Mathematical Programming Study 9, (1978), pp. 7594. [53] Kubiak, W. New results on the completion time variance minimization. Discrete Applied Mathematics 58, (1995), pp. 157168. [54] Laughhunn, D.J. Quadratic binary programming with applications to capital budgeting problems. Operations Research 18, (1970), pp. 454461. [55] Laughhunn, D.J. and D.E. Peterson. Computational experience with capital expenditure programming models under risk. J. Business Finance 3, (1971), pp. 4348. [56] MacWilliams, F.J. and N.J.A. Sloane. The theory of error-correcting codes (NorthHolland, Amsterdam, 1979). [57] McBride, R.D. and J.S. Yormark. An implicit enumeration algorithm for quadratic integer programming. Management Science 26, (1980), pp. 282296. [58] Mehlhorn, K. and S. Nher. The LEDA platform of combinatorial and geometric coma puting (Cambridge University Press, Cambridge, England, 1999). [59] Merz, P. and B. Freisleben. Genetic algorithms for binary quadratic programming. In Proceedings of the 1999 international Genetic and Evolutionary Computation Conference (GECCO99) (1999), pp. 417424.
Page 42
RRR 102006
[60] Merz, P. and B. Freisleben. Greedy and local search heuristics for the unconstrained binary quadratic programming problem. Journal of Heuristics 8(2), (2002), pp. 197 213. [61] Nemhauser, G.L. and L.E. Trotter. Vertex packings: Structural properties and algorithms. Mathematical Programming 8, (1975), pp. 232248. [62] Palubeckis, G. A heuristic-based branch and bound algorithm for unconstrained quadratic zero-one programming. Computing 54, (1995), pp. 283301. [63] Palubeckis, G. Steepest ascent: A generalization of the greedy algorithm for the maximum clique problem. Information Technology and Control 28(3), (2003), pp. 713. [64] Palubeckis, G. Multistart tabu search strategies for the unconstrained binary quadratic optimization problem. Annals of Operations Research 131, (2004), pp. 259282. [65] Palubeckis, G. and V. Krivickiene. Application of multistart tabu search to the max-cut problem. Information Technology and Control 31(2), (2004), pp. 2935. [66] Palubeckis, G. and A. Tomkevi`ius. Grasp implementations for the uncostrained binary e quadratic optimization problem. Information Technology and Control 24(3), (2002), pp. 1420. [67] Papaioannou, S.G. Optimal test generation in combinational networks by pseudoboolean programming. IEEE Transactions on Computers 26, (1977), pp. 553560. [68] ...