ABSTRACT
Interior point methods have been widely used to determine the solution of large-scale linear programming problems. The predictor-corrector method stands out among all variations of interior point methods due to its efficiency and fast convergence. In each iteration it is necessary to solve two linear systems to determine the predictor-corrector direction. Solving such systems corresponds to the step which requires more processing time, and therefore, it should be done efficiently. The most common approach to solve them is the Cholesky factorization. However, Cholesky factorization demands a high computational effort in each iteration. Thus, searching for effort reduction, the continued iteration is proposed. This technique consists in determining a new direction through projection of the search direction and it was inserted into PCx code. The computational results regarding medium and large-scale problems, have indicated a good performance of the proposed approach in comparison with the predictor-corrector method.
Keywords: interior point methods; linear programming; continued iteration
1 INTRODUCTION
Linear programming is an optimization technique widely used due to its robustness, the efficiency of its algorithms and the large number of problems that can be formulated through it. The study of interior point methods for linear programming has been one of the most active areas of research in optimization in recent decades because of great advances described by Gondzio12. Among the variations of the interior point methods for the solution of linear programming problems, the predictor-corrector method presented by Mehrotra15 has great prominence due to its efficiency and fast convergence. This is a primal-dual method that needs to solve two linear systems with the same symmetric positive definite matrix in each iteration. To solve these systems the most used approach is the Cholesky factorization. Thus, interior point methods present a critical step in its computational performance, since this factorization in most cases demands a large computational effort. Vavasis & Ye19 propose an algorithm that follows the central path, using short steps or a layered least squares step, which return an exact optimum after a finite number of steps. Besides Mehrotra’s correction direction15, other technique which determine correction directions have been developed over the last decades to improve interior point methods. Gondzio5), (11 presented the multiple centrality corrections to interior point methods that consist in computing a direction so that xizi belongs to the neighborhood of the central path20, thereby enhancing convergence.
Jarre & Wechs13 propose to generate corrector directions making use of recursion in a modification of Mehrotra’s corrector direction. The idea is to change the residual vector ra for another vector r in the linear system which computes the predictor direction. They question what the best choice would be for r to get a good search direction. In general, it is not easy to find a good value for r. So, they propose to find a subspace spanned by k different directions, (dx 1, dy 1, dz 1), ..., (dx k, dy k, dz k), generated for each r 1, r 2, ..., rk , where each one of these directions is a predictor-corrector type. They formulate a small linear subproblem to find the best combination of weights in the generated directions to increase the stepsize and the reduction of the complementarity gap.
Mehrotra & Li16 follow the approach of Jarre & Wechs13. They found multiple corrector directions by using information generated by a suitable Krylov subspace. The affine-scaling direction k directions (dx 1, dy 1, dz 1), ..., (dx k, dy k, dz k), generated by Krylov subspace are combined with appropriate weights. Thus a linear programming problem is solved to determine the best set of weights in the combined search direction. and the first , the Mehrotra’s corrector direction
Numerical experiments have shown that direct methods provide sufficiently accurate solutions for interior point methods to achieve fast convergence regardless the ill-conditioning involved1. Such results are also supported by theory20. However, iterative methods have good accuracy when applied to well-conditioning systems or when preconditioners, as proposed by Oliveira & Sorensen18, are used for an ill-conditioning system12. The direct methods take advantage of the factorization already computed to determine the solution of additional linear systems, but using iterative methods this advantage is lost.
We present in Section 3 an alternative technique to improve the predictor-corrector interior point method as12), (13), (16, using direct methods. Inspired by the original continued iteration proposal presented on17 for the dual affine-scaling method, we present the continued iteration for the predictor-corrector interior point method. The technique is incorporated in the predictor-corrector method in order to reduce the total number of iterations required to solve a linear programming problem.
The continued iteration consists in determining a new projection direction. Additionally, to take advantage of the Cholesky factorization already computed, the new direction is found by minimizing the deviation from the original rescaled direction. This technique is used when an iteration of the predictor-corrector method is completed.
This paper is organized as follows: In Section 2, we define the linear programming problem and describe, briefly, the primal dual interior point method. In Sections 3 and 4, the continued iteration is presented and incorporated into the predictor-corrector interior point method for linear programming problems, whether with bounded variables or not. Section 5 shows the results of computational experiments conducted with several problems selected from different collections. Finally, in Section 6, the conclusions are presented.
2 INTERIOR POINT METHODS
The linear programming problem in standard form is given by:
in which A∈Rm ∈× n rank (A) = m, x ∈ Rn , c ∈ Rn and b.∈ Rm
The problem (1) is called primal problem. The dual problem, associated with it, is given by:
in which y ∈ Rm represents the vector of free dual variables and z ∈ Rn represents the slack dual variables.
The first order optimality conditions (Karush-Kuhn-Tucker) of problems (1) and (2) are given by20:
in which X=diag(x), Z=diag(z) and e ∈ Rn such that e = (1.1...,1)T
2.1 Primal-Dual Affine-Scaling Interior Point Method
The primal-dual affine-scaling method entails applying Newton’s method to F(x,y,z)=0, which is formed by the optimality conditions (3), ignoring (x,z) ≥ 0, but starting from an interior point (x,z) > 0. Thus, the primal and dual problems are solved simultaneously. The affine-scaling direction is obtained by solving the linear system given below:
in which rp = b - Ax, rd = c - Aty - z and ra = - XZe.
We use variables substitution to determine the solution of (4). By eliminating the variable dX -1(ra - Zx) we obtain the augmented system: =
in which D = XZ -1 and r 1 = rd -X -1 ra . Now eliminating x = D(Aty - r1), the following symmetric and positive definite normal equations system is obtained
As ADAT is symmetric positive definite, the most used approach for solving system (5) is the Cholesky factorization(9).
The primal αP and dual αD stepsize, to keep the interior point in each iteration are given by:
in which τ ∈ (0,1).
Remark 2.1. In the primal-dual affine-scaling method, the products xizi of ra can converge to zero at different speeds. Thus, the method may fail or progress slowly. To avoid this, the parameter (μ) is added to the optimality conditions so that xizi = μ, which is updated every iteration and tends to zero as the method approximates to a solution.
2.2 Predictor-Corrector Interior Point Method
The predictor-corrector method presented by Mehrotra(15) consists in using a direction that contains three components: the affine-scaling direction, the centering direction and the non-linear correction direction. The affine-scaling direction is determined by (4). The centering direction is computed to prevent products xizi from ra converging to zero at different speeds as mentioned in Remark (2.1). Mehrotra(15) presents a proposal for computing μ.
When αP = αD = 1 the primal and dual constraints are satisfied, but the complementary products have xizi = z=diag(xix = diag(z), the corrector direction zi then a non-linear correction direction is computed. Consider, is the combination of the centering and the non-linear correction direction, which is obtained byx) and
Thus, the (d) predictor-corrector direction, given by the sum of two directions found, i.e., d = , can be interpreted as a solution of the linear system: +
in which rs = ra + μe -ze Solving the system (7), the predictor-corrector direction is given by:x
Note that, the most expensive step is solving a linear system with the same matrix ADAT. The Cholesky factorization is used and it has been already computed to obtain the affine-scaling direction.
Therefore the predictor-corrector method is an iterative method that starts from an interior point (x0, y0, z0) by computing a new point:
in which (dxk, dyk, dzk) is given in (8) and k varies from 0 to a maximum number of iterations.
2.2.1 Stopping criteria
Consider εp, εd, εo primal feasibility, dual feasibility and optimality tolerances, respectively, then the stopping criteria is given by:
For further details on the predictor-corrector interior point methods see(14), (15), (19).
3 CONTINUED ITERATION
The continued iteration is proposed in order to reduce the total number of predictor-corrector interior point method iterations, and consequently, reduce the total computational time. This iteration determines a new direction, after computing the predictor-corrector step, without being necessary to perform a new Cholesky factorization. The new direction is determined by setting some of its components to zero and thus, the deviation from the predictor-corrector direction is kept to a minimum. Furthermore, the Cholesky factorization of ADAt already computed in the predictor-corrector method iteration, is used again for solving the linear systems involved. Therefore, the computing of this new direction is not very expensive in comparison with a new predictor-corrector method iteration. This new direction is called continued predictor-corrector direction.
3.1 Computing the continued predictor-corrector direction
The continued predictor-corrector direction is computed in an analogous way to the computing of the predictor-corrector direction. Initially, the continued affine-scaling direction is computed, after that the centering and nonlinear correction direction is computed. More details on these directions will be seen later, because before that we need to define the blocking components.
The components that must be fixed at zero, in the new direction, are those responsible for blocking in the previous direction. The blocking components in the directions (dx,dz) are given by:
As blocking component in any or both directions may not exist, some likely cases can be considered: two, one or none blocking component.
Remark 3.1. If there is more than one blocking component i or j, only the first blocking determined component is considered in this study. When there is no blocking component, the continued iteration is not performed.
In the original proposal of continued iteration to the predictor-corrector method only one of the blocking components was fixed at zero (3). In this work, the two blocking components (i and j) are fixed at zero.
3.1.1 Two blocking components
In the continued affine-scaling direction, in which is to be computed, the components that perform blocking in the previous direction are fixed at zero, thus, we get:
which are additional conditions to the linear system (4). We compute approximately, since linear system (4) admits a unique solution. In addition, we are imposing new conditions (10) to such system.
Similarly to the idea of Dikin7 by developing the primal affine-scaling method, we find the solution of this system using an auxiliary problem to determine x rescaling the new direction and minimizing the deviation in relation to the original rescaled direction, so that the Cholesky factorization is already computed and used. The problem is given by:xj corresponding to x. The auxiliary problem is built using (10) and the first equation of (4). The component j is computed using the third equation of (4), in which the value must be -xj. We formulate a problem that determines
in which D = Z-1X, x, dx ∈ Rn, βa =0, and βb = - xj.
Consider: θi = ATi(ADAT)-1Ai, θj=ATj(ADAT)-1Aj and θij = ATi (ADAT)-1Aj.
Since the problem (11) is quadratic, we find the solution to this problem using the Lagrangian function, which is given by:
in which
From the last two equations of (4), and using least squares in the overdetermined system to compute y, we find the remaining directions given by:
Remark 3.2. The continued affine-scaling direction computation implies solving three linear systems with ADAT. We solve two linear systems: one to obtain v in (12) and another to get y in (13).
Given the continued affine-scaling direction, when performing the centering and non-linear correction, we obtain continued predictor-corrector direction, similarly as done in the predictor-corrector method (2.2). Keeping the components that block at zero, in direction we should have:
Furthermore, the direction βb due to the computing of centering and non-linear correction. If zj = 0 we find the value of the component should be an approximate solution of the linear system (7). We determine xj corresponding to the last equation of (7). Then, the auxiliary problem, in this case, is given by: in a similarly way to determine the continued affine-scaling direction described above. However in auxiliary problem (11) there is an change in the value
in which D = Z-1X, x, dx ∈ Rn, βa = 0 and βb = -xj + μ/zj
Solving the problem (15), it is obtained:
We find the remaining directions in (7) by:
Remark 3.3. To determine the direction ADAt to compute v. In fact, it is necessary to solve only a linear system with matrix ADAt to determine y. above it is not necessary to solve again the linear systems involving the matrix
Remark 3.4. The computing of direction is analogous whether there is only one or two blocking components. That is, in (10) and (14) only one condition is considered. In the auxiliary problems (11) and (15), the restriction that has no blocking component is disregarded.
Remark 3.5. As the direction is determined in such a way that the change from the direction d is as low as possible, the stepsizes P and D adopted are
and
Thus, 3), there is no limitation on the new stepsizes.D + αD ∈ (0,1⌋ In(P + αP ∈ (0, 1⌋ and
3.2 Criteria for using continued iteration
Consider 3)., formed by new residuals produced by the continued iteration e r = (rp,rd, ra)t. Given that the primal infeasibility, the dual infeasibility and/or the complementary products are reduced with continued iteration, then the continued direction only is used if ||||2 < ω1||r||2 in which ω1 ∈ (o, 1) This is an improvement over the implementation described in(
4 CONTINUED ITERATION WITH BOUNDED VARIABLES
Consider the primal linear programming problem with bounded variables in its standard form:
in which s ∈ Rn represents the slack of bounded variables, and u ∈ Rn represents the upper bound of x.
The dual problem in its standard form associated with (16) is given by:
in which w ∈ Rn represents the dual variables associated with the slack variables s.
The first order optimality conditions (Karush-Kuhn-Tucker) of the problems (16) and (17) are as follows:
in which S=diag(s) and W=diag(w).
4.1 Computing the continued predictor-corrector direction
In this case, it is conducted the blocking components and they are given by the directions (dx, ds, dz, dw). The blocking components i and j in directions dx and dz are defined as in (9). Consequently, the blocking components in the directions ds, dw are defined by:
The computing for determining the null components in the new direction is performed as follows.
Consider h = min{- xki/dxki, - skl1/ dskl1} and g = min{- zkj/dzkj, - wkl2/dwkl2}. If h = -x ki/dxki, then in the new direction xi = 0 Otherwise, we assign i ← l1 and si If =-zkj/dzkj then in the new direction zj = 0; otherwise, j ← l2 and wj = 0.
4.1.1 Two blocking components
With the given components, we have computed the continued affine-scaling direction so that:
Additionally, the continued affine-scaling direction should be an approximate solution of the linear system obtained when Newton’s method is applied to the optimality conditions of the problem, being the corresponding linear system given by:
in which s, w ∈ Rn and z, y ∈ Rm.x,
To determine the direction , initially, we compute x analogously to (3.1.1), using an auxiliary problem, which in this case is given by:
in which D = (XZ-1 + SW-1)-1,
and
Problem (20) is similar to the problem (11), in which only the values βa e βb may change. Thus,
For the remaining equations (19), we obtain:
The procedure performing centering and non-linear correction is similar to the case without bounds, thus obtaining the continued predictor-corrector direction.
Remark 4.1.When there is one blocking component, the direction computation is similar to (4.1.1). Furthermore, the equations related to the missing blocking component are disregarded in (18) and (20).
Remark 4.2.The criterion for using the continued iteration with bounded variables is similar to that given in (3.2); in this case, ,and r = (rp, ru, rd, ra,rb )t
5 COMPUTATIONAL EXPERIMENTS
The computational experiments were performed on an Intel Core i7, 8 GB RAM, 1TB HD and Linux Operating System.
The continued iteration was implemented in language C and incorporated into the PCx code6 with multiple centrality corrections turned off11. To analyze the performance of PCx with the proposed approach, the computational experiments were performed with free access problems belonging to Netlib10, Kennington4 and Qaplib2 collection.
The PCx has a good behavior near a solution, thus the continued iteration was not used in the last iterations. Based on criteria defined in (2.2.1), the continued iteration is to be applied if:
in which ω 2 ∈ (0, 1).
For the computational experiments ω 1 = 10-7 in (3.2) and ω 2 = 0.99.
Table 1 shows, in the first column, the name of the selected problems and in the second and third columns, are the dimension of the problems (number of rows and columns) after the preprocessing of the PCx. In the fourth column, there is the density of the Cholesky factor provided by PCx, which performs a reordering to obtain this sparse factor. The fifth column relates the number of variables with bounds (VB) for each problem. The last column indicates the collection to which the problem belongs.
Table 2 shows the results by incorporating the continued iteration into PCx code. The first column shows the tests problems. The second column presents the results obtained by PCx standard version without multiple centrality corrections. The third column has the results obtained by PCx with the continued iteration, which will be referred to as PCx-IC. The columns PCx and PCx-IC indicate the total number of iterations (k) performed by the interior point methods and running time in seconds (t). In addition, PCx-IC has thenumber of performed continued iterations (ic).
In Table 2, we find 18 problems in which the number of iterations is reduced by PCx-IC. The PDS problems present the largest reductions in both of iterations and time by the proposedmethod. The OSA-60 problem has the highest iteration reduction ratio at 45%. Particularly in terms of time, the largest time reduction ratio stands at 7% in the PDS70 problem among the 14 problems displaying time reduction. As for the problems that require more processing time whenever PCx-IC reduces the number of iterations the time is also reduced. However, there is no time reduction in the problems solved by PCx in a few seconds and with low reduction. This occurs due to the increased effort exerted by the continued iteration. Thus, excellent performance is obtained by PCx-IC in most large-scale problems regarding iterations reduction as well as for total computational effort.
Figure 1 shows the performance profile8 by PCx and PCx-IC regarding the total number of iterations. Note that the method PCx-IC offers better performance compared with PCx. It solves 78% of the test problems with a smaller number of iterations.
Figure 2 shows the performance profile of the methods considering the total running time as performance measure. In this case, PCx is the most efficient for approximately 60% of the test problems solved in less time. Moreover, both methods are robust for solving all tested problems.
6 CONCLUSIONS
In this study, the continued iteration has been proposed for the predictor-corrector interior point method with and without bounded variables. The main objective is to reduce the total number of predictor-corrector interior point method iterations and consequently, the total computational time. In the continued iteration, a new direction is determined, considering that some of its components are null according to the blocking component in the last computed direction. The continued predictor-corrector direction is determined by setting to zero some of its components and thus the deviation from the predictor-corrector direction is kept to a minimum. The Cholesky factorization, already computed in the predictor-corrector method iteration, is used again for solving the linear systems involved. Thus, it avoids computing a new Cholesky factorization. Then, the new method is applied on two levels. On the outer level, the Cholesky factorization and the traditional predictor-corrector direction are computed. On the inner level, a new direction employing the existing Cholesky factorization on the outer level is used by the continued iteration. The computational effort of each continued iteration is dominated by the solution of up to four linear systems with the same matrix already factored. The computational results show that the continued iteration reduces the number of iterations in most of the test problems. The problems that the PCx needs more time to solve also have their time reduced by PCx-IC. Although the continued iteration with the proposed direction has reduced the number of iterations in most problems, the extra effort per iteration to compute the new direction in smaller problems does not affect the reduction of total computational time. In future work other continued directions will be searched aiming to reduce the extra effort that such directions add in the predictor-corrector method.
ACKNOWLEDGEMENTS
This research was sponsored by the Brazilian Agencies CAPES, CNPq and FAPESP.
REFERENCES
- 1 ANDERSEN ED, GONDZIO J, MÉSZÁROS C & XU X. 1996. Implementation of interior point methods for large scale linear programming. In: Interior Point Methods in Mathematical Programming, TERLAKY T (Ed.), Kluwer Academic Publishers, pp. 189-252.
- 2 BURKARD RE, KARISCH S & RENDL F. 1991.QAPLIB - A quadratic assignment library problems. European Journal of Operational Research, 55(1): 115-119.
- 3 BERTI LF & OLIVEIRA ARL. 2011. Iterações continuadas aplicadas ao método de pontos interiores preditor corretor. Anais do XLIII Simpósio Brasileiro de Pesquisa Operacional, 2328-2338 (in Portuguese).
- 4 CAROLAN WJ, HILL JE, KENNINGTON JL, NIEMI S & WICHMAN SJ. 1990. An empirical evaluation of the KORBX algoritmos for military airlift applications. Operations Research, 38(2): 240-248.
- 5 COLOMBO M &GONDZIO J. 2008. Further development of multiple centrality correctors for interior point methods. Computational Optimization and Applications, 41: 277-305.
- 6 CZYZYK J, MEHROTRA S, WAGNER M &WRIGHT SJ. 1999. PCx an interior point code for linear programming. Optimization Methods and Software, 11-12: 397-430.
- 7 DIKIN II. 1967. Iterative solution of problems of linear and quadratic programming. Soviets Mathematics Doklady, 8: 674-675.
- 8 DOLAN ED & MORÉ JJ. 2003. Benchmarking optimization software with performance profiles. Math. Program., Serie A, 91(2): 201-213.
- 9 DUFF IS, ERISMAN AM & REID JK. 1986. Direct Methods for SparseMatrices. New York: Oxford University Press, Inc.
- 10 GAY DM. (1985). Electronic mail distribution of linear programming test problems, Mathematical Programming Society Committee on Algorithms COAL Newsletter, 13: 10-12.
- 11 GONDZIO J. 1996. Multiple centrality corrections in a primal-dual methods for linear programming. Computation Optimization and applications, 6: 137-156.
- 12 GONDZIO J. 2012. Interior point methods 25 years later. European Journal of Operational Research, 218: 587-601.
- 13 JARRE F &WECHS M. 1999. ExtendingMehrotra’s corrector for linear programs. Advanced Modeling and Optimization, 1: 38-60.
- 14 LUSTIG IJ, MARSTEN RE & SHANNO DF. 1992. On Implementing Mehrotra’s Predictor-Corrector Interior Point Method for Linear Programming. SIAM Journal on Optimization, 2: 435-449.
- 15 MEHROTRA S. 1992. On implementation of a primal-dual interior point method. SIAM Journal on Optimization, 2(4): 575-601.
- 16 MEHROTRA S & LI Z. 2005. Convergence conditions and Krylov subspace based corrections for primal-dual interior-point method. SIAM Journal on Optimization, 15: 635-653.
- 17 OLIVEIRA ARL, LYRA C & CORREIRA PB. 1988. Implementação computacional de algoritmo polinomial de programacão linear: aplicação ao planejamento de operação de sistemas hidrotérmicos. Annals of the VII Congresso Brasileiro de Automática - CBA, São José dos Campos, Brazil, 928-929 (in Portuguese).
- 18 OLIVEIRA ARL & SORENSEN DC. 2005. A new class of preconditioners for large-scale linear systems from interior point methods for linear programming. Linear Algebra and its Applications, 394: 1-24.
- 19 VAVASIS SA & YE Y. 1996. A primal-dual interior point method whose running time depends only on the constraint matrix. Mathematical Programming, 74: 79-120.
- 20 WRIGHT SJ. 1997. Primal-dual interior-point methods. Philadelphia: SIAM Publications, Inc.
Publication Dates
-
Publication in this collection
Sep-Dec 2016
History
-
Received
27 May 2015 -
Accepted
07 Nov 2016