A hard integer program made easy by lexicography

Noname manuscript No. (will be inserted by the editor)

A hard integer program made easy by lexicography Egon Balas · Matteo Fischetti · Arrigo Zanette

October 12, 2010

Abstract A small but notoriously hard integer program formulated by Donald Knuth fifty years ago is solved by three versions of a lexicographic algorithm using Gomory cuts. The lexicographic cutting plane algorithms are faster than CPLEX on this problem by a factor of at least 20. Keywords Cutting Plane Methods · Gomory Cuts · Degeneracy in Linear Programming · Lexicographic Dual Simplex · Computational Analysis

1 Introduction During a recent visit to Carnegie Mellon University, Donald Knuth mentioned in a conversation with one of the authors of this note that early in his career he formulated an integer program that could not be solved with the technology of the time and remained unsolved for another 35 years. Upon being asked for the details of the case, he then kindly provided the information that follows. As an undergraduate student at Case Institute of Technology, in 1960 Donald Knuth wrote a paper [2] that gave an integer programming formulation to the problem of minimizing drum latency time, i.e. optimizing the allocation of memory positions for the instructions and data of a program stored on a rotating magnetic drum. The formulation used 51 integer variables and 43 inequalities, besides the nonnegativity constraints on each variable. The young Knuth also implemented Gomory’s algorithm, or, to put it in his own words, “I used Gomory’s original all-integer algorithm, based on his preprint dated January 29, 1960; he told me my implementation was probably the E. Balas Carnegie Mellon University, Pittsburgh, PA E-mail: [email protected] M. Fischetti DEI, University of Padova E-mail: [email protected] A. Zanette DEI, University of Padova E-mail: [email protected]

2

first outside of Yorktown Heights” [3]. He then tried to solve the instance on an IBM 650 and of course failed. In spite of its relatively modest size and the obvious historical interest that it presented, the problem remained unsolved for the next 35 years. Finally, in 1995 Dimitris Alevras at ZIB (Berlin) solved it by running it on CPLEX, on a SPARCstation 5. A solution of value of 22,996 was found after generating 8,880 nodes of a branch-and-bound tree, but it took 732,200 nodes of the search tree to prove that solution optimal. Alevras reported that the problem was “amazingly sensitive from a numerical point of view” [3]. Another run, with different CPLEX parameter settings for branching, and using SOS, clique and cover inequalities, needed more than 22,000 nodes just to find a first integer solution.

2 The experiment Having recently discovered [5, 1] that a lexicographic version of Gomory’s cutting plane algorithm for pure integer programs is numerically stable and capable of solving instances of considerable difficulty for other codes, we decided to run it on Dr. Knuth’s problem, which is shown in Appendix 1. We first ran the problem with IBM ILOG CPLEX 11.2, with default settings. Of course, the CPLEX of today is by far superior to the one of 1995, but it still took 91,287 nodes of the branch-and-bound tree to find a solution of the same value of 22,996 as the one found by Alevras, and to prove its optimality. We then ran the problem with our version of the lexicographic cutting plane algorithm in its simplest variant, which generates a sequence of individual cuts and reestablishes the lexicographic optimality of the tableau after adding each cut. To our pleasant surprise, our algorithm found the same optimal solution of 22,996 after generating 5,618 cuts! The contrast to the over 90,000 search tree nodes that CPLEX needed to accomplish the task is striking. A branch and bound tree with 90,000 nodes implies the solution of as many linear programs. On the other hand, generating 5,618 cuts requires practically no time, the cuts being read off the simplex tableau, and reoptimizing after the addition of each cut requires the solution of one linear program. True, the linear program is solved to lexicographic optimality, which involves considerably more pivots than just finding an optimal solution; but still, the effort is proportional to the number of cuts and simply not comparable to the effort of solving 90,000 linear programs. In particular, the total number of pivots required by our lexicographic method (including those due to lexicographic reoptimization) is surprisingly small: just 31,842 instead of the over 167,000 required by CPLEX. In terms of the time needed by the two codes to solve the problem was comparable (6.4 seconds for CPLEX and 6.2 seconds for our code), but the current version of CPLEX is the result of two decades of development by a strong team of programming experts, an “industrial product”, whereas our code is merely a research tool, the “handicraft” work of a single student helped by his advisor. We also solved Knuth’s problem by two additional versions of our lexicographic cutting plane algorithm, introduced in [1], and for comparison, with a standard implementation of Gomory’s cutting plane algorithm, and with CPLEX 11.2. Table 1 reports the outcome of these runs. The entries of the first column are as follows. – TB stands for the “textbook” implementation of Gomory’s 1958 cutting plane method for pure integer programs, in its modern, “multi-cut” variant, which generates a Gomory cut (or GFC) from every fractional basic variable before reoptimiz-

