A corrector–predictor interior-point method with new search direction for linear optimization

We introduce a feasible corrector–predictor interior-point algorithm (CP IPA) for solving linear optimization problems which is based on a new search direction. The search directions are obtained by using the algebraic equivalent transformation (AET) of the Newton system which defines the central path. The AET of the Newton system is based on the map that is a difference of the identity function and square root function. We prove global convergence of the method and derive the iteration bound that matches best iteration bounds known for these types of methods. Furthermore, we prove the practical efficiency of the new algorithm by presenting numerical results. This is the first CP IPA which is based on the above mentioned search direction.

1 Introduction Karmarkar (1984) presented the first projective IPA for solving LO problems with polynomial-time complexity.After that, several results related to this theory have been published.The theory and practice of IPAs can be found in monographs written by Roos et al. (1997), Wright (1997), Ye (1997) and Nesterov and Nemirovski (1994).The IPAs for LO can be classified in multiple ways.One way to classify these algorithms is based on the step length.In this way, we can distinguish between short-and long-step IPAs.In theory, short-update algorithms give usually more efficient theoretical results with simpler analysis, while in practice, the large-step versions perform generally better.Another way to categorize the IPAs is related to feasibility of the iterates, hence, we can consider feasible and infeasible IPAs.Another type of IPAs that have proven to be efficient in practice are predictor-corrector IPAs.These algorithms consist of iterations using two types of steps, one predictor step and one or more corrector steps.For further details on classification of different IPAs see Illés and Terlaky (2002), Roos et al. (1997), Wright (1997) and Ye (1997).
The determination of the search directions plays a key role in the theory of IPAs.The most widely used technique for obtaining the search directions is based on barrier functions.By considering self-regular functions, Peng et al. (2002) reduced the theoretical complexity of large-step IPAs.Darvay (2002) introduced a new technique for finding search directions in case of these algorithms, namely the algebraic equivalent transformation of the system which defines the central path.Central path has been introduced independently by Sonnevend (1985,1986) and by Megiddo (1989).The importance of the central path in the literature of IPAs can be highlighted by the fact that it is unique, see Roos et al. (1997), Terlaky (2001), Wright (1997) and Ye (1997).Sonnevend (1985Sonnevend ( , 1986) proved that the central path led to a unique optimal solution called the analytic center.General idea of IPA is to approximately trace the central path and to compute an interior point, called ε-optimal solution, that well approximates the analytic center.From an ε-optimal solution with small enough ε > 0, using the so-called rounding procedure (Illés and Terlaky 2002;Roos et al. 1997), an optimal solution can be computed in strongly polynomial time.
The function of AET is applied to both sides of the nonlinear equation of the system that defines the unique central path.After the transformation the central path remains unique, and the Newton's method is applied to the transformed system in order to determine the displacements.In the literature, the most widely used function for AET is the identity map, namely in most of IPAs the central path has not been transformed.In the papers of Darvay (2002Darvay ( , 2003) ) (Darvay 2002(Darvay , 2003;;Darvay et al. 2016) achieve usually the best known iteration bounds, alongside many other IPAs (Roos et al. 1997;Wright 1997;Ye 1997).Further research is needed to investigate whether IPAs that use AETs are more efficient than IPAs that are not.It would be also interesting to identify some class of LO problems for which application of AET based IPA may be beneficial.
The first PC IPA has been independently developed by Mehrotra (1992) and Sonnevend et al. (1990).The PC IPAs consist of a predictor and several corrector steps in a main iteration.The aim of the predictor step is to approach the optimal solution of the problem in a greedy way.The usual consequence of the greedy predictor step is that the obtained strictly feasible solution no longer belongs to the given neighborhood of the central path.The goal of the corrector steps is to return the iterate back to the designated neighborhood.Mizuno et al. (1993) proposed the first PC IPA which uses only one corrector step in a main iteration.Darvay (2005Darvay ( , 2009) ) introduced PC IPAs for LO that are based on the AET technique and he used the function ψ(t) = √ t with domain D ψ = (0, ∞) in order to determine the transformed central path and the modified Newton system.The unique solution of this system led to a new search direction.Kheirfam (2016Kheirfam ( , 2015) ) proposed CP IPAs for convex quadratic symmetric cone optimization and second-order cone optimization, respectively.
Before summarizing the structure and results of this paper, it is worthwhile to mention that IPAs for LO have been extensively generalized to linear complementary problems (LCPs) (Cottle et al. 1992;Illés et al. 2010a, b;Kojima et al. 1991;Lešaja and Roos 2010;Potra and Sheng 1996;Yoshise 1996).There are several generalizations of the Mizuno-Todd-Ye PC IPA (Mizuno et al. 1993) from LO to sufficient LCPs like that of Potra (2002), and Illés and Nagy (2007).Recently, Potra (2014) published a new PC IPA for sufficient LCPs using wide neighborhood with optimal iteration complexity.The AET method for determining search directions for IPAs has been also extended to LCPs (Achache 2010;Asadi and Mansouri 2013;Kheirfam 2014;Wang et al. 2009) and LCPs over symmetric cones (Asadi et al. 2017a, b;Mohammadi et al. 2015;Wang 2012).
In this paper, a new CP IPA for LO is introduced.We use the AET method for the system which defines the central path that is based on the function ψ(t) = t − √ t.Newton's method is then applied to the transformed system in order to find the search directions.The analysis of the algorithm is more complicated with this function.Nevertheless, we were able to prove global convergence of the method and derive the iteration bound that matches best-known iteration bound for these types of methods.We also present some numerical results and we compare our CP IPA with the classical primal-dual method, which is based on the same search direction and uses only one step in each iteration.
The paper is organized as follows.In Sect.2, the primal-dual LO problem and the main concepts of the AET of the system defining the central path are given.In the following section the new CP IPA is presented.Section 4 contains the analysis of the proposed algorithm, while in Sect. 5 the iteration bound for the algorithm is derived.In Sect.6, we provide some numerical results that prove the efficiency of this algorithm.in the last Section some concluding remarks are provided.

Preliminaries
Consider the LO problem in the standard form and its dual problem We assume that the interior-point condition (IPC) holds for both problems; that is, there exists (x 0 , y 0 , s 0 ) such that Using the self-dual embedding model presented by Ye et al. (1994), Roos et al. (1997) and Terlaky (2001) we conclude that the IPC can be assumed without loss of generality.
In this case, the all-one vector can be considered as a starting point.Under the IPC, finding an optimal solution of the primal-dual pair is equivalent to solving the following system (1) The main idea of primal-dual IPAs is to replace the third equation in ( 1), the so-called complementarity condition for (P) and (D), by the perturbed equation xs = μe with μ > 0. Hence, we obtain the following system of equations: (2) It is proved in Roos et al. (1997) that there is a unique solution (x(μ), y(μ), s(μ)) to the system (2) for any μ > 0, assuming the IPC holds.The set of all such solutions constructs a homotopy path, which is called the central path (see Megiddo 1989;Sonnevend 1986).If μ tends to zero, then the central path converges to the optimal solution of the problem.
In what follows, we recall the AET introduced by Darvay et al. ( 2016) for LO that leads to calculation of new search direction for IPAs.For this purpose, we consider the continuously differentiable function ψ : R + → R + , and assume that its inverse ψ −1 exists.Note that the system (2) can be rewritten in the following form: where ψ is applied componentwisely.Applying Newton's method to system (3) for a strictly feasible solution (x, y, s) produces the following system for search direction (Δx, Δy, Δs) Defining the scaled search directions as one easily verifies that system (4) can be written in the form where .
For analysis of our algorithm, we define a norm-based proximity measure δ(x, s; μ) as follows: which has been considered for feasible IPAs for the first time in Darvay et al. (2016).Considering (6) we have d T x d s = 0. Thus, the vectors d x and d s are orthogonal.Using (8) and v > 0, one can easily verify that Hence, the value of δ(v) can be considered as an appropriate measure for the distance between the given triple (x, y, s) and (x(μ), y(μ), s(μ)).Moreover, note that if ψ(t) = t then p v = v −1 − v and we obtain the standard proximity measure Roos et al. (1997).From Darvay (2002Darvay ( , 2003)).Let Then, the orthogonality of the vectors d x and d s implies As an effect of this relation, we can also express the proximity measure using q v , thus Furthermore, holds.
The lower and upper bounds on the components of the vector v are given in the following lemma.

Corrector-predictor algorithm
In this section, we present a CP path-following algorithm for LO problems based on Darvay et al.'s idea.For this purpose, we define an τ -neighborhood of the central path as follows: where 0 < τ < 1.The algorithm begins with a given strictly feasible primal-dual solution (x 0 , y 0 , s 0 ) ∈ N (τ ).If for the current iterate (x, y, s), nμ > , then the algorithm calculates a new iterate by performing corrector and predictor steps.In the corrector step, we define v = xs μ , Ā = A diag( x v ) and we obtain the scaled search directions d x and d s by solving (6) with p v given in (7), namely 2v−e . (11) Newton directions of the original system (4), i.e., Δx = x v d x , Δs = s v d s can be expressed easily and the corrector iterate is obtained by a full-Newton step as follows: In the predictor step, we define and we obtain the search directions d (12) Note that the right-hand side of the system (12) is inspired by the predictor step proposed in Darvay (2005).Similarly to Δx and Δs, we define and the predictor iterate is obtained by where θ ∈ (0, 1 2 ) and also μ p = (1 − 2θ)μ.At the beginning of the algorithm, we assume that (x 0 , y 0 , s 0 ) ∈ N (τ ).We would like to determine the values of τ and θ in such a way that after a corrector step (x + , y + , s + ) ∈ N (ω(τ )) (where ω(τ ) < τ will be defined later) and after a predictor step (x p , y p , s p ) ∈ N (τ ).The algorithm repeats corrector and predictor steps alternatively until x T s ≤ is satisfied.A formal description of the algorithm is given in Fig. 1.

Analysis of the algorithm
The following technical lemma introduced by Wright (1997), is a generalization of Lemma C.4 (first u−v Lemma) in Roos et al. (1997).We will use it to estimate the norm of the product of scaled search directions.

The predictor step
The next lemma will prove the strict feasibility after a predictor step.
Proof For each 0 ≤ α ≤ 1, denote x p (α) = x + + αθΔ p x and s p (α) = s + + αθΔ p s. Therefore, using the third equation in (12), we obtain From ( 13), it follows that The second inequality is due to the fact that f (α) := α 2 θ 2 1−2αθ is monotonically increasing with respect to α; that is, f (α) ≤ f (1).The third inequality follows from Lemmas 1 and 2. The second equality can be derived from the third equation of ( 12).The inequality before the last line follows from the upper bound given in Lemma 1.
We define 123 It follows from ( 13), with α = 1, that and Lemma 4 Let (x + , y + , s + ) be a strictly feasible primal-dual solution and μ p = (1 − 2θ)μ with 0 < θ < 1 2 .Moreover, let h(δ + , θ, n) > 1 4 and assume that (x p , y p , s p ) denotes the iterate after a predictor step.Then v p > 1 2 e and , from ( 16) we have min v p 2 ≥ 1 4 , which yields v p > 1 2 e.Moreover, from Lemma 3 we deduce that the predictor step is strictly feasible; x p > 0 and s p > 0. Now, by the definition of proximity measure, we have where the first two inequalities are due to Lemma 5.2 in Darvay et al. ( 2016) and ( 16), respectively.The last equality follows from (15) and the last inequality is due to the triangle inequality.
We will give an upper bound for e − v + 2 .Using the definition of v + = x + s + μ and Eq. ( 10), we have Moreover, the third equation of system (11) yields Using ( 18), ( 19) and the fact that x 2 ≤ x 2 , we have Thus, using (20) and Lemmas 1 and 2, we obtain Substitution of this bound into (17) yields the desired inequality.

The corrector step
In this subsection, we deal with the corrector step.One can observe that the algorithm presented in Fig. 1 performs a full-Newton step as a corrector step, which can be obtained in the same way as the one presented in the primal-dual algorithm of Darvay et al. (2016).Thus, for the analysis of this case the lemmas proved in Darvay et al. ( 2016) can be applied.The next lemma gives a condition for strict feasibility of full-Newton step.
In the next lemma we show local quadratic convergence of the Newton process.
In the next subsection, we analyse the update of the duality gap after a main iteration of the algorithm.

The effect on duality gap after a main iteration
The purpose of the algorithm is to reduce the produced duality gap of the primal-dual pair.Therefore, in order to measure this reduction, in the next lemma we give an upper bound for the gap obtained after performing an iteration of the algorithm.

Determining appropriate values of parameters
In this section, we want to fix the parameters τ and θ with suitable values, which guarantee that after a main iteration, the proximity measure will not exceed τ .
Let (x, y, s) ∈ N (τ ) be the iterate at the start of a main iteration with x > 0 and s > 0 such that δ = δ(x, s; μ) ≤ τ < 1 2 .After a corrector step, by Lemma 6, we have It is obvious that the right-hand side of the above inequality is monotonically increasing with respect to δ, and this implies that Following the predictor step and the μ-update, by Lemma 4, we have where δ + and h(δ + , θ, n) are defined as in Lemma 3. It can be easily verified that h(δ + , θ, n) is decreasing with respect to δ + , so h(δ Using ( 22), ( 23), ( 24) and the fact that ρ is increasing with respect to δ + , we obtain It should be mentioned, that the iterates after the corrector steps are in the N (ω(τ )) neighbourhood, while the iterates after the predictor steps are in the N (τ ) neighbourhood.

Iteration bound
The next lemma gives an upper bound for the number of iterations produced by the algorithm.
Using the above result follows the main result of our paper.
Theorem 1 Let τ = 1 4 and θ = 1 5 √ n .Then, the corrector-predictor interior point algorithm given in Fig. 1 is well defined and the algorithm requires at most O 5 √ n log 5(x 0 ) T s 4 iterations.The output is a strictly feasible primal-dual solution (x, y, s) satisfying x T s ≤ .
It is worth mentioning that this CP algorithm has advantages compared to the onestep method presented in Darvay et al. (2016), because in the case of the one-step algorithm θ = 1 27 √ n , which is smaller than the value of θ used by us.Hence, this paper leads to a slightly better complexity result.

Numerical results
In order to demonstrate the efficiency of our CP algorithm we implemented it in the C++ programming language (Darvay and Takó 2012) in such a way to be able to compare it with the primal-dual (PD) method proposed in Darvay et al. (2016).To do so, we made some changes in the algorithm.In case of both algorithms, we used the normalized duality gap x T s n to obtain the value of the barrier parameter μ in the next iterate.In the implementation of the PD algorithm, we multiplied the normalized duality gap by σ = 0.95 in order to be able to reduce the value of μ in each iterate and to maintain the short-step strategy of the method.As it is usual in the case of the implementation of predictor-corrector algorithms, for our CP variant we applied Mehrotra's (1992) heuristics to get the value of μ for the corrector step.
Our CP algorithm performed one corrector and one predictor step in each iteration.After calculating the search direction, we determined the maximal step size to the boundary of the feasible region.We reduced the obtained step by multiplying it by the parameter ρ = 0.5.For both algorithms, we set the value of the proximity parameter to 10 −5 .We tested both algorithms on two set of problems.The first set contains the classical one-step method.We concluded that our algorithm is more efficient than the classical one.Future theoretical research direction includes extending our CP IPA to more general problems.The computational direction for further research includes implementation with different functions ψ used for AET.We would like to compare the effect of these functions on the computational results of the algorithms on quite large and well selected LO test problems.Moreover, this kind of computational study might give some insight into the good selection strategy of the target point of the search direction for different LO problems.Ideal situation would be, if the algorithm could be able to choose automatically among different ϕ functions depending on the structure of the problem.Furthermore, such computational study on the effect of different AETs could be based on the similar LO test set, like in the study of Illés and Nagy (2014) for pivot algorithms.