0 Preprint (1999) ?{?
1
Adaptive Use of Iterative Methods in Predictor-Corrector Interior Point Methods for Linear Programming
Weichung Wang a
a
Department of Mathematics Education, National Tainan Teachers College, Tainan 700, Taiwan E-mail: [e-mail]
Dianne P. O'Leary b;
b
Dedicated to Richard Varga, for making the use of iterative methods a science, without diminishing the artistry.
In this work we devise e cient algorithms for nding the search directions for interior point methods applied to linear programming problems. There are two innovations. The rst is the use of updating of preconditioners computed for previous barrier parameters. The second is an adaptive automated procedure for determining whether to use a direct or iterative solver, whether to reinitialize or update the preconditioner, and how many updates to apply. These decisions are based on predictions of the cost of using the di erent solvers to determine the next search direction, given costs in determining earlier directions. We summarize earlier results using a modi ed version of the OB1-R code of Lustig, Marsten, and Shanno, and we present results from a predictor-corrector code PCx modi ed to use adaptive iteration. If a direct method is appropriate for the problem, then our procedure chooses it, but when an iterative procedure is helpful, substantial gains in e ciency can be obtained. Keywords: Interior point methods, linear programming, iterative methods for linear systems, adaptive algorithms, self-timing algorithms AMS Subject classi cation: Primary 65K05, 65F10, 90C05
Department of Computer Science and Institute for Advanced Computer Studies, University of Maryland, College Park, MD 20742, U.S.A. E-mail: [e-mail]
This work was supported in part by the National Science Foundation under Grants CCR 95-03126 and CCR 97-32022.
2
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
1. Introduction
minimize c x subject to Ax = b; (1) x 0; where c; x are real n-vectors, b is a real m-vector, and A is a real m n matrix of rank m, with m n. These methods typically solve a sequence of logarithmic
T
Interior point algorithms are now widely used to solve linear programming problems
barrier subproblems with the barrier parameter decreasing to zero. Newton's method is applied to solve the rst order optimality conditions corresponding to each of the logarithmic barrier subproblems. The bulk of the work in such algorithms is the determination of a search direction for each step. Gonzaga 19] and Wright 38] surveyed interior point methods, and many computational issues are addressed by Lustig, Marsten, and Shanno 26] and Andersen, Gondzio, Meszaros, and Xu 1]. Therefore, in this section we focus only on the linear systems arising in interior point methods. For de niteness, we consider the primal-dual formulation of interior point methods, but the linear algebra of primal formulations and dual formulations is similar. The search direction is usually determined by solving either the reduced KKT (Karush-Kuhn-Tucker) system, x = r + Ze X 1 e ; X 1Z A (2) y r A0 or the normal equations, formed by eliminating x from this system. De ning r = b Ax, r = c A y z, and D2 = Z 1 X, we obtain (AD2A ) y = AD2(r + Ze x 1 e) + r : (3) Here z is the vector of dual slack variables, is the barrier parameter, and X and Z are diagonal matrices containing x and z (respectively) on their main diagonals. Once y is determined from the normal equations, x may be easily computed from (X 1Z) x + A y = r + Ze X 1e: (4) Comparing the normal equations (3) and the KKT system (2), we observe that the matrix for the normal equations is positive de nite and symmetric, has smaller size (m m), and may be more dense. In contrast, the KKT matrix is symmetric inde nite and usually more sparse. One nice feature of these problems is that only D and the right hand side of the system change from step to step. Thus, the sparsity structure of the problem remains the same, in contrast to the linear systems arising in the simplex algorithm which di er by exchanges of columns of A. Some interior point algorithms
T d p p d T T d p T d
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
3
(e.g., OB1-R) solve one linear system with each matrix, while others (e.g., PCx) solve multiple systems. The roots of interior point algorithms date back to the algorithms of Fiacco and McCormick 10], but ever since interior point algorithms rst gained prominence in 1984 20], researchers have given attention to speeding up the step time through e cient solution of the linear system. Direct methods that rely on sparse matrix factorizations have been the most popular approaches (e.g., 25], 33]), although iterative methods for solving linear systems have also received a fair amount of attention. Karmarkar and Ramakrishnan reported computational results of Karmarkar's dual projection algorithm using a preconditioned conjugate gradient solver 21]. An incomplete Cholesky factorization of the matrix AD2 A was computed for one interior point step and then used as a preconditioner over several subsequent steps. In their experiments, Cholesky factorization was performed on average every 2 to 3 steps. Mehrotra used preconditioned conjugate gradients to solve the normal equations in a dual a ne scaling interior point algorithm 27]. He addressed issues such as the stopping criterion and the stability of the implementation. At each interior point step, an incomplete Cholesky factor was computed and used as the preconditioner. Carpenter and Shanno used a diagonal preconditioner for a conjugate gradient solver for the normal equations in an interior point method for quadratic and linear programs 3]. They also considered recomputing the preconditioner every other step. Portugal, Resende, Veiga, and Judice introduced a truncated primal-infeasible dual-feasible interior point method, focusing on network ow problems 32]. The preconditioned conjugate gradient algorithm was used to solve the normal equations. They initially used the diagonal of the matrix AD2 A as a preconditioner and replaced it by spanning tree preconditioners in later steps. Mehrotra and Wang 28] used an incomplete Cholesky factor of AD2 A as a preconditioner for conjugate gradients in a dual interior point method for network ow problems. Gill, Murray, Saunders, Tomlin, and Wright established the equivalence between Karmarkar's projected method and their projected Newton barrier method 15]. They used LSQR 31], preconditioned by an approximation to AD2A , to nd the search directions. Goldfarb and Mehrotra developed a relaxed version of Karmarkar's method that allows inexact projection 17]. They applied CGLS 31] to determine the search direction. Nash and Sofer investigated the choice of a preconditioner in the positive de nite system Z GZ where Z is rectangular and G is general symmetric 29]. Chin and Vannelli 5] solved a reduced KKT system using the preconditioned conjugate gradient algorithm and Bi-CGSTAB with incomplete factorization. In a di erent paper 4] they used an incomplete factorization as a preconditioner for the normal equations (3). Freund and Jarre 11] employed a symmetric variant of the quasi-minimal residual (QMR) method to solve the KKT systems. They suggested using inde nite SSOR preconditioners to accelerate the convergence.
T T T T T
4
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
Despite all of this work, the use of iterative methods has so far produced limited success. The obstacles to the use of these methods are considerable. Over the course of the interior point steps, the requirements on accuracy change greatly; approximate solutions can be allowed early in the steps but can cause the algorithm to fail later when the iterates are near the boundary. The matrix D changes quite rapidly and becomes highly ill-conditioned in the nal steps. For these reasons, it is di cult to nd a preconditioning strategy that produces good performance of iterative methods over the entire course of the interior point computation. In this paper we develop an adaptive algorithm that changes strategy over the course of the interior point algorithm. It determines dynamically whether the preconditioner should be held constant, updated, or recomputed, and it switches to a direct method when it predicts that an iterative method will be too expensive. In our experiments, we use a preconditioned conjugate gradient iteration on the linear system involving the matrix ADA , but our ideas could be extended to iterations involving the KKT formulation as well. This idea of choosing among various numeric algorithms depending on the timing performance or timing prediction of algorithmic components for a particular problem on a particular machine architecture was summarized in a 1995 report by O'Leary and Wang 36] and elaborated by Wang in his 1996 thesis 35]. This idea has proved quite useful in other numeric algorithms, such as a 1997 algorithm of Frigo and Johnson for computing Fourier transforms 12] and a 1998 proposal by Whaley and Dongarra for linear algebra computations 37]. In the next section, we discuss the characteristics of direct and iterative methods and present our preconditioner. Section 3 focuses on our algorithm for the adaptive choice of direct vs. iterative methods and the adaptive choice of a preconditioner. Numerical results are presented in x 4. Final comments are made in x 5.
T
2. The linear system solvers
The most expensive part of an interior point algorithm is determining the search direction by solving one or more linear systems. Either direct or iterative methods may be used for these systems. In this section, we focus on the solution of the normal equations (3). This discussion sets the goals to be accomplished in designing an e cient algorithm. We will assume that the columns of A have been permuted using standard techniques in order to improve sparsity in the Cholesky factor of AD2 A (e.g., 9], 24]).
T
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
5
2.1. Direct solvers: Cholesky factorization
Most existing linear programming interior point methods solve the normal equations by direct methods. The implementations OB1-R of Lustig, Marsten, and Shanno 25] and PCx of Czyzyk, Mehrotra, and Wright 7] are representative of these methods, and the iterative methods will be compared with these implementations. To solve equation (3), the OB1-R implementation computes a sparse Cholesky factorization of the matrix K = AD2 A as LPL , where L is a unit lower triangular matrix and P is a diagonal matrix. Forward and backward substitution are then applied to compute the search direction y . The OB1-R algorithm then checks whether A x is close enough to the arti cial variables (b Ax ). If not, iterative re nement using the factored matrix LPL is employed repeatedly until the one-norm of the di erence is su ciently small. To deal with the dense columns in A, the OB1-R algorithm adopts the method suggested by Choi, Monma, and Shanno 6]. The PCx implementation uses a similar strategy for the solution of linear systems, using the Ng-Peyton 30] sparse matrix code, with modi cation by replacing small pivots by a very large number, and again dealing with dense columns separately. The algorithm also performs iterative re nement using the conjugate gradient algorithm with the factorization as a preconditioner. There are three main disadvantages to direct methods. First, the iterative re nement used in the OB1-R code may fail if the matrix K = AD2 A is very ill-conditioned because the factorization may not be accurate enough to produce an iteration matrix with spectral radius less than one. Such a situation can occur when the primal and dual variables are near to the optimal solution, since then the matrix D is quite ill-conditioned. The iteration can also be a ected by ill-conditioning in A. Another potential problem of direct methods is ll-in. Although the dense columns of A can be treated separately, the remaining Cholesky factor may still be rather dense. This might be caused by di culty in detecting \dense" columns or by the nature of the problem. For example, network problems solved by linear programming may lead to a Cholesky factor that is much more dense than AD2 A even though A has no dense columns. Lastly, the direct algorithms must form and factor the matrix K = AD2 A each time is changed. This procedure may be expensive in time, especially when the problem size is large. If m n, the resulting matrix K may be small and easy to factor, but forming it can still be costly.
T T k k k T T T T
2.2. Iterative solvers: preconditioned conjugate gradients
A variety of iterative methods can be used to solve the normal equations or the KKT system. For de niteness, we focus on the preconditioned conjugate
6
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
gradient method for solving equation (3). In this method, we compute a sequence of approximate solutions that converge to the true solution. The work during each iteration involves one product of K with a vector, one solution of a linear system involving the preconditioner, and several vector operations. More details about the method can be found in 18]. The storage requirement for the preconditioned conjugate gradient method is quite low, amounting to a few vectors of length m. Although a matrix-vector multiplication Kv = (AD2A )v is required at each iteration, we may compute Kv as (A(D2(A v))) and thus need only to store the nonzeros of A and the diagonal of D rather than the matrix AD2A , which can be quite dense. The preconditioner should also be chosen to conserve storage. Since accuracy requirements for the search direction in the beginning phase of the interior point algorithm are quite low, only a few conjugate gradient iterations are required. As the primal and dual variables approach the optimal solution, the convergence tolerance must be tightened and more iterations are needed. The crucial issue in the preconditioned conjugate gradient algorithm is to nd a preconditioner for each step of the interior point method. A good preconditioner may dramatically accelerate the convergence rate and gain great computational savings. We consider some strategies for choosing the preconditioners in the next subsection.
T T T
2.3. The preconditioner
Convergence of the conjugate gradient iteration will be rapid if the preconditioned matrix has either a small condition number or great clustering of eigenvalues 18, Chap. 10]. We discuss our strategy for preconditioning. The basic preconditioner is the Cholesky factorization of one of the matrices that has been generated in the course of the interior point method. PCx always uses the sparse piece of the current matrix, but this requires frequent factorizations. An alternative to computing a new Cholesky factorization on every interior point step is to reuse the preconditioner that was computed for one value of the barrier parameter in order to solve systems for several successive values of 3] 21]. This reduces the computational work in forming the factorization. An incomplete Cholesky factorization, originally proposed by Varga 34], could be used in place of the Cholesky if density of the matrix factors is too great, but we do not pursue that idea in our implementations. Rather than keeping the preconditioner xed when changes, though, we can update it by a small-rank change, since the normal equations matrix is a ^ continuous function of . Let D be the current diagonal matrix and D be the ^ one for which we have a factorization AD2 A = LPL . De ne D = D2 D2
T T
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
7
and let a be the i-th column of matrix A. Since
i T T T
X ^ AD2A = AD2A + A DA = LPL + d aa ;
n T ii i T i i=1 T
(5)
^ ^^ we may obtain an improved preconditioner LP L by applying a rank- update to LPL , where n. This update may be computed as in 2] and 8]. We choose large enough to include most of the large magnitude terms in the summation. ^ Then we have factored a matrix that di ers from AD2A by a matrix of rank n . This di erence matrix can be expressed as a matrix of small norm plus one of small rank, and we can hope for rapid convergence of the conjugate gradient iteration. We now turn our attention to criteria for deciding when to keep or update or reinitialize the current preconditioner and how many iterations to perform.
T T
3. The adaptive algorithm
Our interior point algorithm chooses the initial variables, the step lengths, the barrier parameter , and convergence criteria following standard strategies 7,25]. Each step requires the solution of one or more linear systems, and that is where the bulk of the computational work lies. The di erence between our algorithm and standard ones is that for each step of the interior point method (each distinct value of ), we choose an e cient linear equation solver adaptively. We need to specify when to use an iterative method, when to refactor the matrix, how many updates to use in the preconditioner, and how to terminate the iteration.
3.1. When to use the iterative method
In the rst step of the algorithm, the normal equations (3) are solved directly by factoring K = AD2 A = LPL . Starting from the second step, the algorithm uses preconditioned conjugate gradients. The preconditioner for each step is determined by factoring the current matrix K or by updating the current preconditioner. This \factor-update cycle" will be continued up to the \end-game," entered when the relative duality gap is small enough. In the end-game, the iterates are close to the optimal solution and accuracy requirements are high. The elements in matrix D vary signi cantly and make the matrix K = AD2 A very ill-conditioned. The Cholesky factorization of K may not generate a good preconditioner, even if stable methods such as 14] are used. For all of these reasons, a direct method is used to determine the nal search directions. We also switch to a direct method when OB1-R computes a Cholesky factorization with a zero on the diagonal. This contingency could be avoided by using a modi ed Cholesky factor; see, for example, 16, Chap. 4].
T T T
8
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
between refactorization drops to 3 or fewer. At that point (the middle stage), we begin to force a refactorization at least every 3 steps. This continues until the relative gap drops below a user-de ned tolerance (1:0 10 8 in the current implementation) at which time a refactorization is performed at least every other step (the late stage), and then the algorithm proceeds with the end-game as above. While the adaptive algorithm monitors the cost of the iterative method, it separates out problems that are not well suited to iterative methods. If twice in a row the updated preconditioner is ine cient in the step after the preconditioner is reinitialized, then the algorithm will use only the direct method from then on. An example of such a situation is illustrated in Figure 5. In summary, our algorithm uses direct methods for linear systems in the rst step, in the nal (end-game) steps, periodically in the middle and late stages of predictor-corrector methods, and at other times in which the iterative method is estimated to be more expensive than the direct method.
3.2. Deciding whether to refactor or to update the preconditioner
?,?], the \factor-update cycle" is modi ed when the number of di erent values
For interior point methods like PCx that use the predictor-corrector strategy
We make decisions regarding refactorization or update of the preconditioner based on the actual cost incurred in determining previous search directions, as measured in seconds by a system timing program: = the cost of factoring and solving the system directly; = the cost of each rank-one update; = the cost of each conjugate gradient iteration. (For simplicity, we neglect the fact that updates and downdates have slightly di erent costs.) We initialize each of these estimates to zero, but after the rst few steps of the interior point method, we have accurate estimates of each. In order to reduce the e ects of variability from the timer output, though, we suggest that these estimates continue to be updated over many steps. The decision to update the current preconditioner or refactor the matrix AD2 A to obtain a new preconditioner is based on the approximate cost of the preceding iteration, including the cost of any updates that were made to the preconditioner. This cost is prev cost = (updt cost updt nmbr) +(pcgi cost pcgi nmbr) +(overhead); where updt nmbr is the number of updates that were performed and pcgi nmbr is the number of pcg iterations. The overhead includes operations such as initializing the solution to zeros, computing the norm of the right-hand side, deciding on the number of rank-one updates, etc.
drct cost updt cost pcgi cost
T
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
T
9
If the cost of determining the previous search direction was high, we reinitialize the preconditioner by factoring the current matrix K = AD2 A . We take this action when the cost of the previous iteration exceeds 80% of the cost of direct solution:
prev cost
> :8
drct cost:
If the cost of the previous iteration was not that high, then we base our decision on a prediction of the cost of the current iteration, refactoring if the predicted cost is greater than the cost of the direct method. Our prediction method is simple and requires only a few arithmetic operations. We t a straight line to the number of iterations required to determine two preceding search directions. We choose the previous number, and the latest other one that gives a line with positive slope, and use this line to predict the number of iterations, predi nmbr, required to determine the current search direction. If the solver refactored on the previous step, or if we cannot obtain a positive slope with data since the last refactorization, then our predicted number of iterations is one more than the number taken last time, predi nmbr = pcgi nmbr+1. Given this predicted number of iterations, our predicted cost for computing the search direction, neglecting overhead, is
pred cost =
(updt cost
updt nmbr) + (pcgi cost
predi nmbr):
If this cost is less than drct cost, then the preconditioner is obtained by updating the previous one. Otherwise it is obtained by factoring K = AD2 A .
T
3.3. The adaptive updating strategy
We adopt the strategy discussed in x 2.3: we update the Cholesky factors using the updt nbmr= \largest" outer product matrices as determined by j d j. (We could have used j d jka k2 instead.) We change the number of Cholesky updates adaptively over the course of the algorithm in order to improve e ciency. The number is increased if the previous search direction took many iterations, and decreased if it took a very small number. Two parameters sml < lrg are initially set to 20 and 30 respectively. The parameter sml denotes a number of conjugate gradient iterations that takes time much less than drct cost, while lrg denotes a number that requires a more substantial fraction of drct cost. After timing data is available, we set
ii ii i
lrg
= 0:15
; pcgi cost
drct cost
sml
= 0:12
drct cost pcgi cost
:
10
pcgi slope
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
To decide the number of rank-one updates, updt nmbr, to be performed, let be the slope of the line connecting last two pcgi nmbrs. 8 < increased, if lrg pcgi nmbr and pcgi slope > 0, The updt nmbr is : decreased, if pcgi nmbr sml and pcgi slope < 0, unchanged, otherwise: (to increase)
updt nmbr = updt nmbr
Increases or decreases in updt nmbr are proportional to the pcgi slope: max(1:2;
pcgi slope
8:0
);
(to decrease) updt nmbr = (updt nmbr 0:9) + 1 : Note that the sparsity of the Cholesky factors remains the same, no matter how many updates are used.
3.4. Terminating the PCG iteration
After computing the preconditioner, we solve the normal equations using the preconditioned conjugate gradient method. We start from an initial guess of zero, iterate and until the computed residual norm is less than a parameter " times the norm of the right-hand side. We choose the parameter " adaptively: for OB1-R, (8 10 ; if relgap > 10 2; 1 " = 10 8 (relgap) 2 ; otherwise, where relgap is the relative duality gap for the previous value of . This is similar to the stopping criterion in 28]. For PCx, we use 8 5:0 10 3 ; for the beginning stage, < for the middle stage, " = : min(relgap 103, 1:0 10 3 ), min(relgap 104, 1:0 10 4 ), for the late stage,
pcg pcg pcg pcg
If the preconditioned conjugate gradient iteration number exceeds the maximum number of iterations allowed, then the current preconditioner is abandoned and a new preconditioner is determined by Cholesky factorization. If this happens twice, the iterative method is not suitable and we switch to a direct method. Unfortunately, the preconditioned conjugate gradient iteration might be stopped just before convergence, thereby making the refactoring wasteful, but we consider such a safeguard bounding the number of iterations to be important. The maximum number of iterations is set to the number that produces a cost of 1:2 times the cost of a direct method:
max pcg itn
= 1:2
drct cost pcgi cost
:
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
11
Problem: pds10 180 169 160 Reinit. preconditioner Updt. preconditioner Max PCG itnallowed lrg/sml
140
Number of PCG iterations
120
100
80
60
40
21 16 0
2
20
40
60 Outer iteration
80
100
118
Figure 1. Number of PCG iterations for the adaptive algorithm
To summarize, our algorithm solves the normal equations directly to determine the rst search direction, uses a preconditioned conjugate gradient method starting from the second search direction, and switches back to the direct method for the nal search directions. The preconditioned conjugate gradient solver solves the normal equations by rst choosing and computing a preconditioner using an adaptive strategy to decide whether to refactor the matrix and the rank of the update performed. The algorithm automatically sets all parameters expected to in uence performance, based on actual time performance of the components of the algorithm.
4. Numerical results
4.1. OB1-R
We modi ed the code OB1-R to adaptively choose the linear system solver, and we performed numerical experiments comparing the results of this modi ed version of OB1-R to the standard OB1-R code, dated December 1989.
12
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
Problem: pds10 400 Use direct solver Estimated direct solver cost With reinitialized preconditioner With updated preconditioner
350
300
250
Time (seconds)
200
143 114 100
50
0
1
11
20
31
40
60 Outer iteration
80
100
119
Figure 2. Timing performance for the adaptive algorithm
Both OB1-R and the adaptive algorithm are coded in FORTRAN and use double precision arithmetic. Our experiments were performed on a SUN SPARCstation 20 with 64 megabytes of main memory, running SunOS Release 4.1.3. The FORTRAN optimization level was set to -O3. We report CPU time in seconds, omitting the time taken by the preprocessor HPREP since it is the same for both codes. Before comparing the two codes, we illustrate the behavior of the adaptive algorithm on a large problem, pds-10 (with arti cial variables) whose problem characteristics are given in Table 1. Figure 1 shows the number of iterations needed by the preconditioned conjugate gradient method for the values chosen by OB1-R. Conjugate gradients are used for 2 through 118, and then the algorithm chooses to switch to direct solution because it detects a zero on the diagonal of the preconditioner. The horizontal line at 169 marks the maximum number of conjugate gradient iterations allowed (i.e. max pcg itn). The two dashed lines at 21 and 16 indicate lrg and sml, respectively. The Cholesky factorization is recomputed 25 times, marked by circles in the gure. This is a savings of 92 factorizations compared to the OB1-R algorithm. In between refactorizations, the
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
13
number of conjugate gradient iterations generally grows, more quickly for later values of than for earlier ones. Figure 2 displays the time taken by each of these linear system solves. The dashed line is drct cost, the estimated direct solver cost based on its performance for the rst value of . The solid line marks 0:8 times drct cost. We highlight the following observations from the gures. The adaptive algorithm produces signi cant savings in the beginning stage, especially from the 11th to the 31st value of . The frequency of reinitializing the preconditioner grows as is decreased. The preconditioners obtained from refactoring the matrix AD2 A are unsuitable in the late stage. The adaptive algorithm succeeds in keeping the cost near or better than the direct cost on all iterations but three. On those, the predicted number of iterations is too low.
T
We now summarize computational results on various types of linear programming problems. If the total time for solution is small (i.e., 5 minutes or less), then the performance of the two algorithms is similar. On more costly problems, the adaptive method is faster: e.g., 9% faster on pilot87, 16% faster on dfl001, and 28% faster on maros-r7 from the NETLIB collection. More complete results can be found in 36].
4.2. PCx
We modi ed the code PCx to adaptively choose the linear system solver, and we performed numerical experiments comparing the results of this modi ed version of PCx to PCx version 1.1, dated November 1997. Both PCx and the adaptive algorithm are coded by FORTRAN and C language and use double precision arithmetic. Our experiments were performed on a HP C1100/9000 workstation with 128 megabytes of main memory, running HP-UX B.10.20 operating system. The FORTRAN and C optimization level were set to -O. We report CPU time in seconds, including the time taken by the preprocessor which is the same for both codes. We report computational results on various types of linear programming problems chosen from NETLIB, NETLIB's Kennington problems, and some network problems. We omit data for problems taking fewer than 10 seconds, since direct methods are quite suitable for these. Table 1 summarizes the problem characteristics. The numbers of rows and columns indicated in the table refer to the output from the PCx preprocessor and may be di erent from the data in 13]. The tabulated number of nonzero
14
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
Problem: pds20 Use direct solver Estimated direct solver cost With reinitialized preconditioner With updated preconditioner
250
206
Time (seconds)
150
100
50
0
12
10
20
30 33 Outer iteration
40
50
54
Figure 3. Timing performance for the adaptive algorithm on problem pds-20.
elements of the Cholesky factor L include the diagonal part of L. The density of L is computed as b 2 (Number of nonzeroes of L) m ; bm b m b where m is the number of rows of A after presolving. Minimum cost ow network problems may be solved using linear programming algorithms (although it is generally more e cient to use a network algorithm like 22]). We test our algorithm on this class of problems because the matrix AA and its resulting Cholesky factor tend to be much more dense than the original coe cient matrix A, even if there is no dense column in A. Forming and factoring the matrix AD2 A is thus quite expensive. We generated minimum cost ow network problems using NETGEN, developed by Klingman, Napier, and Stutz 23].
T T
Before comparing the two codes, we illustrate the typical behavior of the adaptive algorithm using three examples.
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
15
Problem NETLIB: d2q06c degen3 d 001 t2d greenbea maros-r7 pilot pilot87 stoch3 Kennington: cre-b cre-d ken-11 ken-13 ken-18 osa-07 osa-14 osa-30 osa-60 pds-06 pds-10 pds-20 Network: net0108 net0116 net0408 net0416 net0816 net0832 net0864 net1632
LP size Rows Columns 2132 1503 5984 25 1933 2152 1368 1971 15362 5336 4102 10085 22534 78862 1081 2300 4313 10243 91556 15648 32287 1000 1000 4000 4000 8000 8000 8000 16000 5728 2604 12143 10524 4164 7440 4543 6373 22228 36382 28601 16740 36561 128434 25030 54760 104337 243209 28472 48780 106180 8000 16000 8000 16000 16000 32000 64000 32000
Cholesky factor L Nonzeros Density 137349 120906 1638085 324 49055 534188 200812 425654 177936 248629 212094 102906 298417 1928863 28276 60795 115081 265909 589339 1687660 7089645 207560 280678 556366 1766394 2201390 7000874 12247346 8653616 0.060 0.106 0.091 0.997 0.026 0.230 0.214 0.219 0.001 0.017 0.025 0.002 0.001 0.002 0.047 0.023 0.012 0.005 0.014 0.014 0.014 0.414 0.560 0.069 0.221 0.069 0.219 0.383 0.068
Table 1 Statistics of the test problems.
Figure 3 shows the time taken by linear system solves on problem pds-20. The time for computing predictor and corrector are summed. The dashed line is drct cost, the estimated direct solver cost based on its performance for the rst value of . Preconditioned conjugate gradients are used for 2 through 54 . The algorithm switches to the late stage at 33 and then to the ending stage at 55 because it detects a relative duality gap that is smaller than the parameter ending tol. The Cholesky factorization is recomputed 15 times, marked by cir-
16
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
Problem: dfl001 90 Use direct solver Estimated direct solver cost With reinitialized preconditioner With updated preconditioner
80
70
60
Time (seconds)
50
39.5
30
20
10
0
12
5
10 12
15
20 25 Outer iteration
30
35
39
44
Figure 4. Timing performance for the adaptive algorithm on problem dfl001.
cles in the gure. This is a savings of 38 factorizations compared to the PCx algorithm. In between refactorizations, the number of conjugate gradient iterations generally grows, more quickly for late values of than for earlier ones. We highlight the following observations from the gure, similar to the observations for the OB1-R code. The adaptive algorithm produces signi cant savings in the beginning stage, especially from the 2nd to the 10th value of . The frequency of reinitialization of the preconditioner grows as is decreased. The adaptive algorithm succeeds in keeping the cost at or better than the direct cost on all iterations but two. On those, the timings are close to drct cost. Figure 4 shows that Problem dfl001 has a long middle stage. While computing the predictor at 9 ; the adaptive algorithm does not converge to the prede ned tolerance within the maximum number of iterations allowed. The algorithm thus decides to carry out refactorization, resulting in high cost. The algorithm switches to the middle stage at 12 and remains there through 38 . The adaptive algorithm keeps the cost close to or better than the direct cost on
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
17
Problem: creb Use direct solver Estimated direct solver cost With reinitialized preconditioner With updated preconditioner
3
Time (seconds)
2
1.42
1
0
1
2
3 Outer iteration
4
5
Figure 5. Timing performance for the adaptive algorithm on problem cre-b.
all iterations in the middle stage except for 14 and 38 . The adaptive algorithm switches to late and ending stages at 39 and 45 respectively. Problem cre-b is not suitable for the iterative method. The adaptive algorithm discovers this at 5 and switches to the direct method (Figure 5). This happens because twice in a row the updated preconditioner is ine cient right after the preconditioner is reinitialized. Such behavior occurs in problems like cre-d, ken-11, ken-13, and osa-07. It is vitally important that the algorithm can make this decision automatically. Table 2 shows the computational results on the three problem sets, comparing the number of values needed by the interior point method, the relative duality gap in the nal answer, and the CPU time required by PCx and the adaptive algorithm. The last column is the di erence between the PCx and the adaptive times. A positive di erence means the adaptive algorithm is faster. We summarize the following observations from the results in table 2. Both PCx and the adaptive algorithm converge to solutions satisfying the optimality criteria de ned in PCx except on the problem greenbea, which is well-known to be di cult for interior point methods 33]. The optimal criteria
18
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
The more costly problems (PCx timing > 1000 seconds). PCx algorithm Adaptive algorithm
Time (seconds)
10
4
10
3
22
23
24
25 26 27 Test problems sorted by PCx timing
28
29
10
3
The cheaper problems (PCx timing <= 1000 seconds). PCx algorithm Adaptive algorithm
Time (seconds)
10
2
1
5
10 15 Test problems sorted by PCx timing
20 21
Figure 6. Timing performance comparison for all test problems, sorted by the time taken by the PCx algorithm.
include small duality gap, primal feasibility, and dual feasibility. The algorithms take a similar number of values and achieve similar duality gaps in most of the test problems. On some problems like pds-10, pds-20, degen3, and maros-r7, however, the adaptive algorithm takes 1 or 3 additional steps and achieves duality gaps several orders of magnitude smaller. In dfl001, the adaptive algorithm achieves the optimal criteria in 10 fewer steps, obtains a duality gap 1 order of magnitude smaller, and is faster than PCx by 880 seconds. If the total time for solution is small (i.e., 7 minutes or less), then the performance of the two algorithms is similar. On more costly problems such as dfl001, pilot87, net0832, and pds-20, the adaptive method is faster. Figure 6 compares the timing of the problems.
W. Wang, D. P. O'Leary / Adaptive Iterative Method in IPM
19 Di . -0.65 -2.59 880.04 0.59 -0.38 0.15 -0.60 7.36 -3.74
Problem NETLIB: d2q06c degen3 d 001 t2d greenbea maros-r7 pilot pilot87 stoch3 Kennington: cre-b cre-d ken-11 ken-13 ken-18 osa-07 osa-14 osa-30 osa-60 pds-06 pds-10 pds-20 Network: net0108 net0116 net0408 net0416 net0816 net0832 net0864 net1632
IPM ite. PCx Adap 29 16 58 23 48 18 36 34 31 40 40 21 26 30 25 27 27 30 37 41 55 16 18 20 19 20 20 20 22 29 19 48 23 50 19 36 35 31 40 40 21 26 29 25 27 26 31 36 44 58 16 19 21 20 21 21 21 23
Rel. dul. gap PCx Adap 2.53e-07 1.18e-08 5.63e-08 4.94e-07 1.90e-10 1.31e-08 2.88e-07 1.96e-07 6.13e-08 9.69e-07 1.21e-06 6.44e-08 5.42e-07 1.84e-06 2.55e-07 1....