3 Method TB (multi-cut) lex-GFC (single cut) lex-GFC (multi-cut) Bin-GFC Cplex 11.2 (default) Cplex 11.2 (strong-branching) Cplex 11.2 (aggressive cuts)

Time (sec.s) >8h 6.20 0.30 0.15 6.40 27.70 6.00

#LP’s 2,197 5,618 111 42 51,495 673,076 38,216

#pivots 252,225 31,842 1,941 1,032 167,393 308,437 160,102

#cuts 124,414 5,618 5,524 1,171 19 22 83

#nodes 0 0 0 0 91,287 98,973 43,098

Table 1 Comparison of different solution methods; the first 4 are pure cutting plane methods, while the last 3 are branch-and-bound routines of Cplex under 3 different parameter settings.

–

– –

–

ing. The original, “single-cut” version, which reoptimizes after every cut, breaks down after a few thousand cuts. lex-FGC (single cut) is our lexigraphic cutting plane algorithm using GFC cuts, with (partial) lexicographic reoptimization after every cut. For details see [1], sections 3 and 7. lex-FGC (multi-cut) is the same, but with cuts generated from every fractional basic variable before (partial) lexicographic reoptimization Bin-FGC is a version of the lexicographic algorithm introduced in [1] (section 5), which replaces the objective function with a binary search of its optimal value. Thus lex-GFC (multi-cut) is used repeatedly to explore the existence of a feasible solution within a given cost-interval. CPLEX 11.2 is the integer solver of CPLEX, run under three different parameter settings, namely “default,” “strong branching” and “aggressive cuts.”

According to Table 1, the best performer in terms of both time and number of linear programs solved is Bin-GFC, which is 40 times faster than CPLEX under the best parameter settings. The multi-cut version of lex-GFC is 20 times faster than CPLEX. The growth of the lower bounds given by the LP values with the number of iterations (i.e., number of LP’s solved) is shown for the first three methods of Table 1 in Figure 1. Comparing the data for the best performing version of the lexicographic approach with those for the best CPLEX settings, we find that the former requires less than 1% of the number of pivots required by the latter. Another highly interesting aspect of our experiments is the fact, illustrated in Figure 2, is that the sequence of determinants of the successive basis matrices after each lexicographic reoptimization, instead of increasing as usual when the cuts produce dual degeneracy and the solutions get closer and closer to each other, in fact keeps decreasing throughout the run. As discussed in [5, 1], the main reason for this behavior is that as a sequence of cuts is added to the simplex tableau, dual degeneracy becomes a typical phenomenon, resulting in multiple optimal solutions to the current LP relaxation. When the standard dual simplex method is used to reoptimize after adding a new cut, the algorithm tends to pick among the many optimal solutions one that is closest to the old solution. This leads to the appearance of very large basis determinants, needed to differentiate between points closer and closer to each other. The lexicographic dual simplex on the other hand does not stop the reoptimization upon finding the first dual-feasible solution, but continues to pivot until it reaches a lexicographically dual-feasible solution, which is usually among those furthest from the previous optimum; see Figure 4 for an illustration of the LP optimal solutions. In this way, the determinant of the basis matrix is not forced to increase during the run. Our explanation for the fact that the determinants actually

Gap closed 69% 100% 100% 100% 100% 100% 100%

4

4

2.3

x 10

2.2

TB single lex multi lex

Object function value

