EUROPEAN JOURNAL OF OPERATIONAL RESEARCH ELSEVIER
European Journal of Operational Research 94 (1996) 206-211
OR Software - ORSEP Operations Research Software Exchange Program
Edited by Professor H.W. Hamacher
ACCPM
A library for convex optimization based on an analytic center cutting plane method* J. G o n d z i o * , O. du M e r l e , R. S a r k i s s i a n , J.-R Vial
Logilab, HEC, Section of Management Studies, University of Geneva, 102 Bd Carl Vogt, CH-1211 Gen~ve 4, Switzerland
Keywords: Large scale convex optimization; Cutting plane method; Interior point algorithm; Analytic center
1. Introduction
In this note we are concerned with Goffin, Haurie and Vial's [9] Analytic Center Cutting Plane Method (ACCPM for short) for large-scale convex optimization. Its state-of-the-art implementation [ 14] is now
* This research has been supported by the Fonds National de la Recherche Scientifique Suisse, grant #12-34002.92. * Corresponding author. On leave from the Systems Research Institute, Polish Academy of Sciences, Newelska 6, 01-447 Warsaw, Poland.
Hardware information: Any computer with C-F+ and FORTRAN 77 compilers.
Software information: C + + and FORTRAN 77. Public domain software in OR can be contributed to ORSEP by sending a diskette containing the code and the code description to Professor H.W. Hamacher, Department of Mathematics, Uni-
versity of Kaiserslautern, P.O. Box 3049, 67653 Kaiserslautern, Germany, e-mail:
[email protected]. Copies of published codes can be requested from the same address or directly from the authors by sending a formatted diskette and a self-addressed enevelope. See EJOR 48(1) (1990) 161-162, for details.
available upon request for academic research use. Cutting plane methods for convex optimization have a long history that goes back at least to a fundamental paper by Kelley [ 16]. There exist numerous strategies that can be applied to "solve" subsequent relaxed master problems in the cutting planes optimization scheme. In the Analytic Center Cutting Plane Method, subsequent relaxed master problems are not solved to optimality. Instead, an approximate analytic center of the current localization set is looked for. The theoretical development of ACCPM was started by Goffin and Vial [8]. It was later continued in [9,10] and led to the development of the prototype implementation of the method due to du Merle [6] that was successfully applied to solve several nontriviai convex optimization problems [2,3]. A need to solve very large optimization problems [ 11 ] motivated us to rewrite the program from scratch so as to make it more flexible and easier to use for the general convex optimization problems. The new implementation of the method due to Gondzio, du Merle, Sarkissian and Vial [ 14] responds to these needs. All recent theoretical and practical developments [ 7,15 ] have been incorporated into it, resulting in the creation of a sophisticated optimization tool.
0377-2217/96/$15.00 Copyright (~) 1996 Elsevier Science B.V. All rights reserved. Pll S 0 3 7 7 - 2 2 1 7 ( 9 6 ) 0 0 1 6 9 - 5
J. Gondzio et al./European Journal of Operational Research 94 (1996) 206-211
207
We now make it available to the user, in the form of a callable library ]_±baccpm. a.
of optimality cuts at yt, l E Aopt that define a piecewise linear approximation f : R " H R to the convex function f
2. F u n d a m e n t a l s of A C C P M
f ( y ) = ~ x { f ( y l ) + (X t, Y _ yZ)}.
The algorithm implemented in ACCPM combines two powerful optimization techniques: cutting plane methods [ 16] and interior point algorithms [ 13,17]. Its detailed description is definitely beyond the scope of this short note. The reader interested in the theory and the implementation details should consult [ 14] and the references therein. We shall nevertheless recall very briefly some basic concepts used within ACCPM. Although their understanding would surely be useful, it is not a pre-condition for the application of a black box type A C C P M library. We are concerned with the solution of a convex optimization problem
It also generated a set of feasibility cuts at yt, l E Afsb
minimize f ( y ) subject to y E Y,
(1)
where f : ll~m ~ R is a convex function, and Y E R m is a convex set such that Y C {y : f ( y ) < oc}.
Since the problem is convex, both the epigraph of f and the set Y can be approximated by intersections of hyperplanes. The procedure called oracle first tests whether y E Y. If y E Y, then the oracle generates a subgradient X E Of(y) at y with the property
f ( Y ' ) >1 f ( Y ) + ( X , Y ' - Y),
Vy' E Y
(2)
This inequality is a supporting hyperplane for the function to be optimized; we shall call it an optimality cut. If y ¢( Y, then the oracle generates a hyperplane (~r,o-) E R m × N for Y that separates y from Y, i.e., ( ~ , y ) < o-
and
that define an outer approximation Y of the feasible set Y
Y={y:
(7"rt,y) ~>o-l,
VIEAfsb}.
( ~ , y ' ) ~> o-,
Vy' E Y
(3)
(5)
2.2. The relaxed master problem The program minimize f ( y ) subject to y E Y is an outer relaxation of the master problem ( 1 ). Note that Y is a polyhedron and f is a polyhedral (piecewise linear) function. The above program can thus be reformulated as a linear problem, the relaxed master
program minimize
(
subject to ( /> f ( y t ) + (XI, y _ yt}, (7/, y) /> o -t, Vl E Afsb.
2.1. The oracle
(4)
V1 C Aopt, (6)
Its solution gives a lower bound for the master problem ( 1 ). Observe also that the best feasible solution in the generated sequence provides an upper bound for the master problem 0u = min { f ( y t ) } . l E Aoot
2.3. Localization set For a given upper bound O, let us consider the following polyhedral approximation called a localization set
This hyperplane improves the description of the feasible set Y; we shall call it a feasibility cut. Suppose the oracle has been called at a given sequence of points {yl}, l E A. Some of these points were feasible: yt E Y, l E Aopt; the others were not: yt ~ y, l E Afsb = A \ Aopt. The oracle generated a set
z:(0) = { ( ( , y ) : 0~>(,
(>7 f ( y t ) + (Xl, y _ yt), (Trt, y) ~> o t,
Vl E Afsb}.
Vl C Aopt, (7)
208
J. Gondzio et al, /European Journal of Operational Research 94 (1996) 206-211
It is the best (outer) approximation of the optimal set in (1). We are now ready to give a prototype cutting plane method [ 16]. Its main steps are the following: (i) Compute a point ( ( , y ) E /2(0u) and an associated lower bound O. (ii) Call the oracle at ( ( , y). The oracle returns one or several cuts; if all of them are optimality cuts, the oracle also generates an upper bound f ( y ) . (iii) Update the bounds: (a) i f y is feasible, 0u := m i n { f ( y ) , 0 u } ; (b) Ol := max{0_, 01}. (iv) Update the upper bound 0 in the definition of the localization set (7) and add the new cuts. These steps are repeated until a feasible point to ( 1 ) is found such that 0u - 01 falls below a prescribed optimality tolerance.
2.4. Multiple cuts Let us observe that we have been concerned by now with a fairly general convex optimization problem ( 1 ). In many applications, one can benefit from two special structures: the objective is additive,
2.5. Analytic center There exists a number of possible strategies that can be applied at step 1 of the cutting plane method. The optimal point strategy [4], for example, consists in solving every relaxed master program (6) to optimality. The central point strategy [8] consists in finding some "central point" in a localization set. The discussion of their advantages is beyond the scope of this note. Whatever strategy is used, however, we are concerned with a linear optimization problem (6) or, more precisely, with finding some point in the localization set. Localization set (7) is a polytope in IRm+l. It is defined by a set of n inequalities that can be written as
/2 = {y : ATy <. c},
(9)
where a vector (y, ( ) E •m+l is, for brevity, denoted with y and A E R (m+j)xn. We can represent the localization set in an equivalent form that uses explicit dual slack variables s E ]Rn
/2={y:ATy+s=c,
S/>0}.
(lO)
We associate a (dual) potential function to the set of linear inequalities (9)
k
f(Y) = Z
n
fi(Y),
(11)
~D ----"Z In s j . j=l
i=1
and the feasible set is the intersection of different subsets, P
Y--N j=l
The analytic center is a (unique) maximizer of the dual potential over the interior of/2. It is easy to show that it is an interior point (x > 0, s > 0, y) that satisfies Ax = 0,
In such a case, it is possible to modify the relaxed master program so that multiple cuts are introduced one at time: k
minimize i=1 subject to (i >- f i ( Y t) + (XI,Y - yt),
(7r j ,Iy )
Vl E Aopt, i = 1 . . . . . k, >/ o-~, VI E Afsb, j = 1,. . . , p .
(8)
A T y + s = C,
(12)
XSe = e,
where X and S are diagonal matrices with the elements Xj and sj, respectively and e is the n-vector of all ones. In other words, an analytic center is the primal and dual feasible point such that all complementarity products XjSj, j = 1 . . . . . n, are "centered" to one. Note that e~D = e~-'~j=l~Insj =
ft
j=l
sj.
J. Gondzio et al./European Journal of Operational Research 94 (1996) 206-211
Thus the analytic center is a point that maximizes the product of slacks in I]. This gives it a natural geometrical interpretation: analytic center is a point that maximizes the product of distances to all faces of the polytope L. (Note, by the way, that repeating a given constraint in a definition of 12 does not change the geometry of the polytope, but it does change its analytic center - - the repeated constraint moves the analytic center away.) ACCPM employs the potential reduction algorithm of de Ghellinck and Vial [5] to find an analytic center of the current localization set.
3. H o w to use A C C P M
We are aware that ACCPM represents a fairly complicated optimization technique. Hence to facilitate its usage, we provide the user with four carefully designed examples of its application. They start from a simple unconstrained quadratic optimization and, step by step, get more complicated up to a constrained quadratic programming problem with a separable objective and intersection type definition of the feasible set. We hope to have demonstrated with these examples quite a sophisticated usage of the library in a still comprehensive way. Some potential ACCPM users may have already implemented a variant of the cutting plane method (e.g., Dantzig-Wolfe approach, bundle method, ElzingaMoore approach, etc.). For such users, a replacement of the implemented technique with a routine that finds an analytic center might turn out to be a half-an-hour exercise. The inverse is also true: any ACCPM implementation can be easily transformed into another cutting plane approach. Whenever a cutting plane algorithm is thoughtfully designed, such a modification reduces, in practice, to a replacement of the get_AC routine from l i b a c c p m , a with an appropriate g e t _ type routine. To build a new application, the user has to install the l i b a c c p m , a library and follow the examples.
3.1. Installation o f A C C P M
ACCPM library is written in C-F+. An exception are the Cholesky factorization routines written in
209
FORTRAN 77: the sparse one due to Gondzio [ 12] and the dense one from LAPACK [ 1 ]. The library and examples of its application are placed in a directory accpm. This directory contains a r e a d . m e file, a PostScript version of the most recent revision of ACCPM's user's guide ( u s e a c c p m . p s ) and two subdirectories. The first, called a p p l i c a t i o n s , contains four quadratic programming examples QP1, QP2, QP3, and QP4. The second, called ac, contains the ACCPM library, libaccpm, a. The whole directory accpm has been archived (with the Unix command t a r e v f accpm, t a r accpm), resuiting in the creation of a file accpm, t a r . The library is distributed in this form. There are three versions available of the acepm, t a r file: accpm, t a r . ibmxlC is a version compiled with IBM's xlC compiler for IBM RISC 6000 computers and their successors, a c c p m . t a r , ibmg++ is a version compiled with the gnu C + + compiler for IBM RISC 6000 computers, and accpm, t a r . sung++ is a version compiled with the gnu C + ÷ compiler for SUN SPARC workstations. To install the ACCPM library, the authorized user is first of all expected to get the appropriate file a c c p m . t a r for the computer and compiler available. (For example, on an IBM RISC 6000 computer with the xlC compiler, this would require executing the command cp accpm, t a r . ibmxlC accpm, t a r . ) Once the accpm, t a r file has been placed in the current directory, the user should execute the command t a r x v f accpm, t a r accpm. This should create a subdirectory accpm in the current directory. To complete the installation, the user should change to this directory (cd accpm) and run the make command. If this operation is successfully completed, the applications subdirectories QP1, QP2, QP3, and QP4 will contain executable qp programs. It might happen that when linking the program, the loader will not be able to find appropriate libraries l i b l a p a c k , a and l i b b l a s , a. In such a case, the reader will have to install them either in a public place or in the current (application) directory. The public domain FORTRAN source files for both these libraries can be found in n e t l i b (ftp://netlib. art. com). Note that at the beginning of maker i le, several important parameters are defined that the user may want to change. They define the names of compilers: F77 and CXX, the flags (options) for compilation: FFLAGS
210
i
G o n ~ etaL/European ~urnal ~Opera~onalResearch ~ (199~ 206-211
/ * structure types * / struct Problem; struct Work; / * * * * * * ACCPM Callable Library Routines * * * * * * /
extern void extern extern extern extern
void void double void
extern extern extern extern extern extern extern
void void void void void void void
init_problem(Problem**, Work**, int, char*&, double~, double, double, int, int, double*, double*); phaseI_phaseII(Problem*, double*); initnnytime(void); mytime(void); get_AC(Problem*, Work*, char*, int, double, double&, int, int, int*, double*, int*, int*, double*, double*, int*, double*, double*); get~nfo(Problem*, int&, inta, int~, doublea, double~, int&); get_x_primal(Problem*, double*, int&); change_no_of_subpb(Problem*, int); set_compute_binf (int, Problem*); set_testcolumn(int, Problem*); set_push_bound(int, Problem*); set_weight(int, Problem*); Fig. 1. ACCPM include file.
and CFLAGS, and the names (and paths to) different libraries: LAPACK, BLAS, LM, LC, LF77. All these names are exported to m a k e f i l e s in subdirectories. 3.2. Callable routines of the ACCPM library
ACCPM is a library of routines that handle a part of the logic of the cutting plane algorithm. Those routines update the localization set, find an approximate analytic center to it and execute a number of more technical actions. Their callable subset is defined in an ACCPM include file accpm, i n c that should be included in all user's files which make a call to l i b a c c p m , a routines. The include file contains the prototypes of all callable routines. Its form is shown in Fig. 1. 3.3. Running ACCPM
Assume a new user's application has already been compiled and linked. The executable program needs the ACCPM specifications file i n p u t . d a t to be present in the current directory. This file contains a set of technical parameters that control the execution of ACCPM (e.g., determines the requirements imposed on the quality of approximate analytic center, chooses
appropriate data structures to handle cuts, chooses the preferable factorization technique to be applied in the interior point algorithm, and, finally, provides ACCPM with useful hints to determine its initial storage requirements). Below we give an example of such a file chosen appropriately to the quadratic programming applications supplied with ACCPM.
QTEST 1 Verbose 0 sparseG 0 factor_L 0 explicit~ 1 dense_GUB 1 nonzin_cut 2 densityB 2 densityL 0 exp_iters i00 max_iters i00
4. Handling large-scale problems There are no particular limits to the size of problems to be solved, except the total amount of operational
J. Gondzio et al. /European Journal of Operational Research 94 (1996) 206-211
m e m o r y available on a given computer. The origin of the p r o b l e m has clearly an essential influence on these limits. A 64 M B workstation, for example, can h a n d l e really large problems. If cuts are stored as dense c o l u m n s , the rough limits are: m ~< 500 and the total n u m b e r of cuts n ~< 10,000. However, if cuts display a sparse structure, these limits can be pushed m u c h further. The largest p r o b l e m s solved in [ 15] had m = 5 , 0 0 0 and the n u m b e r o f cuts went over n = 100,000. By the appropriate setting o f parameters in input.dat file, the user can help better storage m a n a g e m e n t within A C C P M .
5. Further reading A C C P M is a sophisticated optimization tool. Although its new i m p l e m e n t a t i o n benefits from all experience gathered so far, we are aware that its application is not always straightforward. Potential users of A C C P M should u n d o u b t e d l y consult at least its user's guide, u s e a c c p m , p s . M o r e advanced usage of the library will require careful reading o f the survey by G o n d z i o et al. [ 14].
References l I I Anderson, E., Bai, Z., Bischof, C., Demmel, J., Dongarra, J., Du Croz, J., Greenbaum, A., Hammarling, S., A. McKenney, Ostrouchov, S., and Sorensen, D., LAPACK Users' Guide, SIAM, Philadelphia, 1992. [2] Bahn, O., Goffin, J.-L., Vial, J.-P., and du Merle, O., "Experimental behaviour of an interior point cutting plane algorithm for convex programming: an application to geometric programming", Discrete Applied Mathematics 49 (1994) 3-23. 131 Bahn ,O., du Merle, O., Goffin, J.-L., and Vial, J.-E, "A cutting plane method from analytic centers for stochastic programming", Mathematical Programming 69 (1995) 4573. 14] Dantzig, G.B., and Wolfe, P., "The decomposition algorithm for linear programming", Econometrica 29/4 (1961) 767778. 151 de Ghellinck, G., and Vial, J.-P., "A polynomial Newton method for linear programming", Algorithmica 1 (1986) 425-453.
211
[61 du Merle, O., "Interior points and cutting palnes: a development and implementation of methods for convex optimization and large scale structured linear programming", Ph.D Thesis, Department of Management Studies, University of Geneva, Geneva, Switzerland, 1995 (in French). 17] du Mede, O., Goffin, J.-L., and Vial, J.-P, "On the comparative behavior of Kelley's Cutting Plane Method and the Analytic Center Cutting Plane Method", Technical Report 1996.4, Department of Management Studies, University of Geneva, Switzerland, March 1996. [8] Goffin, J.-L., and Vial, J.-E, "Cutting planes and column generation techniques with the projective algorithm", Journal of Optimization Theory and Applications 65 (1990) 409429. [9] Goffin, J.-L., Haurie, A., and Vial, J.-P., "'Decomposition and nondifferentiable optimization with the projective algorithm", Management Science 38/2 (1992) 284-302. [ 10l Goffin, J.-L., Haurie, A., Vial, J.-P., and Zhu, D.L., "'Using central prices in the decomposition of linear programs", European Journal of Operational Research 64 (1993) 393409. [11] Goffin, J.-L., Gondzio, J., Sarkissian, R., and Vial, J.-P., "Solving nonlinear multicommodity flow problems by the analytic center cutting plane method", Technical Report 1994.21, Department of Management Studies, University of Geneva, September 1994, revised November 1995, Mathematical Programming (to appear). [12] Gondzio, J., "Implementing Cholesky factorization for interior point methods of linear programming", Optimization 27 (1993) 121-140. [13] Gondzio, J., and Terlaky, T., "A computational view of interior point methods for linear programming", in: J.E. Beasley (ed.), Advances in Linear and Integer Programming, Oxford University Press, Oxford, England, 1996, 106-147. [ 141 Gondzio, J., du Merle, O., Sarkissian, R., and Vial, J.-P., "An advanced implementation of the analytic center cutting plane method and its applications", Technical Report, Department of Management Studies, University of Geneva, December 1995 (in preparation). [ 15] Gondzio, J., Sarkissian, R., and Vial, J.-P., "Using an interior point method for the master problem in a decomposition approach", Technical Report 1995.30, Department of Management Studies, University of Geneva, Switzerland, October 1995, revised in May 1996, to appear in European Journal of Operational Research.
[ 161 Kelley, J.E., "The cutting plane method for solving convex programs", Journal of the SIAM 8 (1960) 703-712. [ 17] Roos, C., and Vial, J.-P., "Interior point methods", in: J.E. Beasley (ed.), Advances in Linear and Integer Programming, Chapter 3, Oxford University Press, Oxford, England, 1996.