A Mehrotra type predictor-corrector interior-point algorithm for linear programming

In this paper, we analyze a feasible predictor-corrector linear programming variant of Mehrotra's algorithm. The analysis is done in the negative infinity neighborhood of the central path. We demonstrate the theoretical efficiency of this algorithm by showing its polynomial complexity. The complexity result establishes an improvement of factor \begin{document}$ n^3 $\end{document} in the theoretical complexity of an earlier presented variant in [ 2 ], which is a huge improvement. We examine the performance of our algorithm by comparing its implementation results to solve some NETLIB problems with the algorithm presented in [ 2 ].


1.
Introduction. Karmarkar's publication [11] which appeared in 1984, meant a great moving forward in the area of optimization algorithms. Following this, a large amount of interior-point methods (IPMs) has been published; for example [4-8, 12, 17, 19, 23-25]. These algorithms have been used in different fields, such as technology, economics, transportation, statistics, etc.
After some first versions were developed either in a primal-only or dual-only settings, primal-dual IPMs were investigated. Meggido [16] proposed the skeleton of a symmetric primal-dual algorithm. Kojima, Mizuno and Yoshise [12] developed the first primal-dual path-following method. Their method works in the way that, at each iteration, the iterates are allowed to be located in a wide neighborhood of the so-called central path. They take only a single Newton step to do that. They were able to obtain an O (nL) complexity bound for a linear programming (LP) with n nonnegative variables and integer data with length L.
Predictor-corrector methods are one of the most studied variants of IPMs. The performance of practical primal-dual algorithms is significantly improved by the incorporation of higher-order information as in Mehrotra's predictor-corrector method [15]. Mehrotra's predictor-corrector algorithm [15] is an implementation of IPMs was proposed in 1989. The idea is to first compute a proper search direction based on a first order term (predictor). The step length that can be taken in this direction evaluates that how much centrality correction is needed. Then, a corrector term is computed which contains both a centrality term and a second order term. The final search direction is the sum of the predictor direction and the corrector direction. Original variants of Mehrotra's predictor-corrector algorithm [14,15] are among the main used algorithms in IPMs based software packages [3,9,28]. Almeida et al. [1], Salahi et al. [21] and Zhang and Zhang [26] proved the polynomial complexity of some Mehrotra type predictor-corrector variants. In [20], the authors have shown by a numerical example that a feasible version of the algorithm may be forced to make many small steps. This motivated them to introduce certain safeguards, that allowed them to prove polynomial iteration complexity. Afterwards, in [21] the author proposed a new modification of the algorithm that also enjoys polynomial iteration complexity and maintains its efficiency in practice. The theoretical study of complexity theory is one of the key points in LP and IPMs allow us to find an upper bound on the time needed by an algorithm to solve a given instance. The complexity of IPMs was analyzed at first by Kojima, Meggido, Mizuno [13] and Zhang [27]. Kojima et al. [13] studied the convergence of large neighborhood primal-dual interior-point algorithms for LP, showing that the gap between theory and practice is very small for those algorithms. Recently, Almeida and Teixeira [2] have presented a predictor-corrector interior-point algorithm for LP, based on the negative infinity wide neighborhood. They established the complexity bound O n 4 |log ε| for their algorithm and they proved its superlinear convergence. Motivated by the mentioned works, in this paper, we analyze a Mehrotra type predictor-corrector interior-point algorithm for LP. Our algorithm works based on the same used neighborhood and the same systems in [2] to determine the predictor and the corrector directions. But, by the system defining the predictor directions, we introduce a step length which fulfills some conditions. Then, we give a different analysis that leads to a complexity bound of O n log 1 ε iterations, which shows a huge improvement in the theoretical complexity bound of [2].
We introduce some notations used throughout the paper. R n + (R n ++ ) is the nonnegative (positive) orthant of R n . M m×n denotes the set of real m × n matrices. The p-norm is denoted by · p , especially the Euclidean norm is denoted as · . If x, s ∈ R n then xs denotes the componentwise product of the vectors x and s, and so is true for other operations, e.g., √ xs, 1 xs and x s . x i denotes the i−th component of the vector x. x min and x max mean the minimum and the maximum components of x, respectively. The nonnegative part of λ ∈ R is denoted as λ + and is defined as λ + := max{λ, 0}. Also, the nonpositive part of λ is λ − := min{λ, 0}. Correspondingly, for x ∈ R n we denote where the superscript T denotes the transpose of a vector in R n . e denotes the all ones vector in R n .
The outline of the remainder is as follows. In Section 2, after introducing the problem and systems defining the predictor and the corrector steps, we present our algorithm and its description. Section 3 is devoted to theoretical investigation of the algorithm. Through this section, we follow some steps to get the upper bound on the number of required iterations to capture an ε-optimal solution of the problem using our presented algorithm. The results will improve the theoretical complexity bound presented in the very similar algorithm in [2], by a factor n 3 . In Section 4, we compare our algorithm with the algorithm in [2] by solving some NETLIB problems via both algorithms. Finally, some concluding remarks are presented in Section 5.
2. The algorithm and its preliminaries. Consider the LP problem in the standard form: with dual variable vector u ∈ R m and dual slack vector v ∈ R n . Denote the closed set F as a collection of all points that satisfy the constraints of LP and its dual programming: and the open set F 0 , as a collection of all points that satisfy the constraints of LP and dual programming, and (x, v) are strictly positive: Throughout the paper, we assume the following: • A is a full rank matrix, Because of K.K.T optimality conditions, finding simultaneous optimal solutions for both above problems is equivalent with to find a nonnegative solution (x, v) of the following system The idea underlying primal-dual IPMs is to replace the complementarity condition xv = 0 by the nonlinear equation xv = µe, where µ > 0. Letting the interior-point condition (IPC) be satisfied for the primal and dual problems, i.e., F 0 = ∅, then, the new obtained system after the mentioned replacement, has a unique solution for a fixed µ > 0, called the µ-center or analytic center (Sonnevend [22]). The set of µ-centers for µ > 0 forms a well-behaved curve, called central path which plays an important role in convergence analysis of the IPMs. The central path converges to the primal-dual optimal solution as µ → 0. Our predictor-corrector algorithm uses the negative infinity norm neighborhood of the central path defined by where γ ∈ (0, 1) is a constant, and µ := x T v n is the normalized duality gap. The predictor direction of this algorithm, (∆x a , ∆u a , ∆v a ), is obtained by solving the Then, we compute the maximum step length α a ∈ [0, 1] which guarantees α a ∆x a ∆v a ≤ xv 2 and (x + ∆x a )(v + ∆v a ) ≥ µe.
Using the outcomes of (2) and (3), we obtain the corrector direction (∆x, ∆u, ∆v) by solving the system where the centering parameter µ is defined adaptively as µ := γ 1−γ µ. With α ∈ [0, 1] be the step length taken along (∆x, ∆u, ∆v), the new iterate is This algorithm can be formalized in Algorithm 1.