2.1 2 1.9 1.8 1.7 1.6 1.5 1.4 0 10

1

10

2

10 LP’s solved

3

10

Fig. 1 Lower bound growth for the pure cutting plane methods.

decreased during the run is that whereas the initial matrix included some relatively large numbers (like 50), the cutting planes generated were numerically “clean”, i.e. with small coefficients (see Figure 3), and as they gradually made the original constraints redundant, they produced basis matrices with smaller numbers and cleaner structure. In his message to us accompanying the description of the problem, Dr. Knuth states that he has recently discovered an error in the original problem statement, which he hopes will not invalidate the optimal solution found by Alevras: one of the inequalities needs to be replaced with two others, also involving a new variable (see the correction at the end of the Appendix, or pages 435–436 of [4]). We also ran the corrected version of the problem, and found that the seemingly minor change makes the problem substantially more difficult: although the optimal solution remains the same as before there is a sharp increase in the computational effort required by either method to solve it. Table 2 reports the results. Although the difference in performance is less sharp than before, CPLEX with the most favorable settings still takes 15 times as long as the lexicographic Bin-GFC.

5

70

10

TB single lex multi lex

60

10

50

Determinant

10

40

10

30

10

20

10

10

10

0

10

0

1

10

10

2

3

10 LP’s solved

10

Fig. 2 Comparison of basis determinants

Method TB (multi-cut) lex-GFC (single cut) lex-GFC (multi-cut) Bin-GFC Cplex 11.2 (default) Cplex 11.2 (strong-branching) Cplex 11.2 (aggressive cuts)

Time (sec.s) >8h 8.1 0.9 0.6 12.5 38.5 9.4

#LP’s 6,407 6,554 316 161 122,316 895,168 66,776

#pivots 691,746 39,528 5,015 5,027 402,009 580,067 257,461

#cuts 422,821 6,554 13,652 1,243 21 22 59

#nodes 0 0 0 0 158,058 133,545 84,574

Table 2 Comparison of different solution methods for the corrected instance.

Appendix (a) Problem statement (1960) 30(50t21 + x2 ) + 10(50(t21 + t26 − t12 ) + x2 ) + (50(t21 + t28 − t14 ) + x2 ) + (50(t21 + t26 − t12 + t28 − t14 ) + x2 ) + 2(50t30 + x2 ) subject to the following constraints: x2 + 2 ≤ 2z1

Gap closed 96% 100% 100% 100% 100% 100% 100%

6

14

10

TB single lex multi lex

12

10

Average cut coefficient value

10

10

8

10

6

10

4

10

2

10

0

10

0

10

1

10

2

10 LP’s solved

3

10

Fig. 3 Maximum cut coefficients

50t2 + x3 − 2z1 − 24 ≥ 0 50(t3 − t2 ) + x4 − 2z2 − 5 ≥ 0 50(t4 − t3 ) − x4 − 6 ≥ 0 50(t5 − t4 ) + x5 − 15 ≥ 0 50(t6 − t5 ) + x6 − 2z3 − 4 ≥ 0 50(t7 − t6 ) + x4 − 2z4 − 12 ≥ 0 50(t8 − t7 ) + x7 − 2z5 − 35 ≥ 0 50(t9 − t8 ) + x8 − 2z6 − 3 ≥ 0 50(t10 − t9 ) + x6 − x8 − 10 ≥ 0 50(t11 − t10 ) + x3 − 2z7 − 8 ≥ 0 50(t12 − t11 ) + x8 − x3 − 3 ≥ 0 50(t13 − t12 ) + x9 − 2z8 − 7 ≥ 0 50(t13 − t12 ) + x9 + w1 − 2z8 − 21 ≥ 0 50(t14 − t13 ) + x5 − 2z9 − 35 ≥ 0 50(t15 − t14 ) + x3 − 2z10 − 34 ≥ 0 50(t16 − t15 ) + x5 − x3 − 6 ≥ 0 50(t17 − t16 ) + x9 − 2z3 − 30 ≥ 0 50(t18 − t17 ) + x4 − 2z9 − 8 ≥ 0 50(t19 − t18 ) + x3 − 2z5 − 34 ≥ 0

x3 + 2 ≤ 2z2

x5 + 3 ≤ 2z3 x6 + 2 ≤ 2z4 x4 ≤ 2z5 + 1 x7 ≤ 2z6 x6 ≤ 2z7 + 1 x8 ≤ 2z8 + 1 x9 + w1 ≤ 2z9 + 1 x5 + 2 ≤ 2z10

7

TB Lex

Lex vs TB solutions trajectory 4

x 10

2.2 2.1

objective

2 1.9 1.8 1.7 1.6 1.5 1.4

Y

X

Fig. 4 LP solution trajectories (X and Y represent a projection of the variable’s space; see [5] for details)

50(t20 − t19 ) − 2z11 − 9 ≥ 0 50(t21 − t20 ) + x2 − 9 ≥ 0 50(t22 − t9 ) + x3 − x8 − 9 ≥ 0 50(t23 − t22 ) + x5 − 2z11 − 8 ≥ 0 50(t24 − t23 ) + x4 − x5 − 6 ≥ 0 50(t25 − t24 ) + x6 − x4 − 6 ≥ 0 50(t26 − t25 ) + x8 − x6 − 3 ≥ 0 50(t27 − t7 ) + x7 − 2z5 − 21 ≥ 0 50(t13 − t12 ) + x9 + w2 − 2z8 − 24 ≥ 0 50(t28 − t27 + t8 − t13 ) + x5 − 2z12 − 35 ≥ 0 50(t29 − t13 ) − 2z9 − 15 ≥ 0 50(t30 − t29 ) + x2 − 3 ≥ 0

x3 ≤ 2z11 + 1

x9 + w2 ≤ 2z12 + 1

in nonnegative integers (x2 , . . . , x9 , t2 , . . . , t30 , z1 , . . . , z12 , w1 , w2 ).

8

Optimal solution of cost 22,996 (nonzero entries) t21 = 10; x2 = 9; t26 = 3; t12 = 1; t28 = 4; t14 = 4; t30 = 5; z1 = 6; x3 = 36; z2 19; x4 = 43; t4 = 1; x5 = 65; z3 = 34; x6 = 72; z4 = 37; t7 = 1; z5 = 21; x7 127; z6 = 64; x8 = 139; t10 = 2; z7 = 36; t11 = 3; x9 = 198; z8 = 69; w1 = 11; z9 104; z10 = 34; t15 = 6; t16 = 6; t17 = 4; t18 = 8; t19 = 9; z11 = 18; t20 = 10; t22 3; t23 = 3; t24 = 4; t25 = 4; w2 = 14; z12 = 106; t29 = 5.

= = = =

(b) Correction (after 1995) Replace inequality 50t18 − 50t17 + x4 − 2z9 ≥ 8 by the pair of inequalities 50t18 − 50t17 + x4 − 2z13 ≥ 8, x9 − 2z13 ≤ 1 involving a new non-negative variable z13 . The optimal solution is the same as for the previous (uncorrected) instance, with z13 = 99.

Acknowledgments The research of the first author was supported by the National Science Foundation through grant #DMI-0352885 and by the Office of Naval Research through contract N00014-03-1-0133. We thank Donald Knuth for calling our attention to this remarkable instance and providing all the information available on its past history.

References 1. E. Balas, M. Fischetti, and A. Zanette. On the enumerative nature of Gomory’s dual cutting plane method. Mathematical Programming, Series B, July 2010. DOI: 10.1007/978-3-54068891-4 29. 2. D. Knuth. Minimizing drum latency time. Journal of the ACM, 8:119–150, April 1961. 3. D. Knuth. An integer programming problem. Manuscript, July 1993. 4. D. Knuth. Selected papers on Design of Algorithms, pages 435–436. CSLI Lecture Notes n. 191, Stanford, California, February 2010. 5. A. Zanette, M. Fischetti, and E. Balas. Lexicography and degeneracy: can a pure cutting plane algorithm work? Mathematical Programming, Series A, December 2009. DOI: 10.1007/s10107-009-0335-0.