AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM FOR LINEAR PROGRAMMING ILAN ADLER, NARENDRA KARMARKAR, MAURICIO G.C. RESENDE, AND GERALDO VEIGA A BSTRACT. This paper describes the implementation of power series dual affine scaling variants of Karmarkar’s algorithm for linear programming. Based on a continuous version of Karmarkar’s algorithm, two variants resulting from first and second order approximations of the continuous trajectory are implemented and tested. Linear programs are expressed in an inequality form, which allows for the inexact computation of the algorithm’s direction of improvement, resulting in a significant computational advantage. Implementation issues particular to this family of algorithms, such as treatment of dense columns, are discussed. The code is tested on several standard linear programming problems and compares favorably with the simplex code MINOS 4.0. 1. I NTRODUCTION We describe in this paper a family of interior point power series affine scaling algorithms based on the linear programming algorithm presented by Karmarkar (1984). Two algorithms from this family, corresponding to first and second order power series approximations, were implemented in FORTRAN over the period November 1985 to March 1986. Both are tested on several publicly available linear programming test problems (Gay 1985). We also test one of the algorithms on randomly generated multi-commodity network flow problems (Ali & Kennington 1977) and on timber harvest scheduling problems (Johnson 1986). Several authors (see, e.g., Aronson, Barr, Helgason, Kennington, Loh & Zaki (1985), Lustig (1985), Tomlin (1987) and Tone (1986)) have compared implementations of interior point algorithms with simplex method codes, but have been unable to obtain competitive solution times. An implementation of a projected Newton’s barrier method reported by Gill, Murray, Saunders, Tomlin & Wright (1986) presents the first extensive computational evidence indicating that an interior point algorithm can be comparable in speed with the simplex method. In our computational experiments, solution times for the interior point implementations are, in most cases, less than those required by MINOS 4.0 (Murtagh & Saunders 1977). Furthermore, we are typically able to achieve 8 digit accuracy in the optimal objective function value without experiencing the numerical difficulties reported in previous implementations. MINOS is a FORTRAN code intended primarily for the solution of constrained nonlinear programming problems, but includes an advanced implementation of the simplex Date: 14 March 1989 (Revised May 25 1995). Key words and phrases. Linear programming, Karmarkar’s algorithm, Interior point methods. This article originally appeared in Mathematical Programming 44 (1989), 297–335. An errata correcting the description of the power series algorithm was published in Mathematical Programming 50 (1991), 415. The version presented here incorporates these corrections to the body of the article. In addition, we updated references to some articles published in final form after the publication of this paper. We request that authors refer to this article by its original bibliographical reference. 1 2 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA method. An updated version, MINOS 5.0 (Murtagh & Saunders 1983), features a scaling option and an improved set of routines for computing and updating sparse LU factors. This latest version of MINOS was not available at the University of California, Berkeley, where the computational tests described in this paper were carried out. We believe that for the purposes of this study, MINOS 4.0 constitutes a reasonable benchmark simplex implementation. Furthermore, as evidenced by the results reported in Gill et al. (1986) we do not expect MINOS 5.0 to perform significantly faster for the test problems considered here. The plan of the paper is as follows. In Section 2, we describe our Algorithm I, a basic interior point method, commonly referred to as the affine scaling algorithm. When Algorithm I takes infinitesimal steps at each iteration, the resulting continuous trajectory is described by a system of differential equations. In Section 3, we discuss the family of algorithms constructed by truncating the Taylor expansion representing the solution of this system of differential equation. A first order approximation to the Taylor expansion results in Algorithm I. We also implement Algorithm II, which is obtained by truncating the Taylor expansion to a second order polynomial. We show in Section 4 how an initial interior solution can be obtained and in Section 5 how the algorithms can be applied to general linear programming problems. In Section 6, we describe the stopping criterion used in the computational experiments. Section 7 discusses some implementation issues, including symbolic and numerical factorizations, using an approximate scaling matrix, reducing fill-in during Gaussian elimination and using the conjugate gradient algorithm with preconditioning to solve problems with dense columns. In Section 8, we report the computational results of running the algorithms on three sets of test problems. The first set is a family of publicly available linear programs from a variety of sources. The others are, respectiv ely, randomly generated multi-commodity network flow problems and timber harvest scheduling problems. We compare these results to those of the simplex code MINOS 4.0 and for the case of the multi-commodity network flow problems to MCNF85 a specialized simplex based code for multi-commodity network flow (Kennington 1979). Conclusions are presented and future research is outlined in Section 9. 2. D ESCRIPTION OF A LGORITHM Consider the linear programming problem: (2.1) P: maximize cT x (2.2) subject to: Ax ≤ b where c and x are n-vectors, b an m-vector and A is a full rank m × n matrix, where m ≥ n and c 6= 0. In assumptions relaxed later, we require P to have an interior feasible solution x0 . Rather than expressing P in the standard equality form, we prefer the inequality formulation (2.1) and (2.2). As discussed later in this section, there are some computational advantages in the selection of search directions under this formulation. In Section 5, we show how to apply this approach to linear programs in standard form. The algorithm described below is a variation of Karmarkar’s original projective algorithm (Karmarkar 1984), substituting an affine transformation for the projective transformation, and the objective function for the potential function. Recently, it came to our attention that this algorithm was first proposed independently by Dikin (1967). Similar algorithms have been discussed by Barnes (1986) and Vanderbei, Meketon & Freedman (1986). They, however, express the linear programming problem in (standard) equality form. Following the taxonomy referred to in Hooker (1986), the algorithms presented in this paper can be AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 3 classified as dual affine scaling algorithms. However, these algorithms are equivalent to their primal counterparts applied to problems in inequality form. Starting at x0 , the algorithm generates a sequence of feasible interior points {x1 , x2 , . . . , xk , . . . } with monotonically increasing objective values, i.e., (2.3) b − Axk > 0 (2.4) cT xk+1 > cT xk , terminating when a stopping criterion to be discussed later is satisfied. Introducing slack variables to the formulation of P, we have: (2.5) P: maximize cT x (2.6) subject to: Ax + v = b (2.7) v ≥ 0, where v is the m-vector of slack variables. Affine variants of Karmarkar’s algorithm consist of a scaling operation applied to P, followed by a search procedure that determines the next iterate. At each iteration k, with vk and xk as the current iterates, a linear transformation is applied to the solution space, (2.8) v̂ = D−1 v v, where (2.9) Dv = diag(vk1 , . . . , vkm ) . The slack variables are scaled so that xk is equidistant to all hyperplanes generating the closed half-spaces whose intersection forms the transformed feasible polyhedral set, (2.10) −1 {x ∈ ℜn | D−1 v Ax ≤ Dv b} . Rewriting equations (2.1) and (2.2) in terms of the scaled slack variables, we have (2.11) P̂: maximize cT x (2.12) subject to Ax + Dv v̂ = b (2.13) v≥0 The set of feasible solutions for P is denoted by (2.14) X = {x ∈ ℜn | Ax ≤ b} and the set of feasible scaled slacks in P̂ is (2.15) V̂ = {v̂ ∈ ℜm | ∃x ∈ X, Ax + Dv v̂ = b} . As observed in Gonzaga (1991), under the full-rank assumption for A, there is a one to one relationship between X and V̂, with (2.16) v̂(x) = D−1 v (b − Ax) and (2.17) −1 T −1 −1 x(v̂) = (AT D−2 v A) A Dv (Dv b − v̂) . There is also a corresponding one to one relationship linking feasible directions hx in X and hv̂ in V̂, with (2.18) hv̂ = −D−1 v Ahx 4 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA and (2.19) −1 T −1 hx = −(AT D−2 v A) A Dv hv̂ . Observe from (2.18) that a feasible direction in V̂ lies on the range space of D−1 v A. As in other presentations of affine variants of Karmarkar’s algorithm, the search direction selected at each iteration is the projected objective function gradient with respect to the scaled variables. Since only the slack variables are changed by the affine transformation, using (2.11) and (2.17), we can compute the gradient of the objective function with respect to v̂, (2.20) T −2 −1 ∇v̂ c(x(v̂)) = (∇v̂ x(v̂))T ∇x c(x) = −D−1 v A(A Dv A) c . The gradient with respect to the scaled slacks lies on the range space of D−1 v A, making a projection unnecessary. Consequently, the search direction in v̂ is (2.21) T −2 −1 hv̂ = −D−1 v A(A Dv A) c , and from (2.19) the corresponding feasible direction in X can be computed, (2.22) −1 hx = (AT D−2 v A) c . Applying the inverse affine transformation to hv̂ , we obtain the corresponding feasible direction for the unscaled slacks, (2.23) −1 hv = −A(AT D−2 v A) c . Under the assumption that c 6= 0, unboundedness is detected if hv ≥ 0. Otherwise, the next iterate is computed by taking the maximum feasible step in the direction hv , and retracting back to an interior point according to a safety factor γ, 0 < γ < 1, i.e., (2.24) xk+1 = xk + αhx , where (2.25) α = γ × min{−vki /(hv )i | (hv )i < 0, i = 1, . . . , m} . The above formulation allows for inexact projections without loss of feasibility. Even if hx is not computed exactly in (2.22), a pair of feasible search direction can still be obtained by computing (2.26) hv = −Ahx . Given an interior feasible solution x0 , a stopping criterion and a safety factor γ, Pseudocode 2.1 describes Algorithm I as outlined in this section. As in subsequent instances in this paper, algorithms are expressed in an Algol-like algorithmic notation described by Tarjan (1983). 3. A FAMILY OF I NTERIOR -P OINT A LGORITHMS U SING P OWER S ERIES A PPROXIMATIONS Consider the continuous trajectory generated by Algorithm I when infinitesimal steps are taken at each iteration (Bayer & Lagarias 1989). Let us denote the path of interior feasible solutions for P by (x̄(τ), v̄(τ)), where the real parameter τ is the continuous counterpart to the iteration counter. For any value of τ, the corresponding search directions hx̄ (τ) and hv̄ (τ) can be computed by expressions (2.22) and (2.23). Alternatively, the search directions can be described by an equivalent system of linear equations, (3.1) −AT D−2 v̄(τ) hv̄(τ) = c AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 5 procedure Algorithm I (A, b, c, x0 , stopping criterion, γ) k := 0; do stopping criterion not satisfied → vk := b − Axk ; Dv := diag(vk1 , . . . , vkm ); −1 hx := (AT D−2 v A) c; hv := −Ahx ; if hv ≥ 0 → return fi; α := γ × min{−vki /(hv )i | (hv )i < 0, i = 1, . . . , m}; xk+1 := xk + αhx ; k := k + 1; od end Algorithm I Pseudo-code 2.1: Algorithm I Ahx̄(τ) + hv̄(τ) = 0 . (3.2) By taking infinitesimal steps, the resulting continuous trajectory is such that d v̄ d x̄ (τ) = hx̄(τ) and (τ) = hv̄(τ) . dτ dτ Given the system of linear equations (3.1) and (3.2), we replace the search directions by the corresponding derivatives with respect to the trajectory parameter, resulting in the system of nonlinear, first order differential equations (3.3) −AT D−2 v̄(τ) (3.4) (3.5) A d v̄ (τ) = c dτ d v̄ d x̄ (τ) + (τ) = 0 , dτ dτ with the boundary conditions (3.6) x̄(0) = x0 and v̄(0) = v0 , where the initial solution (x0 , v0 ) is given and satisfies Ax0 + v0 = b, v0 > 0. By following the trajectory satisfying (3.4)–(3.6) we can, theoretically, obtain the optimal solution to P (Adler & Monteiro 1991). In practice, we build an iterative procedure where, after replacing the current iterate for the initial solution in the boundary condition, an approximate solution to (3.4)–(3.6) is computed, using a truncated Taylor power series expansion. The next iterate is determined through a search on the approximate trajectory. As discussed later in this section, the search procedure requires a suitable reparametrization of the continuous trajectory, τ = ρ(t) , (3.7) where ρ(t) is a monotonically increasing, infinitely differentiable real function such that (3.8) x(t) = x̄(ρ(t)) and v(t) = v̄(ρ(t)). Reparametrizing (3.4)–(3.6), we have (3.9) −AT D−2 v(t) dρ dv (t) = (t)c dt dt 6 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA A (3.10) dx dv (t) + (t) = 0 , dt dt and the boundary conditions x(0) = x0 (3.11) and v(0) = v0 . Since we cannot compute an exact solution to (3.9)–(3.11), we use approximate solutions to the system of differential equations as the backbone of a family of iterative algorithms. At each iteration k, the algorithm restarts the trajectory with an initial solution (x(0), v(0)) = (xk−1 , vk−1 ). A new iterate is generated by moving on the approximate trajectory without violating nonnegativity of the slack variables. The approximate solutions can be computed by means of a truncated Taylor expansion of order r (Karmarkar, Lagarias, Slutsman & Wang 1989), such that for t > 0, r ti dix (0) i i=1 i! dt x(t) ≈ x̃(t) = x(0) + ∑ (3.12) and r ti div (0) . i i=1 i! dt v(t) ≈ ṽ(t) = v(0) + ∑ (3.13) From (3.9) and (3.10), we can compute derivatives of all orders for x(t) and v(t). Consider the functions F(t) = −D−1 v(t) (3.14) and G(t) = Dv(t) , which satisfy F(t)G(t) = −I . (3.15) Applying Leibniz’s differentiation theorem to the product F(t)G(t), we have i (3.16) d i− j F i! d jG ∑ (i − j)! j! dt i− j (t) dt j (t) = 0 for i ≥ 1 . j=0 From (3.16), the derivatives of F(t) can be computed recursively by (3.17) i−1 i diF d i− j F d j G i! −2 d G −1 (t) = D (t) − D (t) j (t) ∑ v(t) dt i v(t) i− j dt i dt j=1 (i − j)! j! dt for i ≥ 1 . Taking derivatives of both sides of (3.9), and rewriting the left hand side of the equation in terms of F(t), we have d iF d iρ (t)e = (t)c for i ≥ 1 . dt i dt i By combining (3.17) with (3.18), and annexing the appropriate derivative of (3.10), we have a system of linear equations that recursively computes the derivatives of x(t) and v(t) evaluated at t = 0, −AT (3.18) (3.19) (3.20) −AT D−2 v(0) i−1 d i− j F i! diρ d jv d iv T −1 (0) = (0)c + A D (0) (0) v(0) ∑ (i − j)! j! dt i− j dt i dt i dt j j=1 A d ix div (0) + i (0) = 0 . i dt dt AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 7 In an implementable reformulation of (3.19) and (3.20), we eliminate the binomial terms by defining (3.21) Fi = 1 d iF (0) , i! dt i (3.22) ρi = 1 diρ (0) , i! dt i 1 div 1 d ix i (0) and z = (0) . v i! dt i i! dt i Solving (3.19) and (3.20), we have (3.23) zix = i−1 (3.24) j −1 T −2 −1 T −1 zix = ρi (AT D−2 v(0) A) c + (A Dv(0) A) A Dv(0) ∑ Fi− j zv j=1 and ziv = −Azix . (3.25) The approximate trajectory based on a power series of order r is rewritten as r (3.26) x̃(t) = x(0) + ∑ t i zix i=1 r and ṽ(t) = v(0) + ∑ t i ziv . i=1 At iteration k, the next iterate is computed as (3.27) xk+1 = x̃(α) and vk+1 = ṽ(α) , with (3.28) α = γ × sup{t | 0 ≤ t ≤ 1, ṽi (t) ≥ 0, i = 1, . . . , m} , where 0 < γ < 1. By selecting a suitable reparametrization ρ(t), the line search on ṽ(t) that computes α in (3.28) can be limited to the interval 0 ≤ t ≤ 1. Since x̃(t) and ṽ(t) depend solely on derivatives of ρ(t) evaluated at t = 0, the desired reparametrization is fully characterized by constants ρ1 , . . . , ρr . Values for ρi are computed by selecting a row index l, and forcing ṽl (1) = 0 (3.29) and (3.30) d i ṽl (0) = 0 for i = 2, . . . , r . dt i Consider the search directions −1 h1x = (AT D−2 v(0) A) c (3.31) and h1v = −Ah1x . (3.32) From (3.24) and (3.25), we have (3.33) z1x = ρ1 h1x and z1v = ρ1 h1v . 8 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA Once again, under the assumptions that A is full rank and c 6= 0, unboundedness is detected if hv ≥ 0. Otherwise, we determine the row index l by performing a ratio test between v(0) and search direction h1v , i.e., l = argmin{−v(0)/(h1v )i | (h1v )i < 0} . (3.34) 1≤i≤m This operation corresponds to searching along the first order approximation of trajectory v(t). From the conditions imposed on the reparametrization by (3.30), the truncated Taylor expansion in (3.26) is such that ṽl (1) = v(0) + ρ1h1v , (3.35) and from (3.29) and (3.33), we compute ρ1 = −vl (0)/(h1v )l . (3.36) Iteratively, for i = 2, . . . , r, we compute ρi , zix and ziv . Based on (3.24) and (3.25), zix = ρi h1x + hix (3.37) and ziv = −Azix . where i−1 (3.38) j −1 T −1 hix = (AT D−2 v(0) A) A Dv(0) ∑ Fi− j zv and hiv = −Ahix . j=1 To satisfy (3.30), we have (3.39) ρi = −(hiv )l /(h1v )l . Pseudo-code 3.1 formalizes the algorithm based on a truncated power series of order r. The computational experiments reported in Section 8 include Algorithm II, which is the second order version of this algorithm. In a practical implementation, the additional computational effort, when compared to Algorithm I, is dominated by the solution of systems of symmetric linear equations to compute hix for i = 2, . . . , r. As in Algorithm I, if x̃(t) is not obtained exactly, a pair of feasible trajectories can still immediately available. By computing (3.40) ziv = −Azix for i = 1, . . . , r , we guarantee that (x̃(t), ṽ(t)) is a feasible for 0 ≤ t ≤ α, where 0 ≤ α ≤ 1 is the maximum feasible step. 4. I NITIAL S OLUTION Algorithm I and the truncated power series algorithms require that an initial interior solution x0 be provided. Since such solution does not necessarily exist for a generic linear program, and a starting solution close to the faces of the feasible polyhedral set can imply in very slow convergence (Megiddo & Shub 1989), we propose a Phase I/Phase II scheme, where we first solve an artificial problem with a single artificial variable having a large cost coefficient assigned to it. Firstly, we compute a tentative solution for P, (4.1) x0 = (kbk2 /kAck2)c . For the computational experiments reported in Section 8, we compute the initial value for the artificial variable as (4.2) x0a = −2 × min{(b − Ax0)i | i = 1, . . . , m} . AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 22. 23. 24. 25. 9 procedure PowerSeries(A, b, c, x0 , stopping criterion, γ, r) k := 0; do stopping criterion not satisfied → x(0) := xk ; v(0) := b − Ax(0); Dv (0) := diag(v1 (0), . . . , vm (0)); h1x := (AT Dv (0)−2 A)−1 c; h1v := −Ah1x ; if h1v ≥ 0 return fi; l := argmin1≤i≤m {−v(0)/(h1v )i | (h1v )i < 0}; ρ1 := −vl (0)/(h1v )l ; z1x := ρ1 h1x ; z1v := ρ1 h1v ; F1 := Dv (0)−2 z1v ; for i = 2, . . . , r → j hix := (AT Dv (0)−2 A)−1 AT Dv (0)−1 ∑i−1 j=1 Fi− j zv ; hiv := −Ahix ; ρi := −(hiv )l /(h1v )l ; zix := ρi h1x + hix ; ziv := −Azix ; j=1 Fi := Dv (0)−2 ziv − Dv (0)−1 ∑i−1 Fi− j zvj ; rof; x̃(t) := x(0) + ∑ri=1 t i zix ; ṽ(t) := v(0) + ∑ri=1 t i ziv ; α := γ × sup{t | 0 ≤ t ≤ 1, ṽi (t) ≥ 0, i = 1, . . . , m}; xk+1 := x̃(α); k := k + 1; od end PowerSeries Pseudo-code 3.1: Power Series Algorithm Although it did not occur in any of the test problems reported in this paper, the computation in (4.2) can be such that x0a ≈ 0. Consequently, for application to generic linear programs, we recommend an alternative computation of the initial value of the artificial variable, (4.3) x0a = 2 × kb − Ax0k2 . The n + 1-vector (x0 , x0a ) is an interior solution of the Phase I linear programming problem (4.4) Pa : maximize cT x − Mxa (4.5) subject to: Ax − exa ≤ b where (4.6) e = (1, 1, . . . , 1)T . The large artificial cost coefficient is computed as a function of the problem data, (4.7) M = µ × cT x0 /x0a , where µ is a large constant. Initially, the algorithm is applied to Pa with a modified stopping criterion. In this Phase I stage, the algorithm either identifies an interior dual feasible solution, or, if no such solution 10 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA exists, finds a solution that satisfies the stopping criterion for problem P. With ε f defined as the feasibility tolerance, the modified stopping criterion for Phase I is formulated as follows: (i) If xka < 0 at iteration k, then xk is an interior feasible solution for problem P. (ii) If the algorithm satisfies the regular stopping criterion and xka > ε f , P is declared infeasible. (iii) If the algorithm satisfies the regular stopping criterion and xka < ε f , either unboundedness is detected or an optimal solution is found. In this case, P has no interior feasible solution. If Phase I terminates according to condition (i), the algorithm is applied to problem P starting with the last iterate of Phase I and using the regular stopping criterion. 5. A PPLICATION TO THE G ENERAL L INEAR P ROGRAMMING P ROBLEM It is not common practice to formulate a linear programming problem as in (2.1) and (2.2). Instead, a standard form is usually preferred, (5.1) LPP: minimize bT y (5.2) subject to: AT y = c (5.3) y≥0 where A is an m × n matrix, c an n-vector and b and y m-vectors. The dual linear programming problem of (5.1)–(5.3), however, has the desired form: (5.4) LPD: maximize cT x (5.5) subject to: Ax ≤ b where x is an n-vector. Note that LPD is identical to P, as defined in (2.1) and (2.2). At iteration k, with current solutions xk and vk , a tentative dual solution is defined as (5.6) T −2 −1 yk = D−2 v A(A Dv A) c , where (5.7) Dv = diag(vk1 , . . . , vkm ) . This computation of the tentative dual solution, similar to the one suggested by Todd & Burrell (1986), can be performed by scaling the first order search direction in each iteration of the algorithm. From (2.23), for iteration k, (5.8) yk = D−2 v hv . The tentative dual solution minimizes the deviation from complementary slackness with respect to the current iterate xk , relaxing the nonnegativity constraints (5.3) (Chandru & Kochar 1985). Formally, consider the problem m (5.9) minimize ∑ (vkj y j )2 j=1 (5.10) subject to: AT y = c . The tentative dual solution in (5.6) is the solution to the Karush-Kuhn-Tucker stationary conditions of minimization problem given by (5.9) and (5.10). Under nondegeneracy, given AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 11 a sequence of feasible primal solutions converging to an optimal solution, the corresponding sequence of tentative dual estimates also converges to the optimal dual solution. 6. S TOPPING C RITERION In the computational experiments reported in this study, both Algorithms I and II are terminated whenever the relative improvement to the objective function is small, i.e., (6.1) |cT xk − cT xk−1 |/ max{1, |cT xk−1 |} < ε , where ε is a given small positive tolerance. The tentative dual solution yk , computed in (5.6), can be used to build an alternative stopping criterion. If yk and vk satisfy (6.2) AT yk = c , (6.3) yk ≥ 0 , (6.4) ykj vkj = 0, j = 1, 2, . . . , m , then, by duality theory of linear programming, yk is optimal for LPP and xk is optimal for LPD. Since (6.2) is automatically satisfied, a stopping criterion, replacing (6.1), is such that the algorithm terminates whenever (6.5) ykj ≥ −ε1 kyk k2 , j = 1, 2, . . . , m , and (6.6) |ykj vkj | ≤ ε2 kyk k2 kvk k2 , j = 1, 2, . . . , m , for given small positive tolerances ε1 and ε2 . The relations in (6.5) and (6.6) serve as a verification that xk and yk are indeed approximate optimal solutions for LPD and LPP, respectively. Unboundedness of LPD is, theoretically, detected by the algorithm whenever hv ≥ 0, i.e., the ratio tests in (2.25) and (3.34), involving the first order search direction fail. In practice, an additional test is required, LPD is declared unbounded whenever the objective function value exceeds a supplied bound. 7. I MPLEMENTATION I SSUES In this section, we briefly discuss some important characteristics of this implementation. A detailed description of the data structures, programming techniques and used in the implementation of Algorithms I and II is the subject of Adler, Karmarkar, Resende & Veiga (1989). 7.1. Computing Search Directions. As in other variants of Karmarkar’s algorithm, the main computational requirement of the algorithms described in this paper consists of the solution of a sequence of sparse symmetric positive definite systems of linear equations determining the search directions for each iteration. In Algorithm I, for each iteration k, a linear trajectory is determined by the feasible direction computed in (7.1.1) Bk h x = c , with (7.1.2) Bk = AT D−2 v A, 12 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA where (7.1.3) Dv = diag(vk1 , . . . , vkm ) . Under the full-rank assumption for A, the system of linear equations in (7.1.1) is symmetric and positive definite. Such systems of linear equations are usually solved by means of the LU factorization, (7.1.4) Bk = LU , where L is an m × m unit lower triangular matrix (a matrix that has exclusively ones in its main diagonal), and U is an m × m upper-triangular matrix. In the case of positive definite matrices, this LU factorization always exists and is unique. Furthermore, if the system is also symmetric, the factor L can be trivially obtained from U with (7.1.5) L = DU−1U T , where DU is a diagonal matrix with elements drawn from the diagonal of U. Rewriting (7.1.1), we have (7.1.6) (LU)hx = c . Direction hx in (7.1.6) can be determined by solving two triangular systems of linear equations, performing a forward substitution (7.1.7) Lz = c , followed by a back substitution (7.1.8) Uhx = z . In each iteration of the higher order algorithms described in Section 3, the additional directions in (3.38) also involve the solution of linear systems identical to (7.1.1), except for different right-hand sides. Solving these additional systems of linear equations involve only the back and forward substitution operations, using the same LU factors. Furthermore, during the execution of the algorithm only Dv changes in (7.1.2), while A remain unaltered. Therefore, the nonzero structure of AT D−2 v A is static throughout the entire solution procedure. An efficient implementation of the Gaussian elimination procedure can take advantage of this property by performing, in the beginning of the algorithm, a single symbolic factorization step, i.e., operations that depend solely on the nonzero structure of the system matrix. For example, at this stage of the algorithm, we determine the nonzero structure of LU factors and build a list of the numerical operations performed during the Gaussian elimination procedure. At each iteration k of the algorithm, the actual numerical values of Bk are computed and incorporated in the symbolic information. Next, the numerical factorization is executed, by traversing the list of operations, yielding the LU factors. 7.2. Using an Approximate Scaling Matrix. Based on the theoretical approach suggested by Karmarkar (1984), a significant reduction in the computational effort can be achieved by using an approximate scaling matrix when computing the numerical values for matrix Bk at each iteration. The matrix for the system of linear equations in (7.1.1) is replaced by (7.2.1) B̃k = AT D̃−2 k A, where D̃k is an approximate scaling matrix, computed by selectively updating the scaling matrix used in the preceding iteration. The approximate scaling matrix at each iteration k AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 13 is computed as follows: (7.2.2) ( D̃k−1 (i, i) D̃k (i, i) = Dv (i, i) if |Dv (i, i) − D̃k−1 (i, i)|/D̃k−1 (i, i) < δ if |Dv (i, i) − D̃k−1 (i, i)|/D̃k−1 (i, i) ≥ δ for a given δ > 0. After computing the current approximate scaling matrix according to 7.2.2, we define (7.2.3) −2 ∆ = D̃−2 k − D̃k−1 . Then, we have the update expression (7.2.4) B̃k = B̃k−1 + AT ∆A. This enables us to update the linear system matrix by rescaling only a reduced set of columns of AT . 7.3. Ordering for Sparsity. When a sparse matrix is factored, fill-in usually occurs. The triangular factors contain nonzero elements in positions where Bk has zeros. Fill-in degrades the performance of the sparse Gaussian elimination used to compute the LU factorization, also affecting the back and forward substitution operations. It is possible to reduce fill-in by performing a permutation of columns and rows in Bk . If P is a permutation matrix, then 7.1.1 and (7.3.1) (PBk PT )Phx = Pc are equivalent systems. Furthermore, there exists a permutation matrix P̄ such that the fillin generated in the triangular factors is minimized. Unfortunately, finding this permutation matrix is an NP-complete problem (Yannakakis 1981). However, the minimum degree and minimum local fill-in ordering heuristics (Rose 1972) have been shown to perform well in practice (Duff, Erisman & Reid 1986). We use either the minimum degree heuristic as implemented in subroutine MD of the Yale Sparse Matrix Package (Eisenstat, Gurshy, Schultz & Sherman 1982), or the minimum local fill-in implementation in deCarvalho (1987). 7.4. Treating Dense Columns in the Coefficient Matrix. In the presence of a few dense columns in AT , AT D−2 v A will be impracticably dense, regardless of the permutation matrix. Consequently, we face prohibitively high computational effort and storage requirements during Gaussian elimination. To remedy this situation, we make use of a hybrid scheme in which we first perform an incomplete factorization of AT D−2 v A. Next, we use the incomplete Cholesky factors as preconditioners for a conjugate gradient method to solve the system of linear equations defined in 7.1.1. Let (N, N̄) be a partition of the column indices of AT , such that the columns of ATN have density smaller than a given parameter λ. At iteration k, the incomplete Cholesky factors L̃k and L̃Tk are such that (7.4.1) T ATN D−2 v AN = L̃k L̃k . Using the conjugate gradient algorithm, we solve the system of linear equations (7.4.2) Qu = f , where (7.4.3) T −2 T −1 Q = L̃−1 , k (A Dv A)(L̃k ) and (7.4.4) f = L̃−1 k c. 14 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA The search direction hx is computed by performing a back substitution operation, solving (7.4.5) u = L̃Tk hx . Given a termination tolerance εcg > 0, the conjugate gradient algorithm is outlined in Pseudo-code 7.4.1. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. 11. 12. 13. 14. 15. procedure ConjugateGradient(Q, f , εcg ) u0 := f ; r0 := Qu0 − f ; p0 := −r0 ; i := 0; do kri k22 ≥ εcg → qi := Qpi ; αi := kri k22 /pTi qi ; ui+1 := ui + αi pi ; ri+1 := Qui+1 − f ; βi := kri+1 k22 /kri k22 ; pi+1 := −ri+1 + βi pi ; i := i + 1; od end ConjugateGradient Pseudo-code 7.4.1: Conjugate Gradient Algorithm 8. T EST P ROBLEMS In this section, we report the computational results of running our implementations of Algorithms I and II on a set of linear programming test problems (Gay 1985) available through NETLIB (Dongarra & Grosse 1987). NETLIB is a system designed to provide efficient distribution of public domain software to the scientific community through computer networks, e.g. Arpanet, UNIX UUCP network, CSNET Telenet and BITNET. We also report the results of running Algorithm I on linear programs generated by two models. The first set is composed of randomly generated multi-commodity network flow problems generated with MNETGN (Ali & Kennington 1977). The other is a collection of timber harvest scheduling problems generated by FORPLAN (Johnson 1986), a linear programming based system used for long-range planning by the US Forest Service. 8.1. The NETLIB Test Problems. This collection of linear programs consists of problems contributed by a variety of sources. It includes linear programming test problems from the Systems Optimization Laboratory at Stanford University, staircase linear programs generated by Ho & Loute (1981), and many real world problems. We considered in our tests all of the available problems that do not have a BOUNDS section in their MPS representation, since the current version of our code cannot handle bounded variables implicitly. Table 8.1.1 presents the statistics for 31 test problems obtained from NETLIB ordered by increasing number of nonzero elements in the coefficient matrix A, after adding slack variables. The dimensions of the problems include null or fixed variables, sometimes removed by a preprocessor before applying a linear programming algorithm. The number of rows does not include the objective function. In column 4, the number of nonzero elements of A is displayed. Column 5 gives the number of nonzero elements of the lower triangular portion of the permuted (PA)T PA matrix, using the minimum degree ordering heuristic. AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 15 TABLE 8.1.1. NETLIB test problem statistics Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israel BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47b Ship08l Ship12l a Using Rows 27 56 129 205 96 117 388 471 300 220 77 174 305 330 223 490 173 147 402 660 402 778 1090 990 1151 397 1480 929 821 778 1151 Columns 51 138 185 317 162 253 466 671 660 303 760 316 472 600 472 1275 295 1350 1506 1200 2166 2467 2500 1800 2869 2750 3340 3340 1876 4363 5533 nonz(A) 102 424 465 665 777 1179 1534 1725 1872 2202 2388 2443 2494 2732 2768 3288 3408 4316 4400 5469 6380 7194 7334 8206 8284 8584 9534 10043 10705 12882 16276 nonz((PA)T PA)a 90 377 606 654 871 967 1915 2370 1686 2190 1133 3545 2929 3143 2683 1953 1720 2099 2827 6306 4147 3552 6595 9469 4233 4280 8866 6265 10919 6772 8959 nonz(L) 107 404 734 1182 1026 1425 2324 2948 2667 2850 1392 3707 4114 4963 3416 5134 1727 2545 3134 9791 4384 4112 14870 14619 5063 5879 19469 6655 33498 7128 9501 the minimum degree ordering heuristic. is also known as BP822. b 25FV47 Column 6 gives the nonzero elements in the Cholesky factor. The fill-in is the difference between columns 6 and 5. The entries in columns 5 and 6 for problem Israel are those used by the algorithm, i.e., after 6 dense columns of AT are dropped. All test runs, with the exception of those of Algorithm II, were carried out on the IBM 3090 at U.C. Berkeley under VM/SP CMS Algorithm II test runs were conducted on the IBM 3081K at U.C.-Berkeley, also under VM/SP CMS FORTRAN programs were compiled using 16 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA the FORTVS compiler with compiler options OPT(3), NOSYM and NOSDUMP The reported CPU times were computed with utility routine DATETM taking into account preprocessing (e.g. input matrix cleanup, ordering and symbolic factorization) and all operations until termination. However, they exclude the effort required to translate the linear programs from MPS format. Execution times for Algorithm I and Algorithm II are compared to those of the simplex code MINOS 4.0 The reported CPU times for MINOS are those of subroutine DRIVER which also excludes the time to read and translate the MPS input file. The results of solving the NETLIB test problems by MINOS 4.0 are displayed in Table 8.1.2. In running MINOS we use default parameters except for LOG FREQ= 200 and ITERATIONS= 7500. The problem size parameters, ROWS, COLUMNS and ELEMENTS were set to the exact number of rows, columns and nonzero elements of the coefficient matrix. In this fashion, MINOS selects the value for the all important PARTIAL PRICE parameter according to its built-in default strategy. Using the default parameters, MINOS 4.0 achieved 7 digit accuracy on all problems, if compared to the optimal objective values reported in Gay (1985). Since the NETLIB test problems are not in the form described by 2.1 and 2.2, Algorithms I and II solve the corresponding dual linear programming problems. A primal optimal solution is obtained from 5.6, involving the last scaling matrix Dv and the last search direction for the dual slacks. By contrast MINOS solves the test problems is its original form. Tables 8.1.3 and 8.1.4 summarize computational results for Algorithm I on the set of test problems, using, respectively, the minimum degree and the minimum local fill-in ordering heuristics. Table 8.1.5 presents results for Algorithm II, which uses the minimum degree ordering heuristic. Both Algorithms I and II were tested with the following parameter settings. The safety factor parameter was set to γ = 0.99 for the first 10 iterations and γ = 0.95 thereafter. Both algorithms were terminated when the relative improvement in the objective function fell below ε = 10−8 . In Phase I, the value of the artificial variable cost defined in 4.7 was determined by the constant µ = 105 . The diagonal update tolerance was set to δ = 0.1 and the column density parameter for building the incomplete Cholesky factorization in the conjugate gradient algorithm was set to λ = 0.3. Table 8.1.6 compares execution times for the three algorithms. CPU times displayed for MINOS and Algorithm I were measured on the IBM 3090 while those for Algorithm II are from the IBM 3081-K For the sake of consistency, when computing the CPU times ratio between MINOS and Algorithm II, the execution times for MINOS are those of the IBM 3081-K Figures 8.1.1 and 8.1.2 illustrate graphically the performances of MINOS and Algorithm I. Table 8.1.7 focuses on five subgroups of problems, where each subgroup is made up of problems having the same structure, or generated by the same model. In this table, ratios of MINOS iterations to Algorithm I iterations, MINOS CPU time per iteration to Algorithm I CPU time per iteration and MINOS total CPU time to Algorithm I total CPU time are given. Table 8.1.8 presents detailed results related to the generation of primal solutions in Algorithm I. Let t be the number of iterates generated by the algorithm. The normalized duality gap (8.1.1) |bT yt − cT xt |/|bt yt | is given in column 2. Column 3 presents the maximum normalized primal infeasibility (8.1.2) max{|ATi yt − ci |/kck2 | i = 1, 2, . . . , n} . AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 40 17 b 35 MINOS Algorithm I b b r 30 b 25 b CPU Time (seconds) 20 b 15 10 5 0 b b b b b b r b b br bb rbr r b brr br rb r br rbrb brrb r r 0 2000 b b rr 4000 b b r r r r r r r r rr r 6000 8000 10000 12000 14000 16000 Number of Nonzero elements F IGURE 8.1.1. MINOS 4.0 and Algorithm I(CPU times, IBM 3090) — NETLIB test problems (minimum local fill-in ordering heuristic in Algorithm I, 25FV47 is not included in this graph). b 10 b 8 b b b Ratio b b b 6 b b b b 2 bb b b 0 b b 4 b 0 b b b b b b b b b b b 2000 b 4000 6000 8000 10000 12000 14000 16000 Number of Nonzero elements F IGURE 8.1.2. MINOS 4.0 Algorithm I (CPU times ratio, IBM 3090) — NETLIB test problems (minimum local fill-in ordering heuristic in Algorithm I) 18 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.1.2. MINOS 4.0 test statistics (IBM 3090) — NETLIB test problems) Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israel BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47 Ship08l Ship12l Phase I Itrs 2 28 38 30 83 238 146 228 184 185 97 72 207 238 69 135 9 189 309 432 441 470 667 646 561 585 696 1035 1262 594 782 Total Itrs 6 119 81 111 119 391 215 614 305 313 261 281 472 437 568 523 41 576 395 792 593 657 1121 1196 718 1590 1273 1955 7157 958 1268 Time (secs.) 0.01 0.23 0.29 0.45 0.26 1.27 1.06 4.54 1.98 1.34 1.02 1.54 2.54 2.04 3.68 5.64 0.18 2.69 3.54 8.18 7.18 9.16 22.02 19.31 12.78 14.53 32.87 39.08 217.67 14.72 26.10 Objective Value -4.6475319e+02 2.2549495e+05 -2.3313897e+06 -5.2201997e+01 -4.1573247e+02 -7.6589255e+04 1.8781247e+03 -1.4753432e+07 1.4122500e+03 1.5185099e+03 8.6666659e+00 -8.9664483e+05 -1.5862800e+02 1.8416758e+04 -1.8638928e+01 9.0429701e+02 3.3592486e+04 5.0500000e+01 1.7987147e+06 3.6660260e+04 1.7933246e+06 1.9200982e+06 1.7248071e+03 5.4901252e+04 1.4892361e+06 9.0499997e+02 1.4240000e+03 2.1851968e+06 5.5018494e+03 1.9090552e+06 1.4701879e+06 Column 4 gives the minimum normalized primal entry, (8.1.3) min{yti /kyt k2 | i = 1, 2, . . . , m} and column 5 the maximum normalized complementarity violation (8.1.4) max{|yti vti |/kyt k2 kvt k2 | i = 1, 2, . . . , m} . From the results above, we make the following observations: AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 19 TABLE 8.1.3. Algorithm I test statistics (IBM 3090) — NETLIB test problems (minimum degree ordering heuristic in Algorithm I) Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israela BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47 Ship08l Ship12l a The Phase I Itrs 1 1 3 4 4 7 5 3 6 38 2 8 7 33 40 39 23 2 5 38 5 5 7 37 4 2 7 3 44 5 5 Total Itrs 20 24 25 29 28 39 25 28 34 38 19 38 33 33 40 39 23 22 31 38 31 34 33 37 35 24 36 46 55 34 34 Time (secs.) 0.04 0.12 0.17 0.28 0.29 0.58 0.51 0.69 0.85 1.52 0.41 4.00 1.54 1.97 1.56 2.30 0.69 0.83 1.26 4.38 2.31 1.58 6.71 6.66 1.91 1.68 9.44 9.76 46.17 4.00 4.89 Objective Value -4.6475315e+02 2.2549496e+05 -2.3313898e+06 -5.2202061e+01 -4.1573224e+02 -7.6589319e+04 1.8781248e+03 -1.4753433e+07 1.4122488e+03 1.5185099e+03 8.6666666e+00 -8.9664483e+05 -1.5862802e+02 1.8416759e+04 -1.8751929e+01 9.0429695e+02 3.3592486e+04 5.0500000e+01 1.7987145e+06 3.6660261e+04 1.7933245e+06 1.9200982e+06 1.7248068e+03 5.4901254e+04 1.4892361e+06 9.0500000e+02 1.4240000e+03 2.1851927e+06 5.5018494e+03 1.9090552e+06 1.4701879e+06 conjugate gradient algorithm was triggered when running this test problem. • Iterations for Algorithm I vary from 19 to 55 (see Tables 8.1.3 and 8.1.4), growing slowly with problem size. • Algorithm I is, in general, faster than MINOS with a speed-up of up to 10.7 (see Table 8.1.6 and Figures 8.1.1 and 8.1.2). The total problem set execution time was 5.14 times faster. MINOS was faster on 5 small problems, all having less than 225 rows and 3500 nonzero matrix elements. 20 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.1.4. Algorithm I test statistics (IBM 3090) — NETLIB test problems (minimum local fill-in ordering heuristic) Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israela BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47 Ship08l Ship12l a The Phase I Itrs 1 1 3 4 4 7 5 3 6 38 2 8 7 33 34 39 23 2 5 38 5 5 7 37 4 2 7 3 44 5 5 Total Itrs 20 24 24 28 29 38 24 29 33 38 19 37 30 33 34 39 23 22 30 39 28 32 34 40 35 23 36 52 54 31 32 Time (secs.) 0.05 0.13 0.18 0.26 0.32 0.47 0.49 0.70 0.80 1.49 0.44 4.01 1.26 1.66 1.58 2.28 0.62 0.84 1.01 3.83 1.39 1.35 5.52 5.87 1.75 1.82 7.78 3.64 31.70 2.44 3.43 Objective Value -4.6475315e+02 2.2549496e+05 -2.3313898e+06 -5.2202061e+01 -4.1573224e+02 -7.6589319e+04 1.8781248e+03 -1.4753433e+07 1.4122488e+03 1.5185099e+03 8.6666666e+00 -8.9664483e+05 -1.5862802e+02 1.8416759e+04 -1.8751929e+01 9.0429695e+02 3.3592486e+04 5.0500000e+01 1.7987145e+06 3.6660261e+04 1.7933245e+06 1.9200982e+06 1.7248068e+03 5.4901254e+04 1.4892361e+06 9.0500000e+02 1.4240000e+03 2.1851927e+06 5.5018494e+03 1.9090552e+06 1.4701879e+06 conjugate gradient algorithm was triggered when running this test problem. • If the test problems are categorized into three groups according to number of nonzero elements: ”small” (Afiro–Ship04s), ”medium” (Scfxm2–Ship12s) and ”large” (Scsd8– Ship12l), Algorithm I is, for each corresponding category, on the average 1.8, 4.0 and 6.8 times faster than MINOS. One may infer a growth in the relative speed of Algorithm I with respect to MINOS as problem sizes increase. AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 21 TABLE 8.1.5. Algorithm II test statistics (IBM 3081-K) — NETLIB test problems (minimum degree ordering heuristic in Algorithm II) Problema Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israel BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb Ship08l Ship12l Phase I Itrs 1 1 3 3 4 5 4 3 5 24 2 6 4 30 30 29 17 2 4 29 4 4 5 30 3 2 6 3 4 4 Total Itrs 15 18 19 20 21 33 19 21 23 24 16 29 24 30 30 29 17 18 22 29 21 21 25 30 23 18 27 35 23 23 Time (secs.) 0.08 0.27 0.43 0.66 0.69 1.47 1.22 1.71 2.00 2.79 1.10 8.22 3.46 4.94 3.54 5.07 1.51 2.13 2.85 9.28 4.71 3.30 13.55 14.25 4.17 4.22 19.59 14.33 8.12 10.07 Objective Value -4.6475315e+02 2.2549496e+05 -2.3313898e+06 -5.2202061e+01 -4.1573224e+02 -7.6589319e+04 1.8781248e+03 -1.4753433e+07 1.4122488e+03 1.5185099e+03 8.6666666e+00 -8.9664483e+05 -1.5862802e+02 1.8416759e+04 -1.8751929e+01 9.0429695e+02 3.3592486e+04 5.0500000e+01 1.7987145e+06 3.6660261e+04 1.7933245e+06 1.9200982e+06 1.7248068e+03 5.4901254e+04 1.4892361e+06 9.0500000e+02 1.4240000e+03 2.1851927e+06 1.9090552e+06 1.4701879e+06 a Problem 25FV47 was not solved with Algorithm II, as its current implementation does not incorporate a dense window data structure (Adler et al. 1989) necessary to solve this problem under a 4 Mbytes memory limit. • In Table 8.1.7, where problems with similar structure are grouped together, one can observe that MINOS’s disadvantage in number of iterations grows with problem size, while its advantage in time per iteration seems to level off for the larger problems. This is so, because MINOS refactors the basis after a fixed number of iterations. Refactorization requires work that is in the same order of one iteration of Algorithm I, dominating the work carried out in the intermediate iterations. Also, as the problem 22 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.1.6. Comparison of run times (IBM 3090) — NETLIB Test Problems (minimum local fill-in ordering heuristic) Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israel BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47 Ship08l Ship12l a IBM MINOS Time (secs.) 0.01 0.23 0.29 0.45 0.26 1.27 1.06 4.54 1.98 1.34 1.02 1.54 2.54 2.04 3.68 5.64 0.18 2.69 3.54 8.18 7.18 9.16 22.02 19.31 12.78 14.53 32.87 39.08 217.67 14.72 26.10 Alg-I Time (secs.) 0.05 0.13 0.18 0.26 0.32 0.47 0.49 0.70 0.80 1.49 0.44 4.01 1.26 1.66 1.58 2.28 0.62 0.84 1.01 3.83 1.39 1.35 5.52 5.87 1.75 1.82 7.78 3.64 31.70 2.44 3.43 Alg-IIa Time (secs.) 0.08 0.27 0.43 0.66 0.69 1.47 1.22 1.71 2.00 2.79 1.10 8.22 3.46 4.94 3.54 5.07 1.51 2.13 2.85 9.28 4.71 3.30 13.52 14.25 4.17 4.22 18.95 14.33 NA 8.12 10.07 MINOS Alg-I Time Ratio MINOS Alg-II Time Ratio 0.2 1.8 1.6 1.7 0.8 2.7 2.2 6.5 2.5 0.9 2.3 0.4 2.0 1.2 2.3 2.5 0.3 3.2 3.5 2.1 5.2 6.8 4.0 3.3 7.3 8.0 4.2 10.7 6.9 6.0 7.6 0.5 1.7 1.2 1.6 0.9 1.8 2.1 6.1 1.7 1.1 2.2 0.4 1.9 1.0 1.7 1.9 0.3 3.9 2.1 1.8 2.5 4.6 2.8 2.5 5.6 7.7 3.0 4.3 NA 4.0 5.5 3081-K CPU times and minimum degree ordering heuristic in Algorithm II. size increases, the overhead incurred by Algorithm I in the preprocessing is absorbed and therefore one should expect the ratio of work per iteration between the two algorithms to be at worst proportional to the inverse of the refactorization frequency. AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 23 TABLE 8.1.7. Iteration, time per iteration and total time ratios for NETLIB subgroups (IBM 3090) (minimum local fill-in ordering heuristic) Problem Rows Columns Scsd1 Scsd6 Scsd8 Scfxm1 Scfxm2 Scfxm3 ScTap1 ScTap2 ScTap3 Ship04s Ship08s Ship12s Ship04l Ship08l Ship12l 77 147 397 330 660 990 300 1090 1480 402 778 1151 402 778 1151 760 1350 2750 600 1200 1800 660 2500 3340 1506 2467 2869 2166 4363 5533 MINOS Alg-I Iter. Ratio 13.7 26.2 69.1 13.2 20.3 29.9 9.2 33.0 35.4 10.3 20.5 20.5 21.2 30.9 39.6 MINOS Alg-I Time/Iter. Ratio 0.169 0.122 0.115 0.093 0.105 0.110 0.268 0.121 0.119 0.340 0.330 0.356 0.244 0.195 0.192 MINOS Alg-I Time Ratio 2.32 3.20 7.98 1.23 2.14 3.29 2.47 3.99 4.22 3.50 6.79 7.30 5.17 6.03 7.61 Thus, relative to simplex based codes, Algorithm I becomes increasingly faster as the problems get larger. • The optimal objective values reported in Table 8.1.3 are accurate within 6 and 8 digits if compared to the values reported in (Gay 1985). The primal solution computed at termination was accurate in most cases with the exception of problems CzProb and 25FV47 (see Table 8.1.8). The primal-dual relation presented for problem Israel are based on an implementation using direct factorization to compute the search directions at each iteration. • All of the above considerations made for Algorithm I are valid for Algorithm II. Furthermore, Algorithm II requires on average 26% less iterations to converge to 8 digit accuracy than Algorithm I. However, this does not translate into any significant gain in solution time, since the extra work per iteration in the current implementation offsets the savings in the number of iterations. We believe that by refining the extra computation, further reduction in work per iteration in Algorithm II is possible, leading to a faster implementation. • The conjugate gradient algorithm, with column density parameter λ = 0.3, was important for solving test problem Israel. Factorization times for that problem were reduced significantly, when compared with an earlier version of Algorithm I, where the conjugate gradient routine was not implemented. Even though dropping a few dense columns did not significantly affect convergence of the dual solution, it did not generate a direction precise enough for computing an accurate primal solution at termination. This, however, can be accomplished by means of a more exact solution to the system defining the last direction. Actually, by using an exact factorization in the last 24 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.1.8. Primal dual relations for Algorithm I (IBM 3090) — NETLIB test problems Problem Afiro ADLittle Scagr7 Sc205 Share2b Share1b Scorpion Scagr25 ScTap1 BrandY Scsd1 Israel BandM Scfxm1 E226 Scrs8 Beaconfd Scsd6 Ship04s Scfxm2 Ship04l Ship08s ScTap2 Scfxm3 Ship12s Scsd8 ScTap3 CzProb 25FV47 Ship08l Ship12l Relative Duality Gap 3.56e-09 3.69e-09 3.16e-09 2.41e-09 6.83e-11 3.33e-09 2.84e-09 1.06e-08 2.25e-09 1.34e-09 3.31e-09 1.22e-09 1.06e-08 1.20e-07 4.10e-09 1.70e-09 1.44e-09 9.81e-10 8.34e-09 1.54e-08 2.15e-09 1.56e-09 1.06e-08 9.29e-09 8.16e-09 4.26e-09 1.73e-09 1.60e-09 1.22e-08 4.11e-10 1.27e-09 Max. Normal. Primal Infeas. 1.58e-15 1.05e-11 2.31e-10 2.23e-13 2.09e-12 3.26e-09 5.27e-12 1.32e-10 1.41e-13 4.00e-06 4.90e-13 5.46e-12 2.12e-09 7.29e-07 3.36e-07 7.28e-06 1.42e-06 4.30e-10 2.41e-09 5.99e-10 2.25e-11 1.09e-11 9.57e-15 4.11e-10 1.78e-09 1.35e-14 1.08e-14 1.11e-05 5.43e-03 1.02e-09 9.38e-10 Min. Normal. Primal Entry -5.76e-24 -9.67e-22 -4.02e-20 -2.61e-20 -1.06e-16 -1.73e-21 -6.07e-13 -9.73e-23 -2.96e-21 -6.91e-04 -1.05e-14 -5.38e-24 -1.28e-16 -8.78e-19 -5.94e-05 -1.67e-09 -6.48e-06 -1.61e-16 -4.95e-14 -1.78e-19 -1.75e-16 -4.24e-18 -6.09e-20 -5.37e-20 -1.74e-15 -4.67e-21 -1.01e-20 -1.02e-19 -3.17e-06 -3.05e-18 -2.19e-17 Max. Normal. Complem. Viol. 7.99e-12 3.06e-12 6.27e-13 2.73e-13 1.84e-12 8.31e-15 1.14e-12 4.50e-13 9.40e-14 5.43e-14 1.17e-11 1.11e-15 1.76e-13 5.03e-13 1.51e-13 3.55e-15 9.23e-12 7.15e-13 1.45e-11 2.19e-13 3.01e-12 1.40e-12 2.00e-13 8.30e-14 2.96e-12 6.15e-13 4.23e-14 1.19e-14 7.54e-15 4.29e-13 2.83e-13 iteration the accuracy of the primal solution for problem Israel was similar to those reported for the other problems. • Comparing Table 8.1.3 and Table 8.1.4 we observe a clear advantage in using the minimum local fill-in ordering heuristic. While the ordering algorithm for this heuristic usually involves longer processing times, the reduction in fill-in results in a faster Gaussian elimination procedure. Since the same ordering is used in every iteration AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 25 TABLE 8.2.1. Multi-commodity network flow test problems statistics Problem MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 Nodes 300 400 500 600 400 500 600 300 400 500 600 Commodities 1 1 1 1 3 3 3 5 5 5 5 LP Rows 1406 1888 2360 2853 2008 2489 2981 2081 2793 3477 4135 LP Columns 2606 2488 4360 5253 4027 4949 5882 4436 5948 7637 8800 Nonz(A) 5212 6976 8720 10506 9670 11876 14126 11196 15068 19182 22140 of the linear programming algorithm, the total savings more than offsets the extra processing in the ordering procedure. The sum of the solution times (including reordering) for all test problems was reduced from 119.10 seconds to 89.11 seconds with minimal local fill-in, corresponding to a 25% reduction. Solution times for Algorithm I with minimal local fill-in were faster on 21 out the 31 test problems and on 11 of the 12 problems having more than 5000 nonzero elements (see deCarvalho (1987)). 8.2. Multi-Commodity Network Flow Problems. In this section, we report on 11 multicommodity network flow problems generated with MNETGN (Ali & Kennington 1977), a random multi-commodity network flow problem generator. MNETGN generates a random network structure based on the number of arcs and nodes supplied by the user. Additional user specifications include the number of commodities and the ranges of arc costs and capacities. Table 8.2.1 displays the basic data for this set of test problems. We generated multi-commodity network flow problem by varying the number of nodes and commodities and keeping the same relative number of arcs. The test problems were solved with Algorithm I, MINOS and MCNF85 (Kennington 1979), a special purpose simplex code for multi-commodity network flow problems. Algorithm I used the minimum degree ordering heuristic and the same parameter selection described in Section 8.1. All runs were carried out on an IBM 3090 Tables 8.2.2–8.2.4 and Figure 8.2.1 present the statistics for these runs. We list below a few observations on the test results reported on this section. • Multicommodity networks provide an important class of test problems on which a general purpose simplex method performs poorly (see Table 8.2.3 and Figure 8.2.1). The behavior of our implementation of Algorithm I confirms our earlier observation that the relative speed-up in relation to the simplex method grows with size. A similar leveling off of the time per iteration ratios between the two simplex based codes and Algorithm I is observed (see Tables 8.2.5 and 8.2.6). The only exception was for the MCNF85 to Algorithm I ratio on the 5-commodity problems. The trend in the data suggests that leveling off should only occur for larger networks, i.e. networks with more than 600 nodes. As before, the number of iterations grows slowly, never exceeding 36. 26 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.2.2. Algorithm I test statistics (IBM 3090) — Multicommodity Networks (minimum degree ordering heuristic in Algorithm I) Problem MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 Phase I iter. 3 2 2 2 2 2 2 2 2 2 2 Total iter. 33 28 33 31 27 30 28 30 30 36 35 Time (secs.) 6.51 11.79 23.07 30.96 24.27 47.84 58.57 55.81 104.34 204.71 309.50 TABLE 8.2.3. MINOS 4.0 test statistics (IBM 3090) — Multicommodity Networks) Problem MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 Phase I Itrs 362 457 620 817 2498 3366 4855 2930 4611 5919 7693 Total Itrs 940 1351 1807 2358 5926 8115 12286 6691 12193 15238 21915 Time (secs.) 12.73 24.11 40.65 63.88 196.38 345.69 677.31 249.71 666.58 1045.67 1885.34 • It is interesting to note that despite being a general purpose implementation, Algorithm I performed comparably to MCNF85. A tailored implementation of Algorithm I (e.g. with integer operations and specialized Gaussian elimination) could improve on current results. Furthermore, the trend in the data (Table 8.2.6) suggests that for larger problems Algorithm I could outperform MCNF85. • The objective values for the solutions obtained with our implementation of Algorithm I achieved 8 digits accuracy when compared to the values obtained by MINOS. Also a feasible primal solution was recovered at the end of the algorithm with the same degree of accuracy. 8.3. Timber Harvest Scheduling Problems. In this section, we report on 11 timber harvest scheduling problems generated with FORPLAN (Johnson 1986). Based on data collected for a United States national forest, a series of timber harvest scheduling models of AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 27 TABLE 8.2.4. MCNF85 test statistics (IBM 3090) — Multicommodity networks Problem MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 Total iter. 931 1087 1892 3082 4983 5158 16983 4132 10558 8862 16624 Time (sec.) 7.42 11.39 24.87 53.45 45.23 55.60 243.35 34.64 111.30 121.66 260.44 TABLE 8.2.5. MINOS /Algorithm I performance ratios — Multicommodity Networks (minimum local fill-in ordering heuristic) Problem Nodes Commodities MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 300 400 500 600 400 500 600 300 400 500 600 1 1 1 1 3 3 3 5 5 5 5 MINOS /Alg-I iter. ratio 28.5 48.3 54.8 76.1 219.5 270.5 438.8 223.0 406.4 423.3 626.1 MINOS /Alg-I time/iter ratio 0.069 0.042 0.032 0.027 0.037 0.027 0.026 0.020 0.016 0.012 0.010 MINOS /Alg-I time ratio 1.96 2.04 1.76 2.06 8.09 7.23 11.56 4.47 6.39 5.11 6.09 increasing sizes were generated. Under the framework of FORPLAN a forest is divided into management units called analysis areas, each comprising a collection of acres from across the forest, sharing similar silvicultural and economic characteristics. Our test problems were created by successively increasing the number of analysis areas, resulting in models that represent subsets of the original forest. FORPLAN is widely used throughout the United States national forest system, yielding very large linear programs that pose a considerable challenge to linear programming codes. Table 8.3.1 presents the basic statistics of this family of test problems. From Table 8.3.1, we observe that the test problems have significantly more variables than constraints. This characteristic is accentuated with problem size. In this situation, most 28 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.2.6. MCNF85 /Algorithm I performance ratios — Multicommodity Networks (minimum local fill-in ordering heuristic) Problem Nodes Commodities MUL031 MUL041 MUL051 MUL061 MUL043 MUL053 MUL063 MUL035 MUL045 MUL055 MUL065 300 400 500 600 400 500 600 300 400 500 600 1 1 1 1 3 3 3 5 5 5 5 MCNF85 /Alg-I iter. ratio 28.2 38.8 57.3 99.4 184.6 171.9 606.5 137.7 351.9 246.2 475.0 MCNF85 /Alg-I time/iter ratio 0.0404 0.0249 0.0188 0.0174 0.0101 0.0068 0.0069 0.0045 0.0030 0.0024 0.0018 MCNF85 /Alg-I time ratio 1.14 0.97 1.08 1.73 1.86 1.16 4.15 0.62 1.07 0.59 0.84 b 1800 1600 b MINOS Algorithm I r MCNF85 × 1400 1200 b CPU Time 1000 (seconds) 800 b b 600 400 200 b b b × r r × r × b r r × ×br ×r ×r × ×r ×rb 0 ×rb 6000 8000 10000 12000 14000 16000 18000 20000 22000 Number of Nonzero elements F IGURE 8.2.1. MINOS 4.0, MCNF85 and Algorithm I (CPU times, IBM 3090) — Multicommodity Flows (minimum degree ordering heuristic in Algorithm I) advanced simplex implementations provide a choice of pricing strategies. In particular, MINOS allows for partial pricing , a scheme in which the columns of the coefficient matrix are partitioned into equal segments. In each pricing operation, the search for the incoming basic variable is limited to one segment. This reduces the work required for each operation, but, of course, has no predictable effect on the total number of iterations. For this class of AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 29 TABLE 8.3.1. Timber harvest scheduling test problems statistics Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 LP Rows 55 58 64 76 87 109 130 177 220 266 316 LP Columns 744 973 1502 2486 3357 4993 6667 9805 12570 16958 19991 Nonz(A) 6021 8168 12765 21617 30104 43774 59057 86573 109697 149553 176346 test problems, the manipulation of the PARTIAL PRICE parameter, which sets the number of segments in the partition, can impact the performance of MINOS enormously. Consequently, we compare MINOS with Algorithm I using the following pricing strategies: • Total Pricing Strategy: Each pricing operation examines the total set of variables. • Default Pricing Strategy: Sets the PARTIAL PRICE parameter to half of the ratio of number of columns to number of rows. • Improved Pricing Strategy: Sets the PARTIAL PRICE parameter to eight times the ratio of number of columns to number of rows. We arrived at this strategy after extensive testing with some of the test problems. Tables 8.3.2–8.3.8 display the results of running the test problems with Algorithm I and MINOS under the three pricing strategies. Figure 8.3.1 illustrates graphically the computational performance of the algorithms. All run were carried out on an IBM 3090 with same characteristics as in tests reported in Section 8.1. We use minimum degree ordering heuristic in Algorithm I. We close this section with a few observations on the results of the runs with the timber harvest scheduling problems. • The theoretical worst case analysis of Karmarkar’s algorithm (Karmarkar 1984) suggests that the number of iterations should grow linearly with the largest dimension of the coefficient matrix. This behavior is not is not confirmed for the variants of the algorithm described here. In this set of test problems, with number of columns ranging from 744 to 19991, the number of iterations until convergence grows sublinearly, varying from 38 to 71, consistent with our earlier observation that the number of iterations grows slowly with problem size. • This collection of test problems present an unusual behavior when solved with MINOS. Decreasing the number of columns considered in each pricing operation reduced not only the computation effort in each simplex operation, but also the total number of iterations. The number of pivots associated with solution paths followed by different versions of the simplex method may vary greatly. In this sense, interior point methods display a more robust behavior. This recommends caution in carrying 30 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.3.2. Algorithm I test statistics — Timber harvest scheduling problems (minimum degree ordering heuristic) Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Phase I iter. 4 4 4 5 5 5 5 6 5 6 6 Total iter. 38 40 41 49 49 49 56 52 43 67 71 Time (sec.) 0.85 1.10 1.64 2.95 4.27 6.47 9.60 14.99 24.81 34.11 43.80 TABLE 8.3.3. MINOS 4.0 test statistics (total pricing) — Timber harvest scheduling problems Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial Pricing 1 1 1 1 1 1 1 1 1 1 1 Phase I iter. 449 579 688 1318 1521 2350 2867 4570 5375 7322 3488 Total iter. 534 736 878 1574 1827 2768 3446 5407 7300 10156 12009 Time (sec.) 3.52 5.63 9.87 28.81 43.18 93.77 154.78 351.82 600.75 1112.68 1582.80 out computational experiments involving the simplex method, as varying a single parameter among different runs of MINOS the execution time for the two largest test problems is reduced by a factor of over 100. • Our implementation of Algorithm I is faster than the most straightforward version of MINOS by a wide margin (see Table 8.3.3 and Figure 8.3.1) and is up to 2.8 times faster for runs using the default pricing strategy . However, MINOS runs up to 3 times faster when using the improved pricing strategy in a pricing strategy obtained after extensive experimentation. • The objective values for the solutions obtained with our implementation of Algorithm I achieved 8 digits accuracy when compared to the values obtained by MINOS. AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 31 TABLE 8.3.4. MINOS (Total Pricing) Algorithm I performance ratios — Timber harvest scheduling problems (minimum degree ordering heuristic in Algorithm I) Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial Pricing 1 1 1 1 1 1 1 1 1 1 1 MINOS Alg-I iter. ratio 14.05 18.40 21.41 32.12 37.29 56.49 61.54 103.98 169.77 151.58 169.14 MINOS Alg-I time/iter. ratio 0.2947 0.2782 0.2810 0.3040 0.2712 0.2566 0.2620 0.2257 0.1426 0.2152 0.2137 MINOS Alg-I time ratio 4.14 5.12 6.02 9.77 10.11 14.49 16.12 23.47 24.21 32.62 36.14 TABLE 8.3.5. MINOS 4.0 test statistics (default pricing) — Timber harvest scheduling problems Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial pricing 1 1 11 16 19 22 25 27 28 31 31 Phase I iter. 449 579 724 973 1460 1860 2707 3244 4536 6934 3421 Total iter. 534 736 902 1259 1832 2352 3369 4370 6035 9732 11364 Time (sec.) 3.52 5.63 2.74 4.48 7.39 11.37 17.94 29.74 49.13 93.89 123.62 Also a feasible primal solution was recovered at the end of the algorithm with the same degree of accuracy. 8.4. Interpretation of Computational Experiments. We make the following interpretation of the results of the computational experiments. • Our implementation of Algorithm I attained 7 and 8 digit accuracy in the objective function for most test problems without numerical difficulties. The data structures and programming techniques that led to these results are described in Adler et al. (1989) 32 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA TABLE 8.3.6. MINOS (Default Pricing)/Algorithm I performance ratios — Timber harvest scheduling problems (minimum degree ordering heuristic in Algorithm I) Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial pricing 1 1 11 16 19 22 25 27 28 31 31 MINOS Alg-I iter. ratio 14.05 18.40 22.00 25.69 37.39 48.00 60.16 84.04 140.35 145.25 160.06 MINOS Alg-I time/iter. Ratio 0.2947 0.2782 0.0759 0.0591 0.0463 0.0366 0.0311 0.0236 0.0141 0.0190 0.0176 MINOS Alg-I time ratio 4.14 5.12 1.67 1.52 1.73 1.76 1.87 1.98 1.98 2.75 2.82 TABLE 8.3.7. MINOS 4.0 test statistics (improved pricing) — Timber harvest scheduling problems Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial pricing 104 128 160 184 360 480 480 440 440 480 480 Phase I iter. 48 137 79 113 77 105 174 455 461 387 729 Total iter. 1738 812 478 333 379 542 558 924 1049 843 1329 Time (sec.) 2.54 1.52 1.04 1.13 1.89 2.96 3.96 6.61 8.41 9.30 13.26 • The selection of an initial interior solution as described in Section 6 plays a significant role in the fast convergence observed in all test problems. In this procedure, we attempt to start the algorithm far from faces of the feasible polyhedral set, avoiding the difficulties described by Megiddo & Shub (1989). • Recovering a dual solution (primal solution in the format of the test problems) is an intricate proposition. Convergence of the dual solution estimates in 5.6 is guaranteed under nondegeneracy only, and examples can be built were it converges to an infeasible solution. However, our test results display good behavior for the dual solutions. AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 33 TABLE 8.3.8. MINOS (improved pricing)/Algorithm I performance ratios — Timber harvest scheduling problems (minimum degree ordering heuristic in Algorithm I) Problem FPK010 FPK040 FPK050 FPK080 FPK100 FPK150 FPK200 FPK300 FPK400 FPK500 FPK600 Partial pricing 104 128 160 184 360 480 480 440 440 480 480 MINOS Alg-I iter. ratio 45.74 20.30 11.66 6.80 7.73 11.06 9.96 17.77 24.40 12.58 18.72 MINOS Alg-I time/iter. ratio 0.0653 0.0681 0.0544 0.0564 0.0572 0.0414 0.0414 0.0248 0.0139 0.0217 0.0162 MINOS Alg-I time ratio 2.99 1.38 0.63 0.38 0.44 0.46 0.41 0.44 0.34 0.27 0.30 1600 b Algorithm I r MINOS(Total pricing) b MINOS(Default pricing) × 1200 MINOS(Improved pricing) ⋆ 1000 1400 b CPU Time (seconds) 800 b 600 400 b 200 b × ×r r × r r r × × ⋆ ⋆ ⋆ × ⋆ ×× ⋆ ⋆⋆ ⋆ × ⋆ ⋆ ⋆ × 0 × 0 20000 40000 60000 80000 100000 120000 140000 160000 180000 Number of Nonzero elements bb rr rb b r b r b r F IGURE 8.3.1. MINOS 4.0 and Algorithm I (CPU times, IBM 3090) — Timber harvest problems (minimum degree ordering heuristic in Algorithm I) • The stopping criterion described in Section 6 resulted in the correct solution within the desired accuracy. An alternative criterion that checks dual complementarity properties could have been used without loss in efficiency. 34 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA • Using δ = 0.1 in updating the AT D−2 v A matrix, Algorithm I obtained savings of over 10% in execution times, when compared to full updating, without degradation of the method’s numerical stability. • Algorithm I is sensitive to the density of the AT D−2 v A matrix. Hence, besides being sparse, the AT matrix should have a small number of dense columns. Maintaining this matrix sparse is essential to the algorithm. 9. C ONCLUSION AND E XTENSIONS The computational results presented in this paper illustrate the potential of interior point methods for linear programming. While the implementation described here was developed in a short period of time, it still outperforms MINOS 4.0 on the majority of problems tested. For test problems from the NETLIB collection, solution time speed-ups 6 or higher are observed on most of mid-sized problems, and increase with problem size. Seven and eight digit accuracy in the objective function value was achieved on all test problems, except for problems ScTap1, CzProb and 25FV47 . An optimal complementary primal-dual pair was obtained for all test problems with exception of CzProb and 25FV47 . Numerical difficulties with the LU towards the end of the algorithm account for these inaccuracies. Among the test problems is a set of multi-commodity network flow problems. The interior algorithm was clearly superior to MINOS and was also competitive with a specialized simplex algorithm, MCNF85. For the set of timber harvest scheduling problems, linear programs with significantly more variables than constraints, the number of iterations required by Algorithm I grows slowly with the number of variables. However, specialized pricing strategies make the simplex method a better tool for solving this type of linear programming problems. The implementations of interior point algorithms described in this paper are still preliminary. Since many improvements are still possible this approach seems promising as a general purpose solver for large real-world linear programming problem s. Of the several extensions planned to our implementation, some are immediately required: • Implicit treatment of upper and lower bounds on the variables. • Feasibility adjustment of the tentative primal solution computed at the end of the algorithm. Other extensions planned are: • Preprocessor for increasing sparsity of input matrices. • Higher order approximations to the solution of system differential equation 3.9–3.11. • Optimal basis identification for early termination. • Bi-directional search for determining the step direction. • Develop primal and primal-dual implementations. • Include steps in potential function gradient direction. • Implementations for special LP structures, e.g. network flows and GUB. • Implementation for parallel computer architectures. A CKNOWLEDGMENT The authors acknowledge the help volunteered by several individuals in coding the numerous modules that make up the backbone of this implementation, which started as a joint class project for I. Adler’s Fall 1985 Large-Scale Systems graduate seminar at Berkeley. In particular, W. Bein, D. Chou, C.C. Chyu, S. Cosares, M. deCarvalho, J. Doucet, J. Harrison, T. Kaplan, R.C. Monteiro, D. Popken, P.C. Samaratunga, A. Svoboda and Z. Wei were involved in the programming effort. We also wish to thank M. Khellaf and A. Sanstad for AN IMPLEMENTATION OF KARMARKAR’S ALGORITHM 35 valuable discussions related to the implementation of Algorithm II. Furthermore, the authors are indebted to one of the anonymous referees for many constructive suggestions. This research was partially funded by the Brazilian Council for the Development of Science and Technology CNPq and Management Sciences Staff, US Forest Service USDA. R EFERENCES Adler, I., Karmarkar, N. K., Resende, M. G. C. & Veiga, G. (1989), ‘Data structures and programming techniques for the implementation of Karmarkar’s algorithm’, ORSA Journal on Computing 1, 84–106. Adler, I. & Monteiro, R. C. (1991), ‘Limiting behaviour of the affine-scaling continuous trajectories for linear programming problems’, Mathematical Programming 50, 29–51. Ali, A. I. & Kennington, J. L. (1977), MNETGN program documentation, Technical Report IEOR 77003, Department of Industrial Engineering and Operations Research, Southern Methodist University, Dallas, TX. Aronson, J., Barr, R., Helgason, R., Kennington, J., Loh, A. & Zaki, H. (1985), The projective transformation algorithm by Karmarkar: Computational experiment with assignment problems, Technical Report 85-OR-3, Department of Operations Research, Southern Methodist University, Dallas, TX. Barnes, E. R. (1986), ‘Variation on Karmarkar’s algorithm for solving linear programming problems’, Mathematical Programming 36, 174–182. Bayer, D. A. & Lagarias, J. C. (1989), ‘The nonlinear geometry of linear programming’, Transactions of the AMS 314(2), 499–526. Chandru, V. & Kochar, B. (1985), A class of algorithms for linear programming, Research Memorandum 85–14, Dept. of Industrial Engineering, Purdue University, West Lafayette, IN 47907, USA. Revised 1986. deCarvalho, M. L. (1987), On the work needed to factor a symmetric positive definite matrix, Technical Report ORC 87-14, Operations Research Center, University of California, Berkeley, CA. Dikin, I. I. (1967), ‘Iterative solution of problems of linear and quadratic programming’, Soviet Mathematics Doklady 8, 674–675. Dongarra, J. J. & Grosse, E. (1987), ‘Distribution of mathematical software via electronic mail’, Communications of the ACM 30, 403–407. Duff, I. S., Erisman, A. M. & Reid, J. K. (1986), Direct Methods for Sparse Matrices, Claredon Press, Oxford. Eisenstat, S. C., Gurshy, M. C., Schultz, M. H. & Sherman, A. H. (1982), ‘The Yale sparse matrix package, I. the symmetric codes’, International Journal of Numerical Methods in Engineering 18, 1145–1151. Gay, D. M. (1985), ‘Electronic mail distribution of linear programming test problems’, COAL Newsletter 13, 10– 12. Gill, P. E., Murray, W., Saunders, M. A., Tomlin, J. A. & Wright, M. H. (1986), ‘On projected Newton barrier methods for linear programming and an equivalence to Karmarkar’s projective method’, Mathematical Programming 36, 183–209. Gonzaga, C. (1991), ‘Interior point algorithms for linear programming problems with inequality constraints’, Mathematical Programming 52, 209–225. Ho, J. & Loute, E. (1981), ‘Set of staircase linear programming test problems’, Mathematical Programming 20, 245–250. Hooker, J. H. (1986), ‘Karmarkar’s linear programming algorithm’, Interfaces 16, 75–90. Johnson, K. N. (1986), FORPLAN version 1: An overview, Technical report, Land Management Planning Section, USDA, Forest Service, Fort Collins, CO. Karmarkar, N. (1984), ‘New polynomial-time algorithm for linear programming’, Combinatorica 4, 373–395. Karmarkar, N. K., Lagarias, J. C., Slutsman, L. & Wang, P. (1989), ‘Power series variants of Karmarkar type algorithms’, AT&T Technical Journal 68, 20–36. Kennington, J. (1979), A primal partitioning code for solving multicommodity flow problems (version 1), Technical Report 79008, Department of Industrial Engineering and Operations Research, Southern Methodist University, Dallas, TX. Lustig, I. J. (1985), Practical approach to Karmarkar’s algorithm, Technical Report SOL SOL 85-5, Systems Optimization Laboratory, Stanford University, Stanford, CA. Megiddo, N. & Shub, M. (1989), ‘Boundary behavior of interior point algorithms in linear programming’, Mathematics of Operations Research 14, 97–146. Murtagh, B. A. & Saunders, M. A. (1977), MINOS user’s guide, Technical Report 77-9R, Systems Optimization Laboratory, Stanford University, Stanford, CA. Murtagh, B. & Saunders, M. (1983), MINOS 5.0 user’s guide, Technical Report 83-20R, Systems Optimization Laboratory, Stanford University, Stanford, CA. 36 I. ADLER, N. KARMARKAR, M.G.C. RESENDE, AND G. VEIGA Rose, D. J. (1972), A graph-theoretical study of the numerical solution of sparse positive definite systems of linear equations, in R. C. Read, ed., ‘Graph Theory and Computing’, Academic Press, New York, NY, USA, pp. 183–217. Tarjan, R. E. (1983), Data Structures and Network Algorithms, SIAM, Philadelphia, PA, USA. Todd, M. J. & Burrell, B. P. (1986), ‘An extension of Karmarkar’s algorithm for linear programming using dual variables’, Algorithmica 1(4), 409–424. Tomlin, J. A. (1987), ‘An experimental approach to Karmarkar’s projective method for linear programming’, Mathematical Programming Study 31, 175–191. Tone, K. (1986), An implementation of a revised Karmarkar method, Interim report, Graduate School for Policy Science, Saitama University, Urawa, Saitama 338, Japan. Vanderbei, R. J., Meketon, M. J. & Freedman, B. (1986), ‘Modification of Karmarkar’s linear programming algorithm’, Algorithmica 1, 395–407. Yannakakis, M. (1981), ‘Computing the minimum fill-in is NP-complete’, SIAM J Alg Disc Meth 2, 77–79. (I. Adler, M. G. C. Resende and G. Veiga) D EPARTMENT OF I NDUSTRIAL E NGINEERING AND O PERATIONS R ESEARCH , U NIVERSITY OF C ALIFORNIA , B ERKELEY, CA 94720 E-mail address, I. Adler: ilan@ieor.berkeley.edu (N. Karmarkar) AT&T B ELL L ABORATORIES , M URRAY H ILL , NJ 07974 USA E-mail address: Narendra.Karmarkar@att.com Current address, M.G.C. Resende: AT&T Bell Laboratories, Murray Hill, NJ 07974, USA E-mail address: Mauricio.Resende@att.com Current address, G. Veiga: AT&T Bell Laboratories, Holmdel, NJ 07733, USA E-mail address: Geraldo.Veiga@att.com