Q U A D R A T I C P R O G R A M M I N G : Q U A N T I T A T I V E A N A L Y S I S A N D P O L Y N O M I A L R U N N I N G T I M E A L G O R I T H M S By J A D R A N K A S K O R I N - K A P O V B. Sc. The University of Zagreb, Yugoslavia, 1977 M . Sc. The University of Zagreb, Yugoslavia, 1983 A THESIS S U B M I T T E D IN PARTIAL F U L F I L L M E N T O F T H E R E Q U I R E M E N T S FOR T H E D E G R E E OF D O C T O R O F PHILOSOPHY in T H E F A C U L T Y OF G R A D U A T E STUDIES (Department of Commerce and Business Administration) We accept this thesis as conforming to the required standard T H E UNIVERSITY OF BRITISH C OLU MBIA August 1987 © Jadranka Skorin-Kapov , 1987 In presenting t h i s thesis i n p a r t i a l f u l f i l m e n t of the requirements for an advanced degree at the University of B r i t i s h Columbia, I agree that the Library s h a l l make i t f r e e l y available for reference and study. I further agree that permission for extensive copying of t h i s thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It i s understood that copying or publication of t h i s thesis for f i n a n c i a l gain s h a l l not be allowed without my written permission. Department of The University of B r i t i s h Columbia 1956 Main Mall Vancouver, Canada V6T 1Y3 Date DE-6 (3/81) Abstract Many problems in economics, statistics and numerical analysis can be formu-lated as the optimization of a convex quadratic function over a polyhedral set. A polynomial algorithm for solving convex quadratic programming problems was first developed by Kozlov at al. (1979). Tardos (1986) was the first to present a poly-nomial algorithm for solving linear programming problems in which the number of arithmetic steps depends only on the size of the numbers in the constraint matrix and is independent of the size of the numbers in the right hand side and the cost coefficients. In the first part of the thesis we extended Tardos' results to strictly convex quadratic programming of the form max {cTx-\xTDx : Ax 0} with D being symmetric positive definite matrix. In our algorithm the number of arithmetic steps is independent of c and b but depends on the size of the entries of the matrices A and D . Another part of the thesis is concerned with proximity and sensitivity of integer and mixed-integer quadratic programs. We have shown that for any optimal solution z for a given separable quadratic integer programming problem there exist an op-timal solution x for its continuous relaxation such that \\z — X^QQ < nA(A) where n is the number of variables and A(A) is the largest absolute sub-determinant of the integer constraint matrix A . We have further shown that for any feasible solution z , which is not optimal for the separable quadratic integer programming problem, there exists a feasible solution z having greater objective function value and with \\z — 2||oo < nA(A). Under some additional assumptions the distance between a pair of optimal solutions to the integer quadratic program-ming problem with right hand side vectors b and b', respectively, depends linearly on ||6 — 6'||i. The extension to the mixed-integer nonseparable quadratic ii case is also given. Some sensitivity analysis results for nonlinear integer programming problems are given. We assume that the nonlinear 0 — 1 problem was solved by implicit enumeration and that some small changes have been made in the right hand side or objective function coefficients. We then established what additional information to keep in the implicit enumeration tree, when solving the original problem, in order to provide us with bounds on the optimal value of a perturbed problem. Also, suppose that after solving the original problem to optimality the problem was enlarged by introducing a new 0 — 1 variable, say xn+i . We determined a lower bound on the added objective function coefficients for which the new integer variable zn-|_i remains at zero level in the optimal solution for the modified integer nonlinear program. We discuss the extensions to the mixed-integer case as well as to the case when integer variables are not resticted to be 0 or 1 . The computational results for an example with quadratic objective function, linear constraints and 0—1 variables are provided. Finally, we have shown how to replace the objective function of a quadratic pro-gram with 0—1 variables ( by an integer objective function whose size is polynomially bounded by the number of variables) without changing the set of optimal solutions. This was done by making use of the algorithm given by Frank and Tardos (1985) which in turn uses the simultaneous approximation algorithm of Lenstra, Lenstra and Lovasz (1982). in T a b l e o f C o n t e n t s Abstract ii List of Figures vi List of Tables vii Acknowledgement viii I. Introduction 1 1.1. Preliminary Definitions from Linear Algebra and Convexity Theory . . 2 1.2. Duality in Convex Nonlinear Programming 8 1.3. Integer Quadratic Programming 11 1.4. Some Transformations of Quadratic Programs 12 1.5. A Review of Polynomial Algorithms for Quadratic Programs . . . . 14 1.6. Basis Reduction Algorithm and Application to the Simultaneous Diophantine Approximation Problem 19 II. Towards a Strongly Polynomial Algor i thm for Strictly Convex Quadratic Programs 23 2.1. Setup of the Problem 24 2.2. The Quadratic Programming Algorithm 27 2.3. Extension to the Inequality Constrained Case • . 37 III. Proximity and Sensitivity Results for Quadratic Integer Programming 41 iv 3.1. Proximity Results for Separable Quadratic Integer Programming . . 42 3.2. Sensitivity for Quadratic Integer Programming : The Right Hand Side Case 49 3.3. Extension to Nonseparable Quadratic Mixed-Integer Programming 53 IV. Nonlinear Integer Programs: Sensitivity Analysis for Branch and Bound 57 4.1. The Pure Nonlinear 0-1 Program 58 4.2. The Mixed-Integer Case 62 4.3. The Quadratic 0-1 Problem: An Example 63 V . Simultaneous Approximation in Quadratic 0-1 Programming . 70 5.1. Setup of the Problem 72 5.2. Simultaneous Approximation of Objective Function Coefficients . . . 74 VI. Areas for Further Research 79 Bibliography 81 V L i s t o f F i g u r e s Figure 4.1 Branch and Bound Tree for Example 4 L i s t o f T a b l e s Table 4.1 Sensitivity Analysis Sample for Example 4.1 69 vii Acknowledgement I would like to express my sincere appreciation to the University of British Columbia and especially to the Faculty of Commerce and Business Administration for providing me with a financial support through : Dean Earle D. McPhee Memorial Fel-lowship in Commerce and Business Administration, University of British Columbia Summer Graduate Fellowship, University of British Columbia Graduate Fellowship and Graduate Travel Grant. My sincere thanks go to the Division of Management Science for giving me the opportunity to be one of their Ph.D. Students. I am deeply indebted to Professor Frieda Granot, my supervisor, for her guid-ance, many valuable comments and suggestions and for her friendliness and financial support. Furthermore, I would like to thank other members of my dissertation com-mittee Dr Richard Anstee, Dr Daniel Granot and Dr Maurice Queyranne for their helpful suggestions. Finally, special thanks to my husband Darko for " sailing in the same boat " , and to our two little beauties Lea and Nina for sharing their mama with many hours of overtime work. V U l To the memory of my mother, Fumica Boljuncic, and to my father, Ivan Boljuntic . Chapter I I N T R O D U C T I O N Many problems in economics, statistics and numerical analysis can be formulated as the optimization of a convex quadratic function over a polyhedral set. Moreover, some algorithms for solving large scale mathematical programming problems mini-mize a quadratic function over a polyhedral set as a subroutine, e.g. Held at al. [27] , Kennington and Shalaby [29] . Several methods that are based on solving a quadratic programming subproblem to determine a direction of search were also suggested for optimization problems with nonlinear constraints, e.g. Biggs [5] , Garcia and Man-gasarian [20] and Gill at al. [21]. The existence of efficient quadratic programming algorithms and the fact that nonlinear functions can be sometimes accurately approx-imated by quadratic functions led to the development of approximation methods that make use of quadratic subproblems, e.g. Fletcher [18]. The above mentioned are just some of the reasons why the quadratic programming arose as a very important part of the rich theory of Mathematical Programming. We start by presenting in Section 1.1. some preliminary definitions from linear algebra and convexity theory to be used in this thesis. Many of the results in this thesis make use of duality. We will review the convex nonlinear programming prob-lem and its dual as stated by Wolfe[49] and, as a special case in Section 1.2., a convex quadratic programming problem and its dual as stated by Dorn [13]. In Section 1.3. l we introduce integer quadratic programming problems. In Section 1.4. some trans-formations of quadratic programs which will be used in subsequent Chapters are summarized. We do not attempt to survey the algorithms suggested for optimiz-ing a convex quadratic function, rather we restrict our attention to polynomially bounded algorithms and review them in Section 1.5. An introduction to lattices and the transformation of simultaneous Diophantine approximation problem to a short lattice vector problem will be given in Section 1.6. 1.1. P R E L I M I N A R Y D E F I N I T I O N S F R O M L I N E A R A L G E B R A A N D In this Section we first review some well known definitions and results the details on which can be found in any text-book of linear algebra. We are considering in this thesis the vector space Rn . The elements of Rn are The elements of R are called scalars and will be denoted by lowercase Greek letters a, /?, 7, etc. The function N(x) : Rn —•> R is called a norm if 1) N[x) > 0 for all x e Rn , N{x) =0 if and only if x = 0 ; 2) N(ax) = \a\N{x) for all x e Rn , a e R ; 3) N(x + y)< N{x) + N{y) for all x , y e Rn . We will use the standard notation || x || to denote the norm of a vector x . For 1 < p < oo the p-norm will be given by C O N V E X I T Y T H E O R Y ordered tuples of real numbers denoted as and refered to as vectors. n xeR •n (1.1) 2 In this thesis we will make use of - ]Cr=i \xi\ t n e ^i-norm , IMh = \Xi\2)^ t n e ^2-norm , and |x||oo = max {|xj| : i = l,...n} the /oo-norm . These three norms are equivalent in the sense that for any x e Rn , IMIoo < \\xh < Vn\\x\\ (1-2) IMloo < ll^lll < ^Iklloo (1-3) N| 2 < ll^lli < V^Nh- (1-4) For a scalar a, \a] (resp., [aj ) will denote the smallest (resp., greatest) integer not smaller (resp., not greater) than a . The constraint sets of optimization problems to be treated in this thesis are sets of linear equalities and (or) inequalities which are usually stated in matrix form, therefore the second algebraic object of our interest are real matrices. Unless other-wise stated, in our thesis vectors are considered as column vectors, i.e. x e Rn is an n x 1 matrix. Transposition applied to a vector x (resp., matrix A ) will have the usual notation xT (resp., AT ). For our purposes , some special square matrices deserve to be mentioned here. 1) Identity matrix / : / = (<5iy)"y=1 , Sl}- = 2) Diagonal matrix D : D = (dtj)"y= 1 , dij = 0 for i ^ j. 3) Nonsingular matrix A : A is nonsingular if and only if there exist a matrix B such that AB = BA = I . B is usually denoted as A-1 and is refered to as the inverse of A . 3 J 1 for i = j \ 0 for i ^ j . 4) Symmetric matrix A : A is symmetric if and only if A = AT . 5) Positive Semidefinite matrix A : A is positive semidefinite if and only if xTAx > 0 V x e Rn . A is positive definite if in addition xT Ax = 0 implies x = 0 . 6) Idempotent matrix P : P is idempotent if and only if P2 = P . The vectors a^, ...,xr are said to be linearly independent if and only if their linear combination vanishes in a trivial way only , i.e. OL\X\ + ... + arxr = 0 implies ai = ... = aT — 0 . Otherwise, the vectors are said to be linearly depen-dent and each vector x; with ^ 0 can be expressed as a linear combination of the remaining vectors. A row ( resp., column ) rank of an n x n matrix A is the number of its linearly independent rows ( resp., columns ). The row and column rank of a given matrix always coincide and will be denoted by r(A) to be refered to as the rank of the matrix. A nonsingular n x n matrix A has a full rank , i.e. r[A) = n . A determinant is a function that assigns to each n x n matrix A with columns A\,...,An a scalar value denoted by det A that hasthe following properties: For each scalar a and each i = l,...,n 1) det(Ai,...,aAi,...,An) = a d e t ( A i , A i , A n ) , 2) det(Ai,.., A{,.., Aj,.., An) = det(Ax,..,Ai + aAj,.., Aj,.., An) for each j ^ i , 3) det(I) = 1 . As a consequence, it can be observed that det A for a singular matrix A is equal to zero. For a nonsingular matrix with all entries integral , det A is not less than one. 4 For n > 3 , determinants are very inefficient as a computational tool, but they are useful to obtain some theoretical results. For example, if A is an n x n nonsingular matrix then the cofactor of any element ars of A is defined as fars = {-l)r+sdetArs, (1.5) co where Ars is the (n — 1) x (n — 1) matrix obtained from A by deleting row r and column s. The inverse matrix A~1 can then be stated as A - 1 = —^—cornA (1.6) de*A v ' where com A is the transpose of the matrix of cofactors of A . The unique solution xT — (x l 5 ...,xn) to a linear system Ax — b where b is an nxl vector is given by detAi . . . X i = -teZ ' t = 1'-'n ( L 7 ) where A{ is the n x n matrix obtained by replacing the i-th column of A by the vector 6 (Cramer's rule ). An upper bound on the value of detA is given by the following inequality \detA\ < ||oi||2---||°n||2- (1-8) where a\,...,an are columns of A (Hadamard's inequality). Combining (1.2) with (1.8) implies that if A = (aiy) with \a.ij\ < a for all i , j , then \detA\ < anni. For an m x n matrix A we will denote by A(A) = max { \detH\ : all square submatrices H of A}. 5 The scalar A(A) , where A is the constraint matrix of a quadratic optimization problem, will play an important role in many results of this thesis. For an m x n matrix A , the matrix norms to be considered in this thesis are m \\A\\i = max > |at-y| , (1-9) 1 R U {+oo, —oo} is convex if i) {x : f{x) = -oo} = 0 , ii) {x : f{x) < + 0 0 } ^ 0 , iii) f{Xx + (1 - A)y) < A/(x) + (1 - X)f[y) for 0 < A < 1 and x, y e Rn . E X A M P L E A quadratic function f(x) = cTx + ^ xTDx with D symmetric positive semidefinite is a convex function. This since {x-y)TD{x-y) >0 for all x, y and therefore xTDx + yTDy > 2xTDy which in turn implies /(Ax + ( l - A ) y ) = c T(Ax + (1 - A)y) + f (Ax + (1 - A)y)T7J(Ax + (l - A)y) = Ac r x + (1 - A)cTy + \\2xTDx + 1(1 - \)2yTDy + A(l - \)xT Dy 6 < XcTx + (1 - \)cTy + \X2xTDx + \(l - X)2yTDy + |A(1 - X)xTDx + \X(l - X)yTDy = XcTx + (1 - X)cTy + \XxTDx + |(1 - X)yTDy = Xf(x) + (1 - X)f(y) . If D is positive definite, then f[x) = cT x-\-\xT Dx is a strictly convex function. With regard to the optimization of a convex function, recall that every local minimum is a global one and that a strictly convex function has a unique minimum. If / is convex, then — / is concave. A set S C Rn is said to be a convex set if it contains all convex combinations of its elements, i.e. for any x 1 , x 2 , x s e S , Yli=i wix% e & f ° r a u Wi > 0 , V i and X2i=i wi ~ 1 • If S contains all nonnegative linear combinations of its elements , i.e. Yli-i W{Xl , tu,: > 0 for all i , then it is called a convex cone . The solution set P of a finite system of linear inequalities Ax < is a convex set refered to as a polyhedron . A cone which is a polyhedron is called a polyhedral cone and it is the solution set of some homogeneous system of linear inequalities Ax < 0 . In Chapter III we will use the following theorem due to Minkowski (which is a special case of Caratheodory's theorem see e.g. [44], page 55): T H E O R E M 1.1.1 Every polyhedral cone has a finite set of generators. P R O O F See Stoerand Witzgall [44], Theorem 2.8.6., page 55. • 7 1.2. D U A L I T Y I N C O N V E X N O N L I N E A R P R O G R A M M I N G Although our thesis is concerned mainly with quadratic programming, the results presented in Chapter IV were shown to be valid for a broader class of problems, namely for convex nonlinear programs in which some variables are restricted to be integral. In this Section we will review a duality for convex nonlinear problems as introduced by Wolfe [49] and will then state a duality theorem of Dorn [13] for a special case of quadratic convex programming problems. Consider the following nonlinear programming problem min {f{x) : gi(x) > 6,-, i = 1 , m } (l-ll) where x is an n x 1 vector, / is (resp., gi ,i = l , . . . , m are) real valued, differentiable convex (resp., concave) functions on Rn . Denote by S a set of feasible solutions to (l.11). In the sequel we assume that on the boundary of the constraint set no singularities will occur, i.e. that some constraint qualification is satisfied. (For detailed discussion on constraint qualifications see e.g. [3] or [36].) For example, we can assume Abadie's constraint qualification (see [4]) of the following form to be valid. Let x e S , then every z satisfying Vgl(x)Tz > 0 for all i such that gi{x) = b{ has to be an element of a cone of tangents T = { z : z = lim Xn(xn - x), xn e S, Xn > 0 for all n, lim xn = x }. n—>oo n—>oo If the constraints are all linear, this constraint qualification is automatically satisfied. (See Bazaraa [4], Lemma 5.1.4., page 164.) The Karush-Kuhn-Tucker necessary optimality conditions (which are also sufficient optimality conditions under suitable convexity assumptions) are as follows. 8 If x is an optimal solution to problem (1.11) under some constraint qualification condition, then there exists a vector u > 0 such that V/(x) = uTVg(x) and uT{g{x) - b) = 0 . For problem (1.11) Wolfe's dual [49] is given by max bTu + f(x) — uTg(x) s.t. uTVg(x) = V/(x) (1.12) u > 0 where bT = (bu 6m) and ff(x)T = (gi{x), ...,gm(x)). If all the constraints are linear, then the objective function of (1.12) can be equivalently written as bTu + f(x) — xTV f(x) . Consider now the convex quadratic programming problem min {cTx + \^xTDx : Ax > b,x > 0} (1-13) where D is a symmetric positive semidefinite n x n matrix , c and x are re x 1 vectors , b is an m x 1 vector and A is an m x re matrix. Positive semidefiniteness of the matrix D implies convexity of the objective function. As stated by Dorn [13] a dual of problem (1.13) can be written as max {bTu - ~xTDx : ATu - Dx < c,u > 0} . (1.14) The Karush-Kuhn-Tucker optimality conditions for a pair of problems (1.13) and 9 (1.14) are the primal and dual feasibility conditions Ax > b x>0 (1.15) ATu — Dx < c u > 0 , and the complementary slackness conditions xT(c - ATu + Dx) = 0 uT{Ax - b) = 0 . (1.16) The existence theorem for quadratic programming states that the feasibility of both the primal and dual programs implies the existence of optimal solutions for each of them . The following theorem is taken from Dorn [13]. T H E O R E M 1.2.1 i) If x = x is a solution to (1.13) , then a solution [u,x) = (u0,x0) exists to problem (1-14) such that Dx = Dx0 . ii) Conversely , if a solution {u,x) = (uo,Xo) to problem (1.14) exists, then a solution x which satisfies Dx = DXQ exists to problem (1.13). In either case the objective function values for (1.13) and (1.14) are equal. Also, if one of the problems (1.13) or (1-14) is feasible while the other is not, then on its constraint set the objective function of the feasible program is unbounded in P R O O F See Dorn [13], page 156. • 10 the direction of extremization (see [13] ). The Fundamental Theorem of linear programming states that if a linear program has an optimal solution, then it has one which is a basic solution of a linear system of constraints. For a quadratic programming problem this is, however, not the case. An optimal solution for a quadratic programming problem may occur everywhere in the feasible region, in the interior as well as on the boundary. Consideration of nonbasic solutions makes quadratic programming more difficult than the linear one. However, if a quadratic program has an optimal pair (x T, uT) of primal and dual solutions satisfying (1.15) and (1.16)), then it has one which is a basic solution for the system of linear equalities and inequalities (1.15) or, equivalently, a solution that is a vertex of a polyhedron defined by (1.15). Combined with Cramer's rule (see (1.7)) this fact gives us a way to bound the values of the primal and dual variables. This will be discussed in more detail in Section 1.5. 1.3. I N T E G E R Q U A D R A T I C P R O G R A M M I N G Many real world problems require a mathematical programming formulation in which all or some of the variables are restricted to be integral. Moreover, a quadratic objective function enables one to take into account the interactions between variables. The applications in, for example, finance [34], capital budgeting [31] or scheduling [40] have natural representations as 0—1 quadratic programming (i.e. integer quadratic programming in which the variables are restricted to be zero or one). In this thesis we will consider a general mixed-integer quadratic programming problem and will discuss some sensitivity aspects of it (see Chapter III and IV) as well as a transformation of the objective function coefficients in Chapter V . 11 1.4. S O M E T R A N S F O R M A T I O N S O F Q U A D R A T I C P R O G R A M S In some cases the form of the objective function cTx + x Dx of a quadratic programming problem is not suitable for our purposes. In this event some transfor-mations are performed to obtain an equivalent quadratic programming problem of suitable form. In this section we will list the transformations of the objective function to be used in subsequent chapters. Consider, for example, the quadratic cost matrix D . Without loss of generality one can assume that D is symmetric since, if not, D = ^(D + DT) is symmetric and replacing D by D in a quadratic programming problem of the form min {cTx + xTDx : Ax > b, x > 0} will not change the objective function value. Suppose, next, that we have a quadratic 0 — 1 minimization problem of the form min f(x)=cTx + xTDx s.t Ax > b (1.17) 0 < x < 1 x integer If we want to solve this problem by implicit enumeration where at each node the continuous relaxation of a corresponding integer subproblem is solved, then we would like to ensure the convexity of the objective function (in order to avoid local minimum points). If D is not positive semidefinite, it is shown in [25] that problem (1.17) can be replaced by an equivalent problem in which the objective function is given by f'(x) = {cT-XeT)x + xT{D + XI)x (1.18) where eT = (l , . . . , l) and A is a positive scalar such that D + XI is positive 12 semidefinite. This is due to the fact that x2 = x for any vector x of zeros and ones. It is often desirable to have a homogeneous quadratic form in the objective func-tion of a quadratic programming problem. This can be achieved by adding a new variable and a new constraint. For example, for problem (1.13) an equivalent homo-geneous problem is min \{xTM(°TCQ)Ca) s.t. Ax > 0 (1.19) a = 1 x > 0 . Note that if D is positive definite adding a constant ~cTD~1c to the objective function results with a convex quadratic program since the matrix cTr)-ic) is positive definite ( observe that ( X T , O J ) ( ^ CTD-1 c)i.'a) c a n ^ e w r n ^ e n as (x + ac7D-l)T'D{x + ac1'D~l) ). An alternative way to homogenize the objective function of a strictly convex quadratic programming problems is to leave the matrix D unchanged, but to change the right hand side vector. This by replacing the vector x in (1.13) by z — D~lc . The resulting problem is ~-cTD~1c + min -zTDz 2 2 s.t. Az >b + AD'h (1.20) Iz- Ix = D'lc x > 0 . 13 We will employ this transformation in Chapter II. Due to the fact that a positive semidefinite matrix can be diagonalized, we shall see in Chapter III how a convex quadratic programming problem can be replaced by an equivalent separable convex quadratic programming problem. Without loss of generality assume that a positive semidefinite matrix C is of the form ( Q ° ) where A = ( 0 . Observe that the length L' of the input data for P' is not greater than 2L . Using Lemma 1.5.2, one can then give an upper bound to the components of (xT,yT) . Namely, any component of (xT,yT) will have the form | where t and s are integers such that \t\ , \s\ < 22L . Since the objective function is quadratic, the smallest rational number will be ^L)2 ~ 2 ^ a n c * n e n c e \dTx+ -xTDx\ < 25L . I 2 l -After checking the compatibility of a system of linear equalities and inequalities (1.27), Kozlov at al. found the exact value /o = | of the objective function. This was done by checking the compatibility of 13L + 2 systems P^ of the form Ax < b f[x) < . where tk and Sk are integers such that \tk\ < 2 1 3 L + 2 and \sk\ < 2 8 i + 2 . Having 17 the objective function value of (1.25), an exact optimal point was obtained. Helgason, Kennington and Lall [28] gave a polynomially bounded algorithm for a class of singly constrained quadratic programming problems of the form min -xTDx — aTx 2 n s.t. YZXJ = c l 1 - 2 8 ) j=i 0 i = l,...,n (1.30) j=i as where Vij are any integers with det(vij) =_ 1 . Then each a t can be written X2y=i wija'j with integral tu,-y . Substituting this expression into (1.30) and using the linear independence of a; , it follows that Y^WMi = {l 1 if i = I , otherwise . Hence, det(wij)det(vij) = 1 implying det(u>ij) = det(vij) =t. 1 • It follows then that det(ai,-...,an) =1 det{a\,...,a'n) . (See for example Cassels [7]). We can, therefore, define the d e t e r m i n a n t of a l a t t i c e , detA = |det(A)| where A = (ai,...,an) is any basis of the lattice. Geometrically, detA denotes the volume of a parallelopiped whose vertices are lattice points and which contains no other lattice point. An upper bound of detA is given by Hadamard's inequality detA < 1 «x 1 i 2 -'" [| 112 A lower bound of detA is due to Hermite (1850) and is given as follows. 20 Every n—dimensional lattice A has a basis ( 6 1 , . . . , b n ) such that IM|2---||M2 < cndetk (1.31) where cn is a constant depending only on n . This result led to the question of finding such a basis. As stated by Minkowski, the solution to the above problem always exists provided cn > ( f^) 2 • To find a basis in the lattice for which the product of the euclidean norms of its vectors is n ( n — 1 ) minimal is NP—hard. However, for a weaker bound, taking c n = 2 4 , Lovasz gave a polynomial algorithm to find a basis satisfying (1.31). (See Lenstra at al. [33]). A related problem is the following (Short Lattice Vector Problem) : Given an n—dimensional lattice A and a number A , find a vector b e A , b 7^ 0 such that H&H2 < A . A classical result of Minkowski implies that for A > 2^J-^-e y/detk such a vector always exists, but no polynomial algorithm to find it is known to date. The shortest vector in the reduced basis obtained by Lovasz's basis reduction algorithm has a length at most 2~*~~ \fdetk . This is not the shortest vector in the lattice, however it is very useful in some applications. Consider, for example, the Simultaneous Diophantine Approximation Problem (see e.g. Lovasz [35]) : Given cci , . . . ,a n e Q , 0 < e < 1 and Q > 0 find integers p i , . . , p n and q such that 0 < q < Q , | a » - - | < - ^ , t = l , . . . , n . (1.32) 21 A classical result due to Dirichlet (1842) is the proof of the existence of integers q and i = 1, ...n such that (1.32) is satisfied, whenever Q > s~n . No polynomial algorithm to find these integers is known (except for the case n = 1 when the method of continued fractions can be applied). However, a weaker approximation IC /"> n.( re+ 1 ) _ (for Q > 2 4 e n ) can be found by transforming this problem into a short lattice vector problem and using Lovasz's basis reduction algorithm as follows [35] : Consider the lattice A(A) generated by the columns of the (n + 1) x (n + 1) nonsingular matrix / l 0 Vo 0 a i \ 1 ocn Q Any vector b eA(A) has the form bT = (px + p n + 1 o : i , . . . ,p n + p n + i a n , p n + 1 - | ) , where pT = (p l 5 . . . , p n + 1 ) e . Suppose that 6 ^ 0 but ||6||2 < e . This implies Pn+\ 7^ 0 . Without loss of generality assume p n +i < 0 and denote q = —pn+i . It follows then that \bi\ = \pi — qoti\ < E , i = l,...,n E >n+l | q < E or q < 0. (1.33) (1.34) The shortest vector in the reduced basis of a lattice A(/l) obtained by Lovasz's algo-rithm satisfies ||6||2 < 2* " + 0} with D being a positive definite matrix. We assume, without loss of generality, that A and D are integral. We develop a polynomially bounded algorithm for solving the strictly convex quadratic problem where the number of arithmetic steps is independent of c and b but depends on the size of the entries of the matrices A and D. If in partic-ular the size of the entries in A and D is polynomially bounded in the dimension of the input, the algorithm is strongly polynomial, e.g., when the quadratic term corresponds to a least squares and A is a node arc incidence matrix of a directed graph. Following Tardos [45] the algorithm presented here finds optimal primal and dual solutions to the quadratic programming problem (if they exist) by solving a sequence of simple quadratic programming problems using the polynomial algorithm 23 for solving quadratic programming problems given in [30] and by checking feasibility of a linear system in time independent of the right hand side using Tardos' feasibility algorithm [45]. 2.1. S E T U P O F T H E P R O B L E M For simplicity of exposition we will first consider the quadratic programming problem of the form max { cTx — -xTDx : Ax = fe, x > 0 } (2.1) where A is an integral m x n matrix with rank (yl) = ra; D is an integral n x n symmetric positive definite matrix; c and x are n-vectors and fe is an m-vector. The validity of the algorithm for quadratic programming problems with inequality constraints will be discussed in Section 2.3. Using the fact that D is nonsingular we can substitute z ~ x — D~lc in (2.1) resulting with the following equivalent problem ^cTD-lc + max {~zTDz : Az = b - AD~xc , lz- Ix= -D~xc , x > 0}. (2.2) Recall that positive definiteness of D implies that the objective function of (2.1) is strictly concave which in turn implies the uniqueness of the optimal solution of (2.1) (if one exists). Moreover, if the set {x : Ax — fe, x > 0} is not empty, (2.1) will be bounded. The uniqueness of the optimal value of x implies also the uniqueness of 24 the optimal value of z; however observe that z is not restricted to be non-negative any longer. Using now Dorn's duality [13], we can state a dual of the maximization problem given in (2.2) as min {yT{b - AD~lc) - vTD~1c + -zTDz : ATy + Iv + Dz = 0, t; < 0}. (2.3) 2 Substituting v = — Dz — ATy results with (see Section 1.4) min {yT{b - AD~lc) + [yTA + zTD)D'1c + -zTDz : ATy + Dz>d} = 2 min {yTb - yTAD~lc + yTAD~lc + zTc + -zTDz : ATy + Dz > 0}. 2 Finally, after adding slack variables 5 , we get a dual of (2.2) of the form -cTD~1c + min {cTz + bTy + -zTDz : ATy + Dz - Is = 0, 5 > 0}. (2.4) 2 2 It is easy to see that replacing z by x — D~lc will give us the following dual of (2.1) min {bTy + ^xTDx : ATy + Dx - Is = c, s > 0}. (2.5) It is important to note that the same slack variables appear both in (2.4) and (2.5). Since the algorithm to be described in the following is significantly simpler when applied to problems with zero right hand side, we will always use the above transformation to replace a pair of primal and dual problems of the form (2.1) and (2.5) by an equivalent pair of the form (2.2) and (2.4). Recall that the Karush-Kuhn-Tucker optimality conditions for a pair of problems (2.1) and (2.5) will have the form 25 (2.6) x,s > 0 and x s = 0 . At each iteration of the algorithm to follow we will detect at least one new slack variable which is equal to zero in all optimal pairs of primal and dual solutions for the pair of problems (2.1) and (2.5). Or equivalently, at each iteration we will detect at least one new dual constraint which is tight at optimality. After each iteration we will add constraints of the form Si = 0 , i E I (where J is the set of slack variables detected to be zero at the current iteration) to the linear system given in (2.6) and perform Tardos' feasibility algorithm which will give us a feasible solution (xT,yT,sT). We will check if xTs = 0. If so, the algorithm will terminate since an optimal pair of primal and dual solutions was determined. If on the other hand xTs 0, we will perform another iteration of our algorithm. In at most n iterations we will find a pair of optimal primal and dual solutions. As stated above the algorithm will be applied to a minimization problem of the form Before stating the algorithm we will give two preliminary lemmas. The first one is a direct generalization of Lemma 0.1 in [45], while the second is a special case of Lemma 2, p. 707 in [15]. L E M M A 2.1.1 Replacing ( c T , 6 T , 0 r ) in (2.7) by {c'T,b,T,a'T) = (cT ,bT ,0T) P R O O F If (zT,yT,sT) solves (2.7), then it also solves the problem obtained from (2.7) by replacing the linear cost coefficient (c T ,6 T ,0 T ) by (c'T, b'T, a'T) and 26 mm {cTz + bT y + ~zTDz : Dz + ATy - Is = 0, s > 0}. (2.7) PT(D,AT, I) forsome n-vector p will not change the set of optimal solutions. vice versa. This since z (c ,b ,a J + \ z T D z - ( e T y , 0 T ) ( z ^ y -±zTDz = -pT(D,AT,-I) y -pT0 = o which is a constant. • L E M M A 2.1.2 If (zT,yT,sT) solve (2.7), then (azT, ayT, asT) solves the quadratic problem in which (c T,b T,0 T) in (2.7) is replaced by a[cT,bT,0T) for any scalar a > 0. P R O O F Let (zT,yT,sT) solve (2.7). Then for any scalar a > 0 , acT(az) + abT(ay) + \{az)TD(az) = o?{cTz + bTy + ~zTDz) and, therefore, Q : ( 2 T , y T , 5 R ) solves (2.7) in which (cT, bT,0T) was replaced by a(cT,bT,0T) . • We are now ready to describe the polynomial quadratic programming algorithm whose number of arithmetic steps is independent of the size of the numbers in the vectors c and b . 2.2. T H E Q U A D R A T I C P R O G R A M M I N G A L G O R I T H M The Q U A D R A T I C P R O G R A M M I N G A L G O R I T H M (QPA) described below is a direct generalization to a class of strictly convex quadratic programming problems of Tardos' linear programming algorithm. It uses as input a strictly convex quadratic 27 program (2.1), Tardos' feasibility algorithm, which is a polynomial algorithm to check the feasibility of a system of linear inequalities in time independent of the right hand side and if feasible it generates a basic solution, and a polynomial algorithm for solving convex quadratic programming problems, e.g., Kozlov et al.'s algorithm. The output from QPA is a pair of optimal primal and dual solutions for (2.1). T H E Q U A D R A T I C P R O G R A M M I N G A L G O R I T H M S T E P 1. Use Tardos' feasibility algorithm to check whether {Ax = b, x > 0} is feasible. If not, terminate since (2.1) is infeasible. If feasible, then the positive defmiteness of D guarantees boundedness which in turn implies that the dual constraint set is feasible. Set K = 0. S T E P 2. Let D°x + AoTy - E°s = c° denote the equality system Dx + ATy-Is = c together with S{ = 0 for i € K and let P° = {(xT,yT,sT) : Ax = b, Dx + ATy — Is = c, S{ = 0 for i G K , x > 0, s > 0} . Use Tardos' feasibility algorithm to find a point say (xT ,yT ,sT) in P°. If xTs — 0 terminate with (xT,yT,sT) as an optimal solution to (2.1) and (2.5). S T E P 3. Find the projection {c'T,b'T,a'T) of (cT,bT,0T) onto the null space of (D°, AoT, — E°). Since the rows are linearly independent, this can be done using Gaussian elimination. Recall that for a matrix G with full row rank, the projection onto its null space {x : Gx = 0} is determined by the idempotent matrix P = I — GT(GGT)~1G. I.e., for some vector i its projection onto the null space of G is Px (G(Px) = G(x - GT(GGT)~1Gx) = Gx - Gx = 0) where Px = x - GT(GGT)~1Gx = x - GTp and p = (GGT)~~1Gx. Applying this to our case we obtain (c>T,b'T,a'T) = (c T , 6 T , 0 T ) - pT(D°,A°T,-E°). 28 If (c'T, b'T, a'T) = 0 then solve the problem min{±zTDz : D°z + A°Ty -E°s = 0, s > 0} by Kozlov et al.'s algorithm to obtain an optimal solution (xT = zT + cTD-1,yT,sT). S T E P 4. Let _ (3ra + m)(2n + m)A tt= 11(^,6 ,^^ ^)1100 where A is a maximum absolute sub determinant of (D, AT, — I) , i = l , . 6i = M l t = l , . .., m i = l , . ..,n S T E P 5. Use Kozlov et al.'s algorithm to find an optimal solution (vT,zT) to max {-^zTDz : DoTv - Dz = c, = 6, - £ o r w < a} which is a dual of min {cTz + bTy + aTs + -zTDz : D°z + AoTy - E°s = 0, 5 > 0}. (2.8) 2 Let I = {i: (-E°Tv)i < aa'i - (2n + m)A}. Add the set I to K and go to Step 2. The following lemmas which are extensions of corresponding lemmas in [45] will be used in order to verify the validity of QPA. L E M M A 2.2.1 Let (zT,yT,sT) be an optimal solution to the problem (2.8) where {cT ~bT,aT) = \(cT,bT,aT)] of some vector {cT,bT,aT) e R n + m + n, D is 29 an integral nx n positive definite matrix and D°,A° and E° are integral matrices of appropriate dimensions. Let (vT,zT) be an optimal dual solution for (2.8). Then for any optimal solution (zT, yT, sT) to min {cTz + bTy + aTs + -zTDz : D°z + AoTy - E°s = 0, s > 0} (2.9) 2 (i.e., "unrounded" problem), we have, {-E°Tv)i 0 for each u, we obtain from (2.13) that cTu + bTul + aTu2 + zTDu < 0. (2.14) We will prove the validity of (2.10) by contradiction. Suppose that (2.10) is not valid, i.e., there exists an index i for which (—E o Tv)i < a; — (2n + m)A and > 0. Now since by the optimality condition s t = 0 we have that u\ > 0. Moreover, u | > 0 for all j with sy = 0. If for some j , u2 < 0, then sy > 0 which on the other hand implies (see (2.11)) that (-£'°rC;)y > ay or -ay - (EoTv)j > 0. Now observe that the vector (uT,u1T,u2T) satisfies the system D°u + AoTul - E°u2 = 0. Thus there also exists an integral basic solution to D °u + AoTul - E°u2 = 0, u2 > 0 and u | > 0 for j with sy = 0 that satisfies (2.14). Denote this solution (uT, u1T, u2T) which by Cramer's rule satisfies || (uT,u1T,u2T) ||oo< A . Note that for the dual constraints DoTv — Dz = c and A°v = b corresponding to unconstrained primal variables z and y we have c < D°Tv — Dz < c + e and b < A°v < b + e which in turn implies || —c + D°Tv — Dz ||oo< 1? II — b + A°v ||oo< 15 a n d for j with u | < 0 also | (—aT — vTE°)j \< 1. Combining the above facts and the conditions on (uT, u1T, u2T) we obtain from (2.14) that 0 < -cTu - bTu1 - aTu2 - zTDu = (-cT - zTD)u - 6 T u J - aTu2) = {-cT - zTD)u - Pu1 -aTu2 + vT{D°u + A^u1 - E°u2) = ( - c T - zTD + vTD°)u + {-bT + vTAoT)ul + {-aT - vTE°)u2 -->n-Now since all the components of d are bounded by 1 from above and, with the exception of components dn+m+i,d2n.+m all others are bounded from below by — 1, the validity of the Lemma will follow if we can show that || d Hoo^ (2n + m)A . To that end note that d can be written as d — a -c' ^ -b' -a' i o ; ( D° A0T -E° 0 ) T ( v ) It is easy to see that the projection of d onto the null space of H ~ [-D 0 -E° 0 \D 0 is given by the vector d = a ( - c / T , -b'T, -a'T, - \\ Dz ||oo c'T) (i.e. H • d = 0). Note also that || (c'T,b'T,a'T) \\oo<\\ {c'T,b'T,a'T,\\ Dz ||oo c'T) ||oo and that for any n-vector x , || x ||oo> ^ || i H2 and || x ||2>|| x Hoc are valid. Next, since —a(c'T,b'T,a'T, \\ Dz ||oo c'T) is a projection of d, dh>a\\ (c'T,b'T\a'T',|| Dz c'T) ||2 . Therefore, d ||„o> d II,> (3n + m) > (2n + m)A (3n + m) (2n + m)A (3n + m)) || (c '^ ,^,a^) 1^ ( C^,fe^,a^,|| l l o o Q || 2 34 || (c'r,6'V3\|| Dz ||ooc")||2 > (2n + m)A. C A S E 2. It z is a zero vector, then define the (2n + m) vector d as follows i = 1 , n dt = -ac'i + vTD°, in+j = -ab'j. + vTA°T, l , . . . ,m + m+k - <, _ a f l , f c _ _ T £ o or, written in matrix form if sk = 0 is a row of D°z + Aoi y - E°s = 0, otherwise for k = 1,.... n . ^ - c ' ^ d = a -b' + [ L>° AoT -E° ] v. The proof now follows along the same lines as for C A S E 1. —a • L E M M A 2.2.3 Every optimal pair of primal and dual variables (xT,yT,sT) satisfies = 0 for i G /. P R O O F By Lemma 2.1.1 replacing (c T , 6 T , 0 T ) by (c'T, b'T, a'T) will not change the set of optimal solutions (zT,yT,sT) for min{cTz + bTy + ^zTDz : D°z + A°Ty - E°s = 0°, s > 0}. By Lemma 2.1.2 multiplying the linear part of the objective function by a positive scalar a, one obtains that the set of optimal solutions (azT, ctyT, asT) and the set of variables that are equal to zero in all optimal solutions is unchanged. Finally, Lemma 2.2.1 holds with c replaced by etc' , b replaced by ab' , and a replaced by aa' . • L E M M A 2.2.4 After at most n iterations of QPA one gets a pair of optimal primal and dual solutions for problem (2.1). 35 P R O O F By Lemma 2.2.3 adding constraints S; = 0,i G / where / was determined in Step 5 of QPA does not affect the set of optimal solutions (xT, yT,sT). Recall that the set of optimal solutions form a face of a polyhedron {Ax = 6, Dx + ATy — Is = c, x > 0, s > 0} for which xTs = 0. By Lemma 2.2.2 no more than n iterations of QPA are possible. (In the worst case, where all X{ > 0 and Tardos' feasibility algorithm in every iteration "missed" the desired face, i.e., face with xTs — 0, and where / is a singleton in every iteration, one will have exactly n iterations.) • Let us now calculate the complexity of QPA. Denote by T(A) the complexity of Tardos' feasiblity algorithm when applied to the system {Ax = 6, x > 0}, and by T(A,D) its complexity when applied to the system {Ax = b, Dx + ATy — Is = c, x > 0, s > 0, Si: = 0 for i G K}. Denote by K(A,D) the complexity of Kozlov et al.'s algorithm. Note that we will apply Kozlov et al.'s algorithm only to quadratic programs for which the linear part of the objective function is integral and polynomially bounded by the matrices A and D and with right hand side vector which is zero. T H E O R E M 2.2.5 The Q U A D R A T I C P R O G R A M M I N G A L G O R I T H M has running time polynomial in the size of the matrices A and D and independent of the sizes of the vectors c and b . It runs in 0(n(2n + m) 3 + n(2n + m)log(2n + m)(3n + m) A + T(A) + nT(A, D) + nK(A, D)). P R O O F Step 1 takes T(A) time (i.e., time of Tardos' feasibility algorithm when 36 applied to a linear system with constraint matrix A). Consequently Step 2 takes T(A,D) plus at most 2n comparisons to verify xTs = 0. The Gaussian elimination in Step 3 takes 0((2n + rn)3) time (see Edmonds [16]). Step 4 takes 0((2n + m)log(2n + ra)(3n + m)A) comparisons to find (cT,bT,aT) since || (c T ,b T,a T) = (2n -f- m)(3n + m)A and one can use binary search to obtain (cT,bT,aT). Step 5 takes K(A,D) time. Finally, we need at most n iterations and therefore the claimed complexity follows. • Recall (see Section 1.5) that an algorithm is termed strongly polynomial if all its operations consist of additions, comparisons, multiplications and divisions and if the number of such steps is polynomially bounded in the dimension of the input, where the dimension of the input is the number of data in the input. Further, when the algorithm is applied to rational input, then the size of the numbers occurring during the algorithm is polynomially bounded in the dimension of the input and the size of the input numbers. Thus the polynomial algorithm described in this paper becomes strongly polyno-mial if the size of the entries in A and D are polynomially bounded in the dimension of the data. This clearly provides a strongly polynomial algorithm for, e.g., problems where one minimizes the norm over flow (transportation) type constraints [1,12]. 2.3. E X T E N S I O N T O T H E I N E Q U A L I T Y C O N S T R A I N E D C A S E The Algorithm can be generalized in a straightforward way to work on strictly convex quadratic programs with inequality constraints. To show this consider the 37 quadratic programming problem of the form. max {cTx - -xTDx : Ax < b, x > 0} (2.15) with the same assumptions on the dimension of the problem and the input data as for problem (2.1). Observe that the only difference between problems (2.15) and (2.1) is the existence of inequality instead of equality constraints. The dual of (2.15) has the form min {bTy + -xTDx : ATy - Dx > c, y > 0} . (2.16) 2 Using the fact that D is positive definite and applying the transformation z = x — D~lc we obtain the equivalent pair of primal and dual problems -cTD~1c + max{--zTDz : Az + Iw = b - AD~1c, 2 1 2 Iz - Ix = - L > - 1 c , £ > 0,w > 0} (2.17) and -cTD~1c + min{cTz + bTy + -zTDz : ATy + Dz - Is = 0, y > 0, s > 0} . (2.18) 2 2 Note that an optimal solution (zT,yT, sT) for (2.18) provides an optimal solution x = z + D~1c for (2.15) and (xT,yT) an optimal solution for (2.16). The Karush-Kuhn-Tucker optimality conditions for the pair of primal-dual prob-lems (2.15), (2.16) will have the form 38 y (A 0 I 0 ) y (b) {D AT 0 -I ) w ~ (c ) w s (2.19) x,y,w,s>0 and x s — 0, y w = 0. In this case, at each iteration the algorithm will detect at least one new variable y or s which has to equal zero in all optimal pairs of primal and dual solutions for (2.15) and (2.16). After each iteration one will perform Tardos' feasibility algorithm to detect a basic point from the linear system (2.19). The conditions xTs = 0 and yTw = 0 will then be checked. If satisfied the algorithm will terminate since an optimal pair of primal and dual solutions were found. Otherwise, another integration of QPA will be performed. Instead of problem (2.7) (in the equality case), we will work with the problem of the form min {cTz + bTy + -zTDz : Dz + ATy - Is = 0, y > 0, s > 0}. (2.20) The algorithm will have the same form, except that now STEP 1 will read "Set Ki = and K2 = 0, s > 0} will be equivalent to the system {Dx + ATy — Is = c, y > 0, s > 0, s t = 0 for i £ Ki and yj = 0 for j £ K2}, and the polyhedron P° will be "• = {(*'.»T. «'.'')= [AD °AT I _°7) ( x \ y w s x > 0, y > 0, w > 0, 5 > 0, Si = 0 for i £ Ki, yj = 0 for j £ if 2}-S T E P 5 will now read: "Use Kozlov et al.'s algorithm to find an optimal solution 39 {vT,zT) to max {--zTDz : DoTv - Dz = c,A°v < b,-EoTv < a} 2 which is the dual of min {cTz + bTy + aTs + -zTDz : D°z + A0Ty - E°s = 0, y > 0, s > 0} . 2 Let I = {i: {-E0Tv)l < aa'i - {2n + m)A}, and J = {j: [A°v)l < ab'j - (2n + m)A}. Add the set I to K\ , the set J to K2 and go to Step 2." The complexity of the algorithm will be affected in the sense that at most (n+ra) iterations are now possible. 40 Chapter III P R O X I M I T Y A N D S E N S I T I V I T Y R E S U L T S F O R Q U A D R A T I C I N T E G E R P R O G R A M M I N G In a recent paper, Cook et al. [8], obtained many proximity results for integer linear programming problems with a fixed constraint matrix and varying objective function and right-hand side vectors. In this Chapter we will extend their main proximity results to quadratic integer programming problems of the form max cTx + xTDx s.t. Ax < b (3.1) x integer where c and x are n-vectors, b is an ra-vector, A is an integral m x n matrix and D is a negative semidefinite n x n matrix. In the sequel we will assume that the set {x : Ax < b,x integer} is non-empty, that max{c rx + xTDx : Ax < b} exists, and that D is rational. As stated before, for any matrix A let A (A) denote the maximum of the absolute values of the determinants of the square submatrices of A . For simplicity of exposition we first consider problem (3.1) with diagonal matrix D , i.e. the separable case. We start by showing in Theorem 3.1.1 that, in this case, for any optimal solution z for (3.1) there exists an optimal solution x* for its 41 continuous relaxation such that || z — x* ||oo < nA(A) . Also, if z is a feasible solution to (3.1) which is not optimal then we show in Theorem 3.1.3 that there exists another feasible solution z to (3.1) having greater objective function value and with II z ~ z \\oo < nA(A) . With some additional assumptions we show in Theorem 3.2.2 that if z and z! are optimal solutions for (3.1) with right hand side vectors 6 and b' respectively then || z — z' ||oo < a \\ b — b' \\i +/3 where a and /3 are parameters which depend only on A, D and n . Finally, we show how to extend the above results to mixed-integer quadratic programs. 3.1. P R O X I M I T Y R E S U L T S F O R S E P A R A B L E Q U A D R A T I C I N T E G E R P R O G R A M M I N G Theorem 3.1.1, to follow, provides an upper bound on the distance between a pair of optimal solutions for problem (3.1) and its continuous relaxation respectively. The bound obtained depends only on the number of variables n and the largest absolute subdeterminant of the matrix A , and is independent of the vectors 6, c and the matrix D . T H E O R E M 3.1.1 Let A be an integral m x n matrix, D a diagonal negative semidefinite n x n matrix, b an m-vector and c an n-vector such that the set {x : Ax < b , x integer} is non-empty and max{c r i + xTDx : Ax < b} exists. Then (i) For each optimal solution x for max{c T i + xTDx : Ax < 6} there exists an optimal solution z* for max{cTx + xTDx : Ax < 6,x integer} with | | 2 * — x | | o o < nA(A) , and (ii) For each optimal solution z for max{cTx + xTDx : Ax < b, x integer} there 42 exists an optimal solution x* for m&x{cTx-\-xTDx : Ax < 6} with ||x* — z \\oo < nA(A) . PROOF Let z and x be optimal solutions for problems (3.1) and its continuous relaxation respectively. Further assume without loss of generality that the first £ components of x — z are nonpositive and the last n — I components are positive. Partition the rows of the matrix A into submatrices A\ and A2 such that A\x > A\z and A2x < A2z . Consider the cone C = {x = (x\,..., xg, X£+i,..., xn) : A\x > 0, A2x < 0,xl = (xl5...,X£) < 0,x2 = ( x £ + 1 , . . . ,xn) > 0} . Clearly, C is nonempty since x — z G C . Let G be a finite set of integral vectors which generate C . By Cramer's rule for each g G G , || a ||oo < A(yl') = A(A) , where A' = ^ - / f o j and 7x (resp., I2 ) is an I x £ (resp., (n — £) x (n — £) ) identity matrix. Now, since x — z G C there exists a finite set {g\,... ,gt} C G and scalars > 0, i = 1,.. . , t such that x — z = Yll=i ai9i • In order to prove part (i) we will define the vector z* = z + 5Zi=i L^t'Jfft • ^ e will show first that z* is feasible and optimal for (3.1). Coupling with the fact that z* = x — Yli=i{ai ~~ [ai\)di we will obtain that || z* - x \\oo< tA(A) < nA(A) . The proof of the feasibility of z* is identical to the proof given in Cook et al., but we will state it for the sake of completeness. Observe that z* is integral and satisfies Aiz* = Ai{x- Ei=i(«i - lai\9i)) < , A2z* = A2{z + Y?i=Aail9i) < A2* and therefore Az* < b . To prove optimality, recall (see e.g. Dorn [13]) that the dual of the continuous relaxation of (3.1) is given by min{u r6 - xTDx : uTA - 2xTD = cT,u> 0} (3.2) with complementary slackness conditions uT(b— Ax) = 0 . Now, if we partition bT to (bi,b2) corresponding to the partition of A to (Ai,A2) we have A2x < A2z < 43 &2 , and thus it follows from the above condition that for a vector uT = (uf, uT) of optimal dual variables for (3.2) we have u2 = 0 and u7Ai — cT + 2xTD . Also, since for every y £ C, vliy > 0 we have c T y + 2x r7Jy > 0. (3.3) Finally, using (3.3) and the fact that Ei=i [ai\9i £ C* we obtain t t cTz* + z*TDz* = cT(z + Yia^9i) + zTDz + 2z T£(]T[a lJff l) i=i i=i t t = cTz + zTDz + cT{Y[oct\gi) +2xTD{J2l cTz + zTDz+ QTQa.-J -2a t-)ffi) r^£l- a«J^)-i=l i=l It remains to show that (Ei=i(l_aiJ — 2a:t)<7t)T.D(E*=1 [ai\9i] > 0 . To that end recall that gi £ C and hence can be partitioned into ( 0 . Diagonality and negative semidefiniteness of D implies gjDgj < 0 for all i,j = Coupled with the fact that [atJ — OLX < 0 and [a,-J > 0 for all i we obtain ( E U ( M - 2a0ft)TD(E!=iL«iJft) = ES=i (W - 2a t ) [« I j 5 f /J f f l + E-=i Ej>,-(([a.-J - 2a,-) Lay J + ([ayJ - 2ay)[a t-J)^I>W > 0 . 44 Thus, cTz* + z*7Dz* > cTz + zTDz and part (i) is valid. Also, the optimality of z and z* implies (3.4) In order to prove (ii) define first the vector x = (3.5) which will be shown to be feasible and optimal for the continuous relaxation of (3.1). The feasibility of x* can be shown in a similar manner to the way it was done for z*. To show optimality observe that cTx* + x*TDx* = cTx - cT(Y2leti\gi) + xTDx - 2xT D(^2[ai\gi) t ' = i i = i t t + (£L«.:I*)t0(£L«.-J*) i = l i = l t t t = cTx + xTDx - ( c r ( E | a i J g t ) + (21 + ^[ctilgif D(Y,[<*i\gi)) i=i i=i i = i t t + 2 ( ^ ( [ a t J - a i ) » i ) T - D ( X ] L a » - J ^ ) = ° T * + i = i i=i t t + 2(^2{[ai\-ai)gi)TD(^2[ai\gi) > cTx + xTDx, where the equalities follow from (3.5) and (3.4) and the last inequality follows from the fact that i = i (£(La,-J -ai)gi)TD(J2l"i\9i) > 0. Finally, from (3.5) || x* - z ||oo< tA(A) < nA(A). • 45 Using the example provided by L. Lovasz for the linear case (see Schrijver [43], page 241) it can be verified that the bound is tight. Observe that if lower and upper bounds on the variables X{ are known one can naturally improve the bound on the difference between an integral and continuous optimal solution. Indeed, if P = {x : Ax < b, ti < Xi < Ui} then denote by B = min { nA(A), max {| Uj — ti |: i — l,...,n} } . Clearly B is a valid bound. Trivially, in the case of a 0 — 1 quadratic programming problem we have B = 1 . As a consequence of Theorem 3.1.1 we obtain Corollary 3.1.2 which provides a bound on the difference between the optimal objective function value of (3.1) and its continuous relaxation. C O R O L L A R Y 3.1.2 Let A be an integral m x n matrix and b an m-vector. Let P = {x : Ax < b} be a nonempty bounded polyhedron having an integral point. Then, , max{cTx + xTDx : Ax < 6} — max{cTx + xTDx : Ax < b, x integer} < n A ( 4 ) ( | | c ||i +2nA(A | 6) || Z? ||oo). P R O O F Denote by z (resp., QI ) an optimal solution (resp., the optimal ob-jective function value) for max{cTx + xTDx : Ax < b, x integer} and by x (resp., QC ) an optimal solution (resp., the optimal objective function value) for its contin-uous relaxation with || x — z | | o o < nA(A) which clearly exists by Theorem 3.1.1. Then, QC -QI = cTx + xTDx - cTz - zTDz = cT(x - z) + (x + z)TD(x - z) <\\ C ]ji|| X - Z || oo + || X + Z || i | | D(x - 2) || oo <|| C || i | | X - Z || oo + || X + Z || i | | D ||oo || Z - Z ||oo 46 where || D = max{|d\i| : i = . Now since for a vector x , || x \\i< n || x ||oo and for a pair of vectors x, y \\ x + y ||oo<|| % ||oo + || V ||oo , we have II x + z ||i< n(\\ x Hoc + || z ||oo) < n(A(A \ fe) + A(A | fe)) by the boundedness of the polyhedrons. Therefore, QC-QI<\\x-z lu (|| c Id +2nA(A | fe) || D IU) < nA(A)(\\ c Id +2nA(A | fe) || D W^) . • Theorem 3.1.3, which follows, will show that for any integral solution z of Ax < fe which is not optimal for (3.1), there exists an integral solution z of Ax < b which has greater objective function value and with || z — z ||oo< nA(A) . T H E O R E M 3.1.3 For each integral solution z of Ax < b , either z is an optimal solution for (3.1) with a diagonal matrix D or there exists an integral solution z of Ax < b with || z — z < nA(A) and cTz + zTDz > cTz + zTDz . P R O O F Let z be an integral solution of Ax < b which is not optimal for (3.1) with a diagonal matrix D . Then there exists an integral solution z* of Ax < fe with cTz* + z*TDz* > cTz'+ zTDz . As in the proof of Theorem 3.1.1, without loss of generality, assume that the first I components of z* — z are nonpositive and that the last n — £ components are positive. Partition the rows of A into submatrices Ax and A2 such that A\z* > A\z ,A2z* < A2z and consider the cone C = {x = (xi,... ,X£,xe+1,... ,xn) : Aix > 0 , A2x < 0 , xl = ( x l s . . . ,xt) < 0 , x2 = (x^+i,... ,xn) > 0} . Since z* — z G C , we can write z* — z = Y A = I ai9i where a t > 0 for i'. = 1,..., t and {gt,...,gt} is a subset of integral generators of C . By Cramer's rule || gi ||oo < A(A) for every i . If oti > 1 for some i , then z + 0{ = z* — 5Zj = i J ? i , ajSj — iai ~ l)°i is a n integer feasible solution of Ax < b . 47 Now, if in addition cTgi + (2z + gi)T Dgi > 0 , then cT(z + gi) + (z + gi)TD(z + gi) > c Tz + zTDz and by choosing I = z + we are done. Thus, let us assume that for all i for which ai > 1 we have cTgi + (2z + g i ) T D g i < 0. (3.6) Next, define the integral vector t t z = z* - 5^|a«j0i = z + X^" 1' _ La*J)0* i=i i'=i which is a feasible solution of Ax < b and satisfies || z — z ||oo < tA(A) < nA(A) . To complete the proof we will show next that the objective function value at z is larger than at z* . Indeed, cTz + zTDz = cTz* + z*TDz* - c r(E*=i L«iJ*) - 2^ T7J(E! = 1 + ( E L i L^JffO^lELi l«i\9i) = cTz* + z*TDz* - ( ' r ( E ! = i L«.-J*) + (2* + 2 £,-=i ^ - E!=i L«.-J*)tJ>(E!=I L«.-J w)) = cTz*+z*TDz* - (cT(El1M9i)+2z'rD(EL1Mgl)) -(E!=i(2ai - l^J)(E!=i |ajtt) = c rz* + z*TDz* -(c r(El=iL«dffO + 22 Ti?(E:=iKJ^)) - E!=i(2a,- - L"d) W ^ . -- E L Ej>i((2a,- - L«iJ)L«>J + (2«y - L « , J ) K J ) ^ ' > ' T * * + z*rI>** -(cr(E!=iLatJj/i) + 2«r£>(E!=iL«.-JffO) - EL i ( 2 ^ - [a^la^gfDgi where the inequality follows from the fact that gi (E C , [a;J > 0 and 2^ — [a,-J > 0 for all i which in turn implies ((2a, — [at'J)LayJ + (2ay — [ayj) [ai\)gf Dgj < 0 . 48 Therefore t cTz + zTDz > cTz* + z*TDz* - cT{Y^[<*i\gi) i=l t - ^[a»J(2a-|- (2at- - [al\)gl)TDgl • i=l Now, since for cc,; < 1 , |_a;j = 0 , the only nonvanishing terms in the above summation correspond to a; > 1 . Further, since D is negative semidefinite and 2a; - [cti\ = a t + («i - [a,-J) > a, > 1 , (2at- - of Dg^ < gf Dgt- . This implies, using (3.6), that t cTz + zTDz > cTz* + z*TDz* - Y, L«d {cT9x + (2* + O i ) T ^ f ) t = i > C T 0 * + 2 * T D 2 * > CTZ + 2 T £ > 2 . • Observe that the above Theorem assures finite "test set" for problem (3.1). In other words, one has to check only a finite number (which depends only on the constraint matrix and the number of variables) of integral vectors in order to obtain a better solution or to verify the optimality of the current solution. Although the set is finite, it might require the exponential (in n ) number of comparisons in order to get better solution. 3.2. S E N S I T I V I T Y F O R Q U A D R A T I C I N T E G E R P R O G R A M M I N G : T H E R I G H T H A N D S I D E C A S E In this Section we show that if z and z' are optimal solutions for (3.1) with negative definite matrix D , constraint matrix A of full row rank and right hand side vectors 6 and b' respectively, then || z — z' ||oo £ a || b — b' \\\ +/? where 49 a and 0 are parameters which depend only on A, D and n . Observe that restricting D to be negative definite implies that the continuous relaxation of the quadratic integer program has a unique optimal solution whenever the polyhedron P — {x : Ax < ft} is nonempty. For simplicity of exposition Theorem 3.2.2, to follow, will be stated and proved for the separable case. The general case will be discussed in Section 3.3. Theorem 3.2.2 is using a special case of Theorem 2.1 in [11] in which changes in the linear cost coefficients are considered. For completeness we will state the part of Theorem 2.1 in [ll] we use. L E M M A 3.2.1 Let A be an m x n matrix, b an m-vector, c and c' n-vectors and D = diag(di) an n x n diagonal negative definite matrix. Let x (resp., x' ) be the optimal solution for max{c rx + xTDx : Ax < ft} (resp., max{c' rx + xTDx : Ax < b}) . Then \\ X - x' Hoo^ — || C - c' ||i where K = — max{d, : i = 1,. . . , n} > 0. P R O O F See Daniel [11], Theorem 2.1 . • T H E O R E M 3.2.2 Let {Ax < ft} = {AEx = bE,Ajx < ft/} and {Ax < ft'} = {AEX — b'E, Ajx < b'j} have integral solutions where A is an integral mxn matrix of full row rank and ft, ft' are m-vectors. Let c be an n-vector and D = diag(d{) a diagonal, negative definite n x n matrix. Then (i) For every optimal solution z for max{cTx + xTDx : Ax < ft, x integral} and 50 every optimal solution z' , for max{c rx + xTDx : Ax < b',x integral} we have II z - z' ||oo < nA{A)(H{D) \\ b - b' ||i +2) where H(D) = max{|a\| : i = 1,.. . , n}/ min{|o!i| : i = l,...,n}. (ii) Assume, in addition, that D is integral. If f(b) (resp., f(b') ) is the optimal value for max{c rx + xTDx : Ax < b , x integral} (resp., max{cTx -f xTDx : Ax < b' , x integral} ), then |/(6) - /(6')| < nA{A){\\ c ||x +n \\ D {2nA{A) + A{A\bc) + A{A\bc'))){H{D) | |6 -6' | | i+2) where A = PROOF We will start by determining an upper bound on the difference between the optimal solutions to the corresponding continuous quadratic programs. To that end let us first write the dual of the continuous relaxation of (3.1) in primal form - m a x { ( 0 ^ ) Q + ( * V ) ( £ o)(o)> (3.7) -2D AT ) ' c ' S.t. 2D -AT —c . o -I , \u ) . 0 , Now, we will treat changes in the right hand side of (3.1) by considering changes in part of the linear cost coefficient of (3.7). To that end, recall that if x is an optimal solution for max{c rx -+- xTDx : Ax < b} then x is also an optimal solution for the linear programming problem max{(c r + 2xTD)x : Ax < 6} see e.g. Dorn [13]. Thus if (x,u) and (x',u') are, respectively, optimal of (3.7) and the quadratic programming problem obtained from (3.7) by replacing b by b' then, for every feasible solution (x, u) of (3.7) we have - uTb + 2xTDx < -uTb + 2xTDx - uTb' + 2x'TDx < -u'Tb' + 2x'TDx'. (3.8) 51 ( A 0 I - 2 D AT I Substituting (x',u') (resp., (x, u) ) in the first (resp., second) inequality of (3.8) adding and rearranging terms results with 2{x' - x){-D){x' -x)< ( « ' - u)T{b - b') <|| s' - fi lull b-b' HJ. (3.9) Now, since both (x,u) and (x', tt') are feasible for (3.7), (x' — x, u' — u) = (v,w) satisfies 2Dv — ATw = 0 and hence AT(u' -u) =2D{x'-x). Using the fact that AT is integral and has full column rank, Cramer's rule, and a property of determinants we obtain || « ' - fi | | o o < A(AT | 2D(x' - x)) < 2nA(A) \\ D H^ H i ' - z ||oo • (3.10) Now let K = min {—d{ : i = l,...,n} = —max {di : i = l,...,n} . Substituting (3.10) into (3.9) and observing that K > 0 results with 2K || x - x ||^ < 2K || x' - x \\\< 2nA(A) \\ D | |oo| | x' - x W^W b' - b ||i (3.11) or II *' - x lU^ n A ^ D H°° || b> _ 6 ||1= nA(i)ff(i?) || 6' - 6 ||i . Using the fact that the continuous relaxation of our quadratic integer program has a unique solution and applying part (ii) of Theorem 3.1 it follows that for every pair of optimal solutions z and z! for the quadratic integer programming problem with right hand side vectors b and b' respectively || 2 2 | | o o ^ | | 2 X \\oo || X X \\oo + || X Z | | o o ^ (2nA(A) + A ( i |*|) + A ( i |f |))(tf(£>) || 6 - fc' ^ +2). • 3.3. E X T E N S I O N T O N O N S E P A R A B L E Q U A D R A T I C M I X E D - I N T E G E R P R O G R A M M I N G The results given so far can be easily extended to separable quadratic mixed-integer problems of the form 53 T * T 1 f max c1 x + c 2 y + x D\x + y D2y s.t. Ax + By b (4.1) 0 < x < 1 x integer where x is an n x l vector, b is an m x 1 vector, g(x)T = (gi(z),...,om(x)) , / and gi , i = l , . . . , m are real-valued, differentiable functions and, furthermore, / is convex and gi , i = 1, ...ra are concave on Rn . We will also assume that the continuous relaxation of (4.1) satisfies the Kuhn-Tucker constraint qualification (see e.g. [49]), which is automatically satisfied if the constraints are all linear. Explicitly assume the nonlinear integer programming problem (4.1) was solved by implicit enumeration and some small changes have been made in the right hand side or objective function coefficients of (4.1). The question we would like to answer is what information from the implicit enumeration tree, if at hand, will provide us with bounds on the optimal value of the perturbed problem. Before attempting to answer the above question observe that the continuous relaxation of (4.1) is a convex nonlinear optimization problem for which we can state the dual (see Wolfe [49]) max b u + e v +f(x) — u g(x) — (v +v ) x 58 s.t. uTVg(x) + v1TI + v2TI = Vf(x) (4.2) u , v1 > 0 , v2 < 0 where eT = (1,..., l) and u (resp., v1 , v2 ) denotes the vector of dual variables associated with the constraints g(x) > b (resp., x > 0 , — x > —1 ). Solving the primal (i.e. the continuous relaxation of (4.1)) with some readily available computer code, one can obtain the value of these dual variables as a byproduct (e.g. using MINOS -Modular In-core Nonlinear Optimization System). R E M A R K 4.1.1 If all the constraints of (4.1) are linear, then the objective function of (4.2) can be written as (see e.g. Dorn [14]) bTu + eTv2 + f(x) — xTVf(x) . Let us start by solving (4.1) using implicit enumeration. In doing so we con-struct a tree with node 1 corresponding to the continuous relaxation of the original problem. At each node t of the tree one solves R\b) = min f{x) s.t. g(x) > b L} 0 otherwise v y < 0 . Let (uJ, xf, vf) be the associated dual solution obtained by solving (4.3). Using (4.2) the optimal objective function value of (4.3) equals R\b) = bTut + Y, {vt)j + Y min<0' > + f ^ - ufg{xt) - Y (vt)i - Y m i n { ° 5 (vt)j) jeF* jcF\(FluF<) = bTut + f(xt) - ufg(xt) . R E M A R K 4.1.2 If the constraints in (4.1) are all linear the optimal objective function value of (4.3) can be alternatively written as R'ib) = bTut + Y + Y m m < 0 ' + ttXt) - x?Vf{xt) . jeFl jeF\{F*UF*) Assume now that (4.1) is perturbed by replacing b by a new vector d . We would like to use the information from the implicit enumeration tree to derive a lower bound on the value of the perturbed problem. To this end notice that the dual variables ut, Xf, i>t derived at node t of the tree remain feasible, but not necessarily optimal, for the perturbed problem. Therefore, a lower bound on the objective function of the perturbed problem is given by Rt{d) = dTut + f{xt)-ujg{xt) . (4.4) Let I (d) denote a lower bounding function on the objective function value 60 /*(d) . For terminal nodes with feasible solutions for .#'(6) we can set f (d) = R*(d) . For terminal nodes with no feasible solutions set where 14 and Xt used in this case in the evaluation of R (d) are part of the dual variables (uJ, xf, vj) associated with the minimization of the sum of the infea-sibilities. (Having a problem with only linear constraints and using Remark 4.1.2, both Ut and vt and xt = 0 can be used in the evaluation of R*(d) .) For each nonterminal node t define where L(i) and R(i) are the two offsprings of t. Theorem 4.1.3 below provides us with a lower bound on the objective function value of the perturbed problem obtained from (4.1). T H E O R E M 4.1.3 lld) is a lower bound for N(d) . P R O O F We show first that for any node t and any d , I*(d) < /*(d) is valid. If t is a terminal node then the inequality follows from the definition of l'(d) and the convention that an infeasible minimization problem has objective function value +oo . Now, if t is not a terminal node, / 4(d) > Rt(d) > R*(d) . Further, from the implicit enumeration J*(d) = min {ILM{d) , IRW{d) } > min (I L ( t ) (d) , fW{d) } . Therefore, /*(d) > max (R*(d) , min {IL ( f )(^) , lm{d)} and by (4.5), J*(d) > I^d) -oo if R'(d) < 0 + oo if R*(d) > 0 I*(d) = max {R*(d) ,min (I L ( t ) (d) , lm{d)}} (4.5) 61 ^(d) . By induction towards node 1 of the tree we get N(d) = > I 1^) • • Theorem 4.1.4 to follow provides an answer to the following question. Assume we augment problem (4.1) by adding a new 0 — 1 variable, say x n +i , resulting in the addition of new linear terms in the constraints and some terms (not necessarily linear) in the objective function. Under what conditions will x n +i remain at zero level in the optimal solution to the modified integer nonlinear program? T H E O R E M 4.1.4. Suppose after solving (4.1) to optimality the problem was enlarged by introducing a new 0—1 variable, say x n + 1 , resulting in the addition of a new linear term, say a(Xn+i , to each constraint i = 1, .. . ,m , and a number of new terms in the objective function given by / ( x ) x „ + i . Then there exist an optimal solution to the new problem with x n + i = 0 if IF, >N{b)-l\b-a) where a T = (aj,. . . ,am) , fpl = /(x) and x is the optimal solution to the initial problem (4.1). P R O O F Suppose x n + i = 1 is in an optimal solution at node 1 . The remaining optimal values can thus be found by solving problem (4.1) with right hand side b—a . By Theorem 4.1.3, N(b — a) > I 1 (6 — a) and therefore a solution to the enlarged problem with xn+\ = 1 has objective function not less than N = fFl + I 1 ( 6 - a ) . Now, if N > N(b) the solution to the initial problem remains optimal. • 62 4.2. T H E M I X E D - I N T E G E R C A S E Consider now problem (4.1) where only a subset of variables is restricted to be integer. Theorem 4.1.3 can be carried over without any changes. For Theorem 4.1.4 to be valid the assumptions remain unchanged, while the result should read : "Then there exist an optimal solution to the new problem with xn+i — 0 if f{x) > N{b) - I 1 ( 6 - a) where f(x) denotes the sum of the new terms in the objective function evaluated at x , the optimal mixed-integer solution to the original problem, and x n +i = 1 •" Next, consider the nonlinear mixed-integer program of the form (4.1) except that the integer variables are not necessarily restricted to be zero or one. The bound (4.4) can be derived in a straightforward manner as for the 0 — 1 case and Theorem 4.1.3 will be valid without any changes. However, if xn+\ is an integer variable restricted to the interval [0,17], then we will restrict ourselves in Theorem 4.1.4 to the case in which the added terms in the objective function are linear in . In this event the bound can be improved as follows. If min{ / (x)x n + 1 + I1 (6 - axn+l) : x n + i = 1,...,U} > N(b) then the solution to the initial problem remains optimal. 4.3 T H E Q U A D R A T I C 0-1 P R O B L E M : A N E X A M P L E Consider the quadratic integer programming problem Q[b) — min c x-\—x Dx 63 s.t. Ax > b (4.6) 0 b,0 < x < 1} is nonempty. Indeed, as stated in Section 1.4, if D is not of the desired form (4.6) can be replaced by an equivalent problem in which the objective function of (4.6) is replaced by [ J - \ \ e T ) x + \xT{\(D + DT) + \I)x where A is a positive scalar such that \(D + DT) + AJ is positive semidefinite, see for example [25]. Taking into account Remark 4.1.1, the dual of the continuous relaxation of (4.6) can be written as (see e.g. Dorn [13]) max o u + e v x Dx 2 s.t. ATu - Dx + Iv < c (4.7) u>0, -v > 0 . This since the continuous relaxation of problem (4.6) can be written as min {-xTDx + cTx : Ax >b,x>0} 64 where A = ^ A ^ and b = with eT = (!,...,!) . Following Dorn [13] Type I, page 60, its dual is 1 max {— w Dw + b z : A z — Dw < c, z > 0 } 2 which is equivalent to l /T /T /T max it; Dw + b z\ — e zo 2 rp S.t. A 2 ! - Iz2 — Dw < c Z\ > 0 , 2 2 > 0 . Now, replacing w by x , z\ by u and 22' by —v results in (4.7). From Remark 4.1.2, the bound (4.4) in this case becomes Rt{d) = dTut+J2{vt)j+ m i n ( 0 , {vt)j} ~ \xfDxt . (4.8) It is easy to see that R t(d) can be further improved using information obtained from other nodes in the tree (see also [42] for the linear case). Indeed if the dual prices of all nodes, say N , of the implicit enumeration tree are used, R f(d) can be improved to R*(d) = max {dTus + V (va)j + Y] min{0, (us)y> - \xTsDxs } . seN ^ L jeF* jeF\(F*UF*) For problem (4.6), Theorem 4.1.2 specializes to " Suppose after solving (4.6) to optimality the problem was enlarged by adding a new column, say a n + i to A . Then, there exist an optimal solution to the new problem with x n + i = 0 if cn+l + ^ ( d n + l , n + l + rf*>+1) - Q(b) ~ ^ (b ~ an+l) , 65 where Fi is the set of variables fixed to 1 in the optimal solution for (4.6) with Q(b) its optimal objective function value." E X A M P L E 4.1 form Consider the quadratic integer programming problem of the Q — min 65xi - 10x2 + 7x 3 + 58x4 — 8x 5 + 23x6 — 8x^ 4 + 16xxx6 + 4x 4 x 6 s.t. 70X] - 20x2 + 30x3 - 20x4 + 90x5 + 100x6 > 200 lOOxj + 30x2 - 30x3 + 80x4 - 5x5 + 70x6 > 100 0 < X{ < 1 , X{ integer , i' = 1, ...,6 . The equivalent objective function of the form cTx+\xTDx with positive semidefinite matrix D has (15,-10,7,50,-8,5) and D = /100 0 0 0 0 0 0 0 0 - 8 0 0 16 0 0 0 0 V 16 0 0 -8 0 16 \ 0 0 0 0 0 0 4 0 0 0 4 0 36/ The optimal solution to the continuous relaxation equals (0.20588,1,0.5196,0,1,1) with objective function value of 17.139 . The corresponding 0 — 1 optimal solution equals (0,0,1,1,1,1) with objective function value of 84 . The continuous relax-ation subproblems were solved using QPSOL, a F O R T R A N package for Quadratic Programming developed at Systems Optimization Laboratory, Department of Oper-ations Research, Stanford University, and implemented on A M D A H L 470 V-6 com-puter model II. QPSOL minimizes an arbitrary quadratic function subject to linear constraints where upper and lower bounds on the variables are handled separately. It requires an initial estimate of the solution and a subroutine to define the quadratic part of the objective function. Among output arguments, the Lagrangian multipliers 66 for each constraint are given. In the case of an infeasible constraint set , the minimum of the sum of the infeasibilities was determined using the LINDO package. CR X i = 1 xi = 0 XQ = 1 xT = (1,1,0,0,1,1) , z = 86 xe - 0 infeasible X4 = 1 xT = (0,0,1,1,1,1) , 2 - 8 4 X4 = 0 infeasible Figure 4.1: Branch and Bound Tree for Example 4.1 Figure 4.1 describes the enumeration tree associated with example 4.1. The number above each node corresponds to the node index while the entry in each node represents the branching choice. For terminal nodes with feasible solution an optimal solution and the optimal objective function value are denoted by x and 2 , respectively. The lower bounding functions R*(d) for the continuous subproblems solved at each node are given by R7{d) = dx + d2 - 265 , Re{d) = 0.5di - 16 , 67 R5{d) =di+ 0.25d2 - 206.25 , R4{d) = 86 , R3{d) = 2.756<2X + 1.504d2 - 640.88 , R2(d) = 0.318**! - 1.68 , R :(d) = 0A41di + 0.207d2 - 91.75 . The lower bounding functions for the corresponding integer problems are T7(J\ _ / - ° ° i f R7(d) < 0 " 1 ' ~ \ +oo if R7{d) > 0 ' f{d) = R6{d) , T5(J\ - I -oo if R5[d) < 0 - [ > ~ \+oo if R5{d) >0 ' I4(<2) = R4{d) , f{d) = max {R3{d) , min {I6(b (5.1) 0 < x < 1 x, integer where A is an m x n matrix, D is an n x n symmetric matrix, c and x are n-vectors and b is an m-vector. Problem (5.1) is a natural representation of many problems in, for example, finance [34] and capital budgeting [31]. Different approaches for solving the above problem can be found in the literature, e.g., lin-earization methods where the quadratic problem is transformed into a linear 0—1 or a mixed-integer program can be found, respectively, in Watters [47] and Glover [22]. Algorithms based on a branch and bound method have been proposed by many authors, e.g., Mao and Wallingford [37], Laughhunn [31] and Hansen [26]. McBride and Yormark [38] gave an implicit enumeration algorithm in which at each node they solve a quadratic programming relaxation of a corresponding integer subproblem using Lemke-Howston's complementary pivoting algorithm. It is conceivable that a 70 success of such an implicit enumeration algorithm depends greatly on the efficiency of the quadratic programming algorithm used. Although, at present, the polynomiality of an algorithm can not always be identified with real world computational efficiency or practicality, it is an important theoretical result which leads the research efforts in the direction of constructing efficient problem oriented polynomial algorithms. As stated in Section 1.5. Kozlov, Tarasov and Hacijan [30] provided the first polynomial time algorithm for convex quadratic programming problems. For a class of strictly convex quadratic programming problems, in Chapter II we proposed a polynomially bounded algorithm in which the number of arithmetic steps is independent on the size of the numbers in the linear cost coefficients and in the right hand side vector. We show in this Chapter how to replace the objective function of a quadratic 0 — 1 programming problem with n variables by an objective function with integral coefficients whose size is polynomially bounded by n , without changing the set of optimal solutions. We will use Frank and Tardos' [19] algorithm which in turn uses the simultaneous approximation algorithm from Lenstra at al. [33]. The above result assures that the running time of any algorithm for solving quadratic 0 — 1 programming problems can be made independent of the size of the objective function coefficients. This since the equivalent problem can then be solved by e.g. an implicit enumeration algorithm in which at each node the continuous relaxation of the corresponding integer subproblem is solved in polynomial time independent of the size of the objective function coefficients. Observe that since (5.1) is a 0—1 programming problem then a constraint i with b{ > Ey=i \aij\ 15 clearly infeasible. This since Yl^=iaijxj — Ey=i \aij\ < ^» • Therefore, assuming that (5.1) is feasible automatically assures that the entries of b can be bounded by the entries of the constraint matrix A . Note, however, that this 71 does not imply that the input size of 6 is polynomially bounded by the input size of the size of p and (or) q is not polynomially bounded by the size of the entries of A . In any event, if the entries in the constraint matrix are polynomially bounded by the number of variables and/or constraints, then the continuous relaxations of the integer subproblems can be solved in strongly polynomial time using the algorithm presented in Chapter II. Recall that in a strongly polynomial algorithm the number of elementary arithmetic operations (i.e., additions, comparisons, multiplications and divisions) is independent of the size of the input and is polynomially bounded in the dimension of the input (i.e., the number of data in the input). In Section 5.1. we state the problem and give some preliminary definitions. The extension of Frank and Tardos' preprocesing algorithm to quadratic 0 — 1 problems is given in Section 5.2. For simplicity of exposition we will consider a problem with an objective function in homogeneous quadratic form A since 6 can be a rational vector ^ ( pT = (p,, ...,pn) , o r = (gi, qn) ) where 5.1. S E T U P O F T H E P R O B L E M mm 2 s.t. Ax>b (5.2) 0 < x < 1 x integer . This can be done without loss of generality since the transformation of the objective function of (5.1) to the homogeneous quadratic form given in (5.2) can be easily 72 achieved as stated in Section 1.4. For the sake of completeness we will give here some details. For example, by adding a new variable y = 1 problem (5.1) can be restated as (5.2) with (;)< ' A 0 ' ' b ' A = 0 1 and b = 1 , o -1 . ; -1 . = I , -.,n we h ave cTx + ^ xTDx = \xT{D + 2C)i where C is a diagonal matrix with ca = c t, i = 1, ...,n. In the sequel we will use some vector and matrix norms as defined in Section 1.1. Let S = {x £ f?2 : Ax > 6}, where B2 = {0,1} . A vector v G S is said to be a feasible solution of (5.2) while a vector z G S for which zTDz < vTDv for every v G 5 is said to be an optimal solution for (5.2). The following lemma will justify the algorithm to follow. L E M M A 5.1.1 If for every u,v G S we have sign(u — v)TD(u + v) = sign(u — v)TD(u + v) for some symmetric matrices D and D, then problems (5.2) with matrices D and D , respectively, have the same set of optimal solutions. P R O O F We will show that every optimal solution of (5.2) with matrix D (resp., D ) in the objective function is optimal for problem (5.2) with D (resp., D ) in the objective function. To that end suppose that for some u G S u1Du < 73 vTDv is valid for every o £ S. Then it follows from the symmetricity of D that uTDu - vTDv = (u- v)TD(u + v) < 0 for all v e S. Now s ign( t t - i , ) r D(u + t,) = {- 1 !J TT^TT < ^ 1.0 if u J Du — v Dv . By the assumption, sign(u — v)TD(u + v) = sign(u — v)TD(u + v) . This means that uTDu < vTDv for every v G 5, which in turn implies optimality of u for problem (5.2) with the matrix D in the objective function. • 5.2. S I M U L T A N E O U S A P P R O X I M A T I O N O F O B J E C T I V E F U N C T I O N C O E F F I C I E N T S Frank and Tardos [19] presented an algorithm which replaces a rational cost coefficient vector w of a linear programming problem with an integral vector w, without changing the set of optimal solutions. Their algorithm uses a revised ver-sion of Lenstra, Lenstra and Lovasz's simultaneous approximation algorithm (LLL algorithm) which is strongly polynomial. For the sake of completeness we will state Frank and Tardos' algorithm (F-T algorithm) [19] : I N P U T w = (w(l), ...,w(n)) rational vector and an integer TV with 1 < TV < 2n! O U T P U T w = (td(l),...,w(n)) integral such that || w < 2n'+2n2+2nJVn(n+2) and sign (w,b) = sign (w,b) whenever b is an integral vector with || b \\\< TV. 0. Let M = 2" 2 +"+ 1 TV n + 1 , wi=w, w = 0 and i = 1. 1. Let Ul' = T T — T i — 1 I!*".- Moo 1 2. Apply the revised L L L algorithm to TV and w'^l), ..^w^n). Let p t = 74 ( p i ( l ) , P i { n ) ) and g; denote the output. Then || qiw\ — pi | | o o < 1/-W and 1 < qx < 2 n 2 + n N n . 3. Let Wi+\ = q{Wi — pi and w = Mw + pi. If it>;+i = 0, H A L T . Otherwise let i = t + 1 and G O T O 1. E N D . The algorithm presented above can be generalized into a preprocessing algorithm that will transform the objective function coefficients of a 0—1 quadratic program-ming problem into integer coefficients whose size will be bounded by a polynomial function of n and for which the set of optimal solutions remains unchanged. As stated above, without loss of generality, we will assume the homogeneous form in the objective function. P R E P R O C E S S I N G A L G O R I T H M F O R Q U A D R A T I C 0-1 P R O B L E M S I N P U T D = (dij) an n x n symmetric rational matrix and an integer N = An2. O U T P U T D = (dij) an nxn symmetric integer matrix with 6 < 2^+2h~+2hNh(h+2\ where h = n( n+ l) a n d £ — max{dij i,j = 1, ...,n} and for which sign(u — v) TD(u + v) = sign(u — v) TD(u + v) for every integral u,v with || u | | o o < 1 and || v | | oo< 1-S T E P 1. Construct a rational vector d = ( d u , d l n : d 2 2 , d2n-, •••> dnn) where dij is the entry in the i th row and j th column of the matrix D. Recall that since D is a symmetric matrix, we only need to approximate n ^ n 2 + 1 ^ elements of D. S T E P 2. Apply F - T algorithm to the vector d and integer N obtaining the integral vector d. 75 S T E P 3. Construct the integer, symmetric n x n matrix D entries of the output vector d. E N D . Theorem 5.2.1 to follow is a generalization of Theorem 3.1 in [19]. T H E O R E M 5.2.1 The matrix D satisfies the output criteria. P R O O F Using the entries of the vector d{ (resp., d'{, pi ) from F - T algorithm, construct a symmetric matrix (resp., D'{, P; ). Denote by <5t (resp., 7T; ) the largest absolute value of the entries in JDt (resp., P, ) and by r the number of iterations in F - T algorithm. The first part of the output criteria (i.e., the bound on the entries of D ) is satisfied by construction of the matrix D and since the F - T algorithm is valid. The validity of the second assertion can be shown in a similar way it was done in [19] as follows. Recall that for a matrix D its max-norm is given by || D \\oo— m a x K K n E j = i Now, we first show that (u — v)TD{(u + v) > 0 implies (u — v)TPi(u + v) > 0 . Suppose, on the contrary, that (u — v)TPi(u + v) < 0. Then, from the integrality of u,v and P t , (u — v)TPi(u + v) < —1. Therefore (u - v)TDi(u + v) = Si(u - t;)r7J^(tt + v) = 6i{u - v)T{fP{ + ^{qijy. - Pt)}(tt + v) < 6i{=± + i ( u - v)T{qtD' - Pt)(tt + v)} < < $ i { ^ + f{ || u - V ||i|| qiD\ - Pi | |oo|| « + V |]oo} < + j-n || tt - v ||oo n || qid\ - p{ |]oo|| « + v ||oo} 76 = (d{j) using the <^{^r + i^(IMIoo + IIHIoo)2} which is a contradiction. Interchanging the roles of u and v will result in a reversed inequality which in turn proves that ( i t—v) TDi(u+v) = 0 implies (u—v)TP{(u+v) = 0. From F - T algorithm and the construction of the vector d and the matrices D,Pi,...,Pr it follows that D as well as D are linearly dependent on Pi,...,Pr. Therefore, if (u — v)TP{(u -j- v) = 0 for each i , then (u — v)TD(u + v) = (u — v)TD(u + v) = 0 and in this case the theorem is proved. Now, suppose that this is not the case and that j is the smallest index such that (u — v ) T P j ( u + v) ^ 0. F - T algorithm implies that sign(u — v)TD(u + v) = sign(u — v)TDj(u + v). Since sign(w — v)TDj(u + v) is equal to sign(u — v ) T P j ( u + v), it remains to show that sign(u — v ) T P j ( u + v) = sign(u — v)TD(u + v) . To that end recall that D = J ^ = 1 Mr~lPi , where M is given in F - T algorithm. By induction on k we will prove that for j 1. i=i Now, k J2 M^iu - v ) T P i ( u + v)=MY, M * - 1 - * ^ - v)TPt{u + v) i=l i=l + (u - v)TPk[u + v)>M- 4n2 || pk ||oo > M - 2h2+hNhN > 0. 77 The last inequality follows from F - T algorithm since || a1'- \ \ o o — 1 implies || p; \\oo < Q i < 2 h 2 + h N f i . This completes the proof. • The preprocessing algorithm described above can precede, for example, an im-plicit enumeration algorithm for solving quadratic 0—1 programming problems. At each node, due to the above transformation, the continuous relaxation of the corre-sponding integer subproblem can then be solved in time independent of the objective function coefficients. Observe that we can always assume, without loss of generality, that the continuous subproblems are convex quadratic programming problems, i.e. the matrix associated with the quadratic terms is positive semidefinite (see e.g. [25]). 78 Chapter VI A R E A S F O R F U R T H E R R E S E A R C H In this thesis we extended a number of recent results for linear programming problems to quadratic programming problems. Moreover, the results from Chap-ter IV were shown to be valid for a broader class of problems, namely for nonlinear integer programming problems whose convex continuous relaxations satisfy a given constraint qualification. One possible avenue of further research is to try and extend the results obtained in Chapter III to a broader class of problems, for example to separable convex problems in which some or all of the variables are restricted to be integral. As far as Chapter V is concerned, one might try to extend the given result to quadratic integer programming problems in which the integral variables are not necessarily restricted to be 0 or 1 . In Chapter II a polynomial algorithm (whose running time is independent of the size of the linear cost coefficients and the right hand side vectors) was proposed for a class of strictly convex quadratic programming problems. It is an open question whether there exists a strongly polynomial algorithm for the above class of problems, as well as whether there exists such an algorithm for a class of linear programs. Finally, although the results in this thesis have mainly theoretical significance, one might investigate the practical benefits in some cases. For example, the calculation of 79 a lower bound on the objective function value of a problem with perturbed right hand side vector (see Chapter IV) might help a decision maker in deciding on a suitable changes of the initial right hand side vector. 80 B i b l i o g r a p h y [l] A.Bachem and B.Korte, An Algorithm for Quadratic Optimization Over Trans-portation Polytopes , Z. Angew. Math. Mech. 58 (1978) 456-461 . [2] B. Bank and R. Hansel, Stability of Mixed-Integer Quadratic Programming Problems , Mathematical Programming Study 21 (1984) 1-17 . [3] M.S. Bazaraa, J.J.Goode and C M . Shetty, Constraint Qualifications Revisited, Management Science, 18 (1972) 567-573 . [4] M.S. Bazaraa and C M . Shetty, Nonlinear Programming, Theory and Algo-rithms, John Willey, New York (1979) . [5] M . C . Biggs, Constrained Minimization Using Recursive Equality Quadratic Pro-gramming , Numerical Methods for Nonlinear Optimization, ed. F . A . Lootsma, Academic Press (1972) 411-428 . [6] J . C . G . Boot, Quadratic Programming : Algorithms, Anomalies, Applications, Studies in Mathematical and Managerial Economics, Vol 2, ed. H. Theil, North-Holland Publishing Company (1964) . [7] J.W.S. Cassels, An Introduction to the Geometry of Numbers, Die Grundlehren Der Mathematischen Wissenschaften In Einzeldarstellungen, Band 99, Springer-Verlag (1959). [8] W. Cook, A . M . H . Gerards, A. Schrijver and E . Tardos, Sensitivity Theorems in Integer Linear Programming , Mathematical Programming 34 (1986) 251-264 . [9] M . W . Cooper, Postoptimality Analysis in Nonlinear Integer Programming: The Right-Hand Side Case , Naval Res. Logist. Quart. 28 (1981) 301-307 . [10] R.W. Cottle, Symmetric Dual Quadratic Programs , Quart. Appl. Math. 21 (1963) 237-243 . [11] J .W. Daniel, Stability of the Solution of Definite Quadratic Programs , Mathe-matical Programming 5 (1973) 41-53. [12] J .L . Debiesse and G. Matignon, Comparison of Different Methods for Calcula-tion of Traffic Matrices , Ann. Telecommunic. 35 (1980) 91-102 . [13] W.S. Dorn, Duality in Quadratic Programming , Quart. Appl. Math. 18 (1960) 155-162 . [14] _, A Duality Theorem for Convex Programs, IBM J. of Res. and Devel. 4 (1960) 407-413. 81 [15] B .C . Eaves, On Quadratic Programming , Management Science Vol. 17 No. 11 (1971) 698-711 . [16] J . Edmonds, System of Distinct Representatives and Linear Algebra , J. of Res. Nat. Bur. Standards 71-B (1967) 241-245 . [17] D .T . Finkbeiner II, Introduction to Matrices and Linear Transformations, A Series of Books in the Mathematical Sciences, Third Edition, ed. V. Klee, W.H. Freedman and Company (1978) . [18] R. Fletcher, An Algorithm for Solving Linearly Constrained Optimization Prob-lems , Mathematical Programming 2 (1972) 133-165 . [19] A. Frank and E . Tardos, An Application of Simultaneous Approximation in Combinatorial Optimization , 26-th Annual Symposium on Foundations of Com-puter Science, IEEE, New York (1985) 459-463 . [20] U . M . Garcia-Palomares and O.L Mangasarian , Superlinearly Convergent Quasi-Newton Algorithms for Nonlinearly Constrained Optimization Problems , Math-ematical Programming 11 (1976) 1-13 . [21] P.E. Gill, W. Murray, M . A . Saunders and M . H . Wright, QP-Based Methods for Large-Scale Nonlinearly Constrained Optimization , Nonlinear Programming 4, ed. Mangasarian, Meyer and Robinson , Academic Press (1981) 57-98 . [22] F. Glover, Improved Linear Integer Programming Formulations of Nonlinear Integer Programs, Management Science 22 (1975) 455-460. [23] M . Grotschel, L . Lovasz and A. Schrijver, The Ellipsoid Method and Combina-torial Optimization, Springer, Heidelberg, to appear. [24] L . G . Hacijan, A Polynomial Algorithm in Linear Programming , Soviet Math. Dokl. 20 (1979) 191-194 . [25] P.L.Hammer and A . A . Rubin, Some Remarks on Quadratic Programming with 0-1 variables , Revue Francaise d'Informatique et de Recherche Operationalle 4 (1970) 67-79 . [26] P. Hansen, Quadratic Zero-One Programming by Implicit Enumeration, Nu-merical Methods in Nonlinear Optimization, ed. F .A. Lootsma, Academic Press, New York (1972) 265-278. [27] M . Held, P. Wolfe and H. Crowder, Validation of Subgradient Optimization , Mathematical Programming 6 (1974) 62-88 . [28] R. Helgason, J . Kennington and H. Lall, A Polynomially Bounded Algorithm for a Singly Constrained Quadratic Program, Mathematical Programming 18 (1980) 338-343 . 82 [29] J. Kennington and M . Shalaby, An Effective Subgradient Procedure for Minimal Cost Multy-Commodity Flow Problems , Management Science 23 (1977) 994-1004 . [30] M . K . Kozlov, S.P. Tarasov and L . G . Hacijan, Polynomial Solvability of Convex Quadratic Programming , Soviet Math. Dokl. 20 (1979) 1108-1111 . [31] D.J. Laughhunn, Quadratic Binary Programming with Applications to Capital Budgeting Problems , Operations Research 10 (1970) 454-467 . [32] R. Lazimy, Improved Algorithm for Mixed-Integer Quadratic Programs and a Computational Study , Mathematical Programming 32 (1985) 100-113 . [33] A . K . Lenstra, H.W. Lenstra,Jr and L.Lovasz, Factoring Polynomials with Ra-tional Coefficients , Math. Ann. 261 (1982) 515-534 . [34] J. Lintner, The Valuation of Risk Assets and the Selection of Risky Investments in Stock Portfolios and Capital Budgets , Rev. Econ. Statist. 47 (1965) 13-37 . [35] L. Lovasz, An Algoritmic Theory of Numbers, Graphs and Convexity , Report No.85368-OR, Institut fur Okonometrie und Operations Research, Universitat Bonn, (1985) . [36] O.L. Mangasarian, Nonlinear Programming, McGraw Hill Series in Systems Science (1969) . [37] J .C . Mao and B.A. Wallingford, An Extension of Lawler and Bell's Method of Discrete Optimization with Examples from Capital Budgeting, Management Science 15 (1969) 851-860. [38] R.D. McBride and J.S. Yormark, Finding all Solutions for a Class of Parametric Quadratic Integer Programming Problems, Management Science 26 (1980) 784-795. [39] M . Minoux, A Polynomial Algorithm for Minimum Quadratic Cost Flow Prob-lems , European Journal of Operational Research 18 (1984) 377-387 . [40] J .J. Moder and C R . Phillips, Project Management with CPM and PERT , Reinhold Ind. Eng. and Man. Sci. Texbook Series, New York (1964) . [41] M . A . Radke, Sensitivity Analysis in Discrete Optimization , Ph.D. Disserta-tion, University of California , Los Angeles (1974) . [42] L. Schrage and L . Wolsey, Sensitivity Analysis for Branch and Bound Integer Programming , Operations Research 33 (1985) 1008-1023 . [43] A. Schrijver, Theory of Linear and Integer Programming , Wiley-Interscience 83 Series in Discrete Mathematics and Optimization (1986) . [44] J . Stoer and C. Witzgall, Convexity and Optimization in Finite Dimensions I, Springer-Verlag (1970) . [45] E . Tardos, A Strongly Polynomial Algorithm to Solve Combinatorial Linear Programs , Operations Research 34 (1986) 250-256 . [46] C. Van de Panne, Methods for Linear and Quadratic Programming, Studies in Mathematical and Managerial Economics (ed. H. Theil), North-Holland Publish-ing Company (1975) . [47] L . G . Watters, Reduction of Integer Polynomial Problems to Zero-One Linear Programming Problems , Operations Research 15 (1967) 1171-1174. [48] J . Wilkinson and C. Reinsch, Linear Algebra, Handbook for Automatic Compu-tation, Vol.2 , Springer-Verlag, New York (1971) . [49] P. Wolfe, A Duality Theorem for Nonlinear Programming, Quart. Appl. Math. 19 (1961) 239-244. 84