Algorithm 1.
Require: (x 0 , u 0 , v 0 ) ∈ N − ∞ (γ) with γ ∈ (0, 1 3 ) and an accuracy parameter ε > 0; Set k=0; Set Obtain the predictor direction (∆x a , ∆u a , ∆v a ) by solving (2); (b) Obtain the corrector direction (∆x, ∆u, ∆v) by solving (4); (c) Compute the step lengthᾱ k , We remark the following: 1. Note that using the self-dual embedding technique we can always construct an LP problem in such a way that the IPC holds. So, the IPC can be assumed without loss of generality. Furthermore, the self-dual embedding model yields x 0 = v 0 = e and y 0 = 0. 2. Note that our algorithm differs to the one presented in [2], in determining the corrector search directions. In system (4), we involve a factor α a , equipped with (3) for the quantity ∆x a ∆v a . We show that considering such a modified system to define the corrector directions, improves the complexity bound of the algorithm in [2] by a factor n 3 .
Also, noticing system (4), the normalized duality gap related to the new iterates is 3. Analysis of the algorithm. In this section, we investigate the analysis of Algorithm 1 to obtain its computational complexity. Proof. We have which proves the lemma.
Proof. We may write where the first inequality is due to the assumption that (x, v) ∈ N − ∞ (γ). This concludes the result of the lemma.
Proof. The proof is similar to the proof of Lemma 3.2 and therefore is omitted.
Proof. From (3) and noting that α a ∈ [0, 1], we have α 2 a ∆x a ∆v a ≤ xv. For the other inequality, using (3), we may write which concludes that xv + α 2 a ∆x a ∆v a ≥ 0 and completes the proof.
Proof. Defining D = diag ( x v ) and multiplying both sides of the third equation in (2) by ( √ xv) −1 , we obtain Using lemma 3.4, we have where the second inequality is derived using Lemma 3.1. The proof is completed.
, and D be as defined in Lemma 3.5, then where, Proof. By multiplying both sides of the third equation in system (4), similar to (8), we have Therefore, we may write where the second inequality is due to Lemmas 3.2 and 3.5. Using µ = γ 1−γ µ and noticing that γ ∈ (0, 1 3 ), we have Λ ≥ 2. The proof is completed. Now, based on the above obtained results, we determine a value for the required step length α in Algorithm 1.
Lemma 3.7. Let α be the step length taken in the algorithm, then α ≥ γ nΛ 2 .
4. Numerical results. In this section, we solve some LP problems from NETLIB by using our presented algorithm and the algorithm in [2]. In all cases, the initialization points for the algorithms are selected as x = e, v = e and y = 0 after applying the self-dual embedding scheme. The accuracy parameter ε is set to 10 −3 and we have chosen γ = 0.3344. Table 1 compares the number of iterations to solve the given problems via these algorithms. In this table, Alg 1 means our algorithm given in this paper and Alg 2 is the presented algorithm in [2]. It. stands for the number of required iterations to obtain an ε−optimal solution of the problems via each of the algorithms. The fifth and the seventh columns of the table show the value of x T v at the obtained ε−optimal solutions. As the fourth and the sixth columns of the table exhibit, our algorithm solves most of the given problems in a less number of iterations compared to the algorithm in [2]. 5. Concluding remarks. In this paper, we presented a Mehrotra type predictorcorrector algorithm for LP which is a variant of the same algorithm presented in [2]. Both algorithms follow the same neighborhood that is the negative infinity wide neighborhood of the central path for LP. We modified the system defining the corrector directions in [2] by introducing a step length α a equipped with (3). We analyzed the algorithm to obtain a suitable step length for deriving the new iterates of the algorithm. As is shown by Theorem 3.8, the analysis leads to an improvement of factor n 3 in the complexity analysis, which is a significant improvement in the theoretical complexity of the algorithm given in [2]. Besides this, Table 1 illustrates the results of a comparison between our algorithm and the algorithm in [2] to solve some sample problems from NETLIB. This comparison reveals that for most of the selected problems our proposed algorithm works practically better, as well.