USE OF CONTINUED ITERATION ON THE REDUCTION OF ITERATIONS OF THE INTERIOR POINT METHOD

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.


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 Gondzio [12].Among the variations of the interior point methods for the solution of linear programming problems, the predictor-corrector method presented by Mehrotra [15] 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 & Ye [19] 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 direction [15], other technique which determine correction directions have been developed over the last decades to improve interior point methods.Gondzio [5,11] presented the multiple centrality corrections to interior point methods that consist in computing a direction so that x i z i belongs to the neighborhood of the central path [20], thereby enhancing convergence.
Jarre & Wechs [13] 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 r a 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, (d x 1 , dy 1 , d z 1 ), . . ., (d x k , dy k , d z k ), generated for each r 1 , r 2 , . . ., r k , 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 & Li [16] follow the approach of Jarre & Wechs [13].They found multiple corrector directions by using information generated by a suitable Krylov subspace.The affine-scaling direction ( dx, d y, dz), the Mehrotra's corrector direction ( dx, d y, dz) and the first k directions (d x 1 , dy 1 , d z 1 ), . . ., (d x k , dy k , d z 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.
Numerical experiments have shown that direct methods provide sufficiently accurate solutions for interior point methods to achieve fast convergence regardless the ill-conditioning involved [1].Such results are also supported by theory [20].However, iterative methods have good accuracy when applied to wellconditioning systems or when preconditioners, as proposed by Oliveira & Sorensen [18], are used for an ill-conditioning system [12].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 as [12,13,16], using direct methods.Inspired by the original continued iteration proposal presented on [17] 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.

INTERIOR POINT METHODS
The linear programming problem in standard form is given by: The problem (1) is called primal problem.The dual problem, associated with it, is given by: in which y ∈ R m represents the vector of free dual variables and z ∈ R n represents the slack dual variables.

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: We use variables substitution to determine the solution of (4).By eliminating the variable d z = X −1 (r a − Z d x), we obtain the augmented system: , the following symmetric and positive definite normal equations system is obtained As AD A T 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 x i z i of r a 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 x i z i = μ, which is updated every iteration and tends to zero as the method approximates to a solution.

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 x i z i from r a 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 x i z i = dx i dz i , then a non-linear correction direction is computed.Consider, Dx = diag( dx) and Dz = diag( dz), the corrector direction d is the combination of the centering and the non-linear correction direction, which is obtained by solving: Thus, the (d) predictor-corrector direction, given by the sum of two directions found, i.e., d = d + d, can be interpreted as a solution of the linear system: in which r s = r a + μe − Dx Dze.Solving the system (7), the predictor-corrector direction is given by: Note that, the most expensive step is solving a linear system with the same matrix AD A T .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 (x 0 , y 0 , z 0 ) by computing a new point: (8) and k varies from 0 to a maximum number of iterations.

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].

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 AD A t , 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.

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 (d x, d z) 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.

Two blocking components
In the continued affine-scaling direction, in which d 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 d 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 Dikin [7] by developing the primal affine-scaling method, we find the solution of this system using an auxiliary problem to determine dx.The auxiliary problem is built using (10) and the first equation of ( 4).The component dx j corresponding to d z j is computed using the third equation of ( 4), in which the value must be −x j .We formulate a problem that determines d 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: in which D = Z −1 X, d x, d x ∈ R n , β a = 0, and β b = −x j .
Consider: (11) is quadratic, we find the solution to this problem using the Lagrangian function, which is given by: d Pesquisa Operacional, Vol.36( 3), 2016 From the last two equations of (4), and using least squares in the overdetermined system to compute d y, we find the remaining directions given by: Remark 3.2.The continued affine-scaling direction computation implies solving three linear systems with AD A T .We solve two linear systems: one to obtain v in (12) and another to get d 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 d we should have: Furthermore, the direction d should be an approximate solution of the linear system (7).We determine d 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 β b due to the computing of centering and non-linear correction.
If dz j = 0, we find the value of the component dx j corresponding to the last equation of (7).Then, the auxiliary problem, in this case, is given by: Solving the problem (15), it is obtained: We find the remaining directions in (7) by: Remark 3.3.To determine the direction d above it is not necessary to solve again the linear systems involving the matrix AD A t to compute v.In fact, it is necessary to solve only a linear system with matrix AD A t to determine d y.
Remark 3.4.The computing of direction d 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 d 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 Thus, αP + α P ∈ (0, 1] and αD + α D ∈ (0, 1].In [3], there is no limitation on the new stepsizes.

Criteria for using continued iteration
Consider r = (r p , rd , ra ) t , formed by new residuals produced by the continued iteration e r = (r p , r d , r a ) 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 ||r|| 2 < ω 1 ||r|| 2 in which ω 1 ∈ (0, 1).This is an improvement over the implementation described in [3].

CONTINUED ITERATION WITH BOUNDED VARIABLES
Consider the primal linear programming problem with bounded variables in its standard form: in which s ∈ R n represents the slack of bounded variables, and u ∈ R n represents the upper bound of x.
The dual problem in its standard form associated with ( 16) is given by: in which w ∈ R n 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).

Computing the continued predictor-corrector direction
In this case, it is conducted the blocking components and they are given by the directions (d x, ds, d z, dw).
The blocking components i and j in directions d x and d z 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 , then in the new direction d x i = 0. Otherwise, we assign i ← l 1 and ds i = 0.If g = −z k j /d z k j , then in the new direction dz j = 0; otherwise, j ← l 2 and dw j = 0.

Two blocking components
With the given components, we have computed the continued affine-scaling direction d 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 dx, ds, dz, dw ∈ R n and d y ∈ R m .
To determine the direction d, initially, we compute d x analogously to (3.1.1),using an auxiliary problem, which in this case is given by: and 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 proposed method.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 profile [8] 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.

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.

Figure 2 -
Figure 2 -Total time for solving problems.

Table 2 -
Comparison between PCx and PCx-IC.