Infeasible Interior-Point Methods for Linear Optimization Based on Large Neighborhood

In this paper, we design a class of infeasible interior-point methods for linear optimization based on large neighborhood. The algorithm is inspired by a full-Newton step infeasible algorithm with a linear convergence rate in problem dimension that was recently proposed by the second author. Unfortunately, despite its good numerical behavior, the theoretical convergence rate of our algorithm is worse up to square root of problem dimension.


Introduction
We consider the linear optimization (LO) problem in standard format, given by (P), and its dual problem given by (D) (see Sect. 3). We assume that the problems (P) and (D) are feasible. As is well known, this implies that both problems have optimal solutions and the same optimal value. In Sect. 8, we discuss how infeasibility and/or unboundedness of (P) and (D) can be detected. Interior-point methods (IPMs) for LO are divided into two classes: feasible interiorpoint methods (FIPMs) and infeasible interior-point methods (IIPMs). FIPMs assume that a primal-dual strictly feasible point is available from which the algorithm can immediately start. In order to get such a solution several initialization methods have been presented, e.g., by Megiddo [1] and Anstreicher [2]. The so-called big M method of Megiddo has not been received well due to the numerical instabilities caused by the use of huge coefficients. A disadvantage of the so-called combined phase-I and phase-II method of Anstreicher is that the feasible region must be nonempty and a lower bound on the optimal objective value must be known.
IIPMs have the advantage that they can start with an arbitrary point and try to achieve feasibility and optimality, simultaneously. It is worth to mention that the bestknown iteration bound for IIPMs is due to Mizuno [3]. We refer to, e.g., [4] for a survey of IIPMs.
Recently, the second author proposed an IIPM with full-Newton steps [4]. Each main iteration of the algorithm improves optimality and feasibility simultaneously, by a constant factor until a solution with a given accuracy is obtained. The amount of reduction is too small for practical purposes. The aim of this paper is to design a more aggressive variant of the algorithm which improves optimality and feasibility faster. We emphasize that this is our aim and also what happens in practice. As we will see, however, our algorithm suffers the same irony that occurs for FIPMs, namely that the theoretical convergence rate of large-update methods is worse than for full-Newton methods.
In the analysis of our algorithm, we use so-called kernel function-based barrier function. Any such barrier function is based on a univariate function, named its kernel function. 1 Such functions can be found in [5] and are closely related to the so-called self-regular functions introduced in [6]. In these references only FIPMs are considered, and it is shown that these functions are much more efficient for the process of recentering, which is a crucial part in every FIPM, especially when an iterate is far from the central path. Not surprising, it turned out that these functions are also useful in our algorithm, where re-centering is also a crucial ingredient.
The paper is organized as follows. In Sect. 2, we describe the notations which we use throughout the paper. In Sect. 3, we explain the properties of a kernel function-based barrier function. In Sect. 4, as a preparation to our large neighborhood-based IIPM, we briefly recall the use of kernel-based barrier functions in large-update FIPMs, as presented in [5]. It will become clear in this section that the convergence rate highly depends on the underlying kernel function.
In Sect. 5, we describe our algorithm in detail. In our description we use a search direction based on a general kernel function. The algorithm uses two types of damped Newton steps: a so-called feasibility step and some centering steps. The feasibility step serves to reduce the infeasibility, whereas the centering steps keep the infeasibility fixed, but improve the optimality. This procedure is repeated until a solution with enough feasibility and optimality is obtained. Though many parts of our analysis are valid for general kernel function, at some places we restrict ourselves to a specific kernel function. We show that the complexity of our algorithm based on this kernel function is slightly worse than the complexity of the algorithm introduced and discussed by Salahi et al. [7], who used another variant of this function. Note that the best-known iteration bound for the IIPMs that use a large neighborhood of the central path is due to Potra and Stoer [8] for a class of superlinearly convergent IIPMs for sufficient linear complementarity problems (LCP).

Notations
We use the following notational conventions. If x, s ∈ R n , then xs denotes the coordinatewise (or Hadamard) product of the vectors x and s. The nonnegative orthant and positive orthant are denoted as R n + and R n ++ , respectively. If z ∈ R n + and f : for some positive constant c. For z ∈ R n , we denote the l 1 -norm by z 1 , and the Euclidean norm by z .

Barrier Functions Based on Kernel Functions
Let the primal LO problem be defined as and its dual problem as Here, A is a full row rank matrix in R m×n , b, y ∈ R m , c, x, s ∈ R n ; x, y and s are vectors of variables. It is well known that the triple (x, y, s) is optimal for (P) and (D) if and only if The first two equations require primal and dual feasibility, respectively, and xs = 0 is the so-called complementarity condition. Note that since x and s are nonnegative, xs = 0 holds if and only if x T s = 0. Therefore, since the feasibility conditions imply x T s = c T x −b T y, the complementarity condition resembles the fact that at optimality the primal and dual objective values coincide.
Interior-point methods (IPMs) replace the complementarity condition by xs = μe, where μ > 0. This gives rise to the following system: where xs = μe is called the centering condition. It has been established (see, e.g., [9,Theorem II.4]) that the system (2) has a solution for some μ > 0 if and only if it has a solution for every μ > 0. Moreover, this holds if and only if (P) and (D) satisfy the interior-point condition (IPC), i.e., there exist a feasible x > 0 for (P) and a feasible (y, s) with s > 0 for (D). If the matrix A has full row rank, then this solution is unique. It is denoted by (x(μ), y(μ), s(μ)) and called the μ-center of (P) and (D). The μ-centers form a virtual path inside the feasibility region leading to an optimal solution of (P) and (D). This path is called the central path of (P) and (D). When driving μ to zero, the central path converges to an optimal triple (x, y, s) for (P) and (D). We proceed by showing that the μ-center can be characterized as the minimizer of a suitably chosen primal-dual barrier function. In fact we will define a wide class of such barrier functions, each of which is determined by a kernel function. A kernel function is just a univariate nonnegative function ψ(t), where t > 0, which is strictly convex, minimal at t = 1 and such that ψ(1) = 0, whereas ψ(t) goes to infinity both when t goes to zero and when t goes to infinity. Now let (x, y, s) be a feasible triple with x > 0 and s > 0. We call such a triple strictly feasible. We define the variance vector for this triple with respect to μ as follows: v := xs μ .
Observe that v = e holds if and only if (x, y, s) is the μ-center of (P) and (D). Given any kernel function ψ we extend its definition to R n ++ according to It is obvious that Ψ (v) is nonnegative everywhere, and Ψ (e) = 0. Yet we can define a barrier function Φ(x, s, μ) as follows: where v is the variance vector of (x, y, s) with respect to μ. It is now obvious that Φ(x, s, μ) is well-defined, nonnegative for every strictly feasible triple, and moreover, This implies that (x(μ), y(μ), s(μ)) is the (unique) minimizer of Φ(x, s, μ). As in [5], we call the kernel function ψ eligible iff it satisfies the following technical conditions.
In the sequel it is always assumed that ψ is an eligible kernel function. Properties of eligible kernel functions will be recalled from [5] without repeating their proofs. We would like to mention that a kernel function ψ, introduced and discussed in [5], has the following format The term t 2 −1 2 is called the growth term which goes to infinity as t goes to infinity. ψ b has the property that it goes to infinity as t tends to zero. This term is called the barrier term of ψ.

A Class of Large-update FIPMs for LO
In this section we recall from [5] some results for a large-update FIPM for solving (P) and (D) using a kernel function-based barrier function.
Let us start with a definition. We call a triple (x, y, s), with x > 0, s > 0, an ε-solution of (P) and (D) iff the duality gap x T s and the norms of the residual vectors r b := b − Ax, and r c := c − A T y − s, do not exceed ε. In other words, defining (x, y, s) := max{x T s, r b , r c }, we say that (x, y, s) is an ε-solution if (x, y, s) ≤ ε. We assume, without any loss of generality, that the triple is a primal-dual feasible solution. 2 We then have x 0 s 0 = μ 0 e for μ 0 = 1. This means that (x 0 , y 0 , s 0 ) is the 1-center, and hence Φ(x 0 , s 0 , μ 0 ) = 0. We use this triple to initialize our algorithm. Each main (or outer) iteration of the algorithm starts with a strictly feasible triple (x, y, s) that satisfies Φ(x, s, μ) ≤ τ for some μ ∈]0, 1], where τ is a fixed positive constant. It then constructs a new triple (x + , y + , s + ) such that Φ(x + , s + , μ + ) ≤ τ with μ + < μ. When taking τ small enough, we obtain in this way a sequence of strictly feasible triples that belong to small neighborhoods of a sequence of μ-centers, for a decreasing sequence of μ's. As a consequence, the sequence of constructed triples (x, y, s) converges to an optimal solution of (P) and (D).
It is easy to compute the number the outer iterations. Using μ 0 = 1 one has the following result.
outer iterations to generate a strictly feasible ε-solution.
The main task is therefore to get a sharp upper estimate for the number of inner iterations during an outer iteration. We now describe how such an estimate is obtained. We go into some detail, though without repeating proofs, because the results that we recall below are relevant for the IIPM that we discuss in the next section. As said before, at the start of each outer iteration we have a strictly feasible triple (x, y, s) and μ > 0 such that Φ(x, s, μ) ≤ τ . We first need to estimate the increase in Φ when μ is updated to μ + = (1 − θ)μ. For this we need the following Lemma.
Now let v be the variance vector of (x, y, s) with respect to μ. Then one easily understands that the variance vector v + of (x, y, s) with respect to μ + is given by where the last inequality holds because ρ is monotonically increasing and Ψ (v) := Φ(x, s, μ) ≤ τ . Hence the numberτ , defined bȳ is an upper bound for the value of Ψ after a μ-update. Note that this bound depends not on the triple (x, y, s), but only on the kernel function ψ and the parameters n, τ and θ .
To simplify the notation we redefine μ according to μ := μ + . Thus we need to deal with the following question: given a triple (x, y, s) such that Φ(x, s, μ) ≤τ , how many inner iterations are needed to generate a triple (x, y, s) such that Φ(x, s, μ) ≤ τ . To answer this question we have to describe an inner iteration. It has been argued in Section 2.2. of [5] that it is natural to define search directions (Δx, Δy, Δs) by the system This system has a unique solution. It may be worth pointing out that if ψ(t) = ψ 1 (t), then −μv∇Ψ (v) = μe − xs, and hence the resulting direction is the primal-dual Newton direction that is used in all primal-dual FIPMs. By doing a line search in this direction with respect to Ψ we get new iterates where α is the step size. According to [5,Lemma 4.4], we use below the following default step size: where ρ is the inverse function of − 1 2 ψ (t), and The closeness of (x, y, s) to the μ-center is measured by Ψ (v), where v is the variance vector of (x, y, s) with respect to the current value of μ. The initial triple (x, y, s) is as given by (5) and μ = 1. So we then have Ψ (v) = 0 ≤ τ . After a

Input:
a threshold parameter τ > 0; an accuracy parameter ε > 0; a fixed barrier update parameter θ , 0 < θ < 1; begin x := e; y := 0; s := e; μ := 1; Then a sequence of inner iterations is performed to restore the inequality Ψ (v) ≤ τ . Then μ is updated again and so on. This process is repeated until nμ falls below the accuracy parameter ε after which we have obtained an ε-solution.
To estimate the number of inner iterations we proceed as follows. Denoting the decrease in the value of Ψ as ΔΨ , it was shown in [5,Theorem 4.6] that Since the kernel function ψ is eligible, the last expression is increasing in δ(v) [5,Lemma 4.7]. Besides, by [5,Theorem 4.9], δ(v) is bounded from below as follows: Combining (7) and (8), we arrive at Following [5], let γ be the smallest number such that for some positive constant κ, whenever Ψ (v) ≥ τ . In [10], it is established that such constants κ and γ exist for each kernel functions. When denoting the value of the barrier function after the μ-update as Ψ 0 and the value after the k-th inner iteration as Ψ k , it follows from (9) and (10) that withτ as in (6). At this stage we may point out why the use of kernel functions other than the logarithmic kernel function may be advantageous. The reason is that if ψ is the logarithmic kernel function, then γ = 1, whence we obtain This resembles the well-known fact that the best lower bound for the decrease in the logarithmic barrier function is a fixed constant, no matter what the value of Ψ (v) is. As we will see, smaller values of γ can be obtained for other kernel functions, which leads to larger reductions in the barrier function value, and hence lower iteration numbers. By [5, Lemma 5.1], (11) implies that the number of inner iterations will not exceed Multiplying this number by the number of outer iterations, as given by Lemma 4.1, we obtain the following upper bound for the total number of iterations: Given a kernel function ψ, it is now straightforward to compute the resulting iteration bound from this expression.
This function is introduced and discussed in [11] whereby it is established that γ = q+1 2q and κ = 1 124 q . As a result, the iteration number turns out to be The expression qn q+1 2q is minimal at q = 1 2 log n and then it is equal to 1 2 e log n √ n. Hence we obtain the iteration bound for the algorithm, which is the best-known iteration bound for large-update FIPMs based on kernel functions. It should be mentioned that the best-known iteration bound for the FIPMs that use a large neighborhood of the central path is due to Potra [12] who designed a superlinearly convergent predictor-corrector algorithm for linear complementarity problems that has an O( √ nL) iteration bound, with L denoting the length of a binary string encoding the problem data.
In this paper, we consider a IIPM based on the use of a kernel function. Although many of the results below hold for any eligible kernel function, we will concentrate of the kernel ψ q .

A Class of Large-Update IIPMs for LO
As usual in the theoretical analysis of IIPMs, we take the initial iterates as follows: where e denotes the all-one vector in R n , and ζ is a number such that for some optimal solutions (x * , y * , s * ) of (P) and (D). It is worth noting that if the data A, b, c are integral and L denotes the length of a binary string encoding the triple (A, b, c), then it is well known that there exist optimal solutions x * of (P) and (y * , s * ) of (D) that satisfy (14) with ζ = 2 L [4, Section 4.7]. Following [4], for any 0 < ν ≤ 1 we consider the primal-dual perturbed pair (P ν ) and (D ν ), defined as follows: where r 0 b and r 0 c denote the initial primal and dual residual vectors, respectively. Note that, if ν = 0, then the perturbed pair (P ν ) and (D ν ) coincides with the original pair (P) and (D). Moreover, (x 0 , y 0 , s 0 ) is strictly feasible for the perturbed problems (P ν ) and (D ν ) if ν = 1. In other words, if ν = 1, then (P ν ) and (D ν ) satisfy the IPC and since x 0 s 0 = μ 0 e, with μ 0 = ζ 2 , the μ 0 -center of (P 1 ) and (D 1 ) is given by x 0 , y 0 , s 0 . The following Lemma is crucial. Hence, for any ν ∈]0, 1], the problems (P ν ) and (D ν ) satisfy the IPC. This implies that the central path of (P ν ) and (D ν ) exists. The μ-center of (P ν ) and (D ν ) is uniquely determined by the system In the sequel, the parameters ν and μ will always satisfy μ = νζ 2 . Therefore, we feel free to denote μ-center of (P ν ) and (D ν ) as (x(ν), y(ν), s(ν)). The set of these μ-centers, as ν runs through the interval ]0, 1], is called the homotopy path of (P) and (D) with respect to ζ . By driving ν to zero, the homotopy path converges to an optimal solution of (P) and (D) [14]. Our algorithm starts at the μ-center for ν = 1 and then follows this homotopy path to obtain an ε-solution of (P) and (D). Note that along the homotopy path the residual vectors are given by νr 0 b and νr 0 c , whereas the duality gap by νnζ 2 . So these quantities converge to zero with the same speed. As a consequence we have (x(ν), y(ν), s(ν)) = ν (x(1), y(1), s(1)) = ν (ζ e, 0, ζ e) , ν ∈]0, 1].

An Outer Iteration of the Algorithm
Just as in the case of large-update FIPMs we use a primal-dual barrier function to measure closeness of the iterates to the μ-center of (P ν ) and (D ν ). So, just as in Sect. 4, Ψ (v) will denote the barrier function based on the kernel function ψ(t), as given in (3). Here, v denotes the variance vector of a triple (x, y, s) with respect to μ > 0, and we define Φ(x, s, μ) as in (4). The algorithm is designed in such a way that at the start of each outer iteration we have Ψ (v) ≤ τ for some threshold value τ = O(1). As Ψ (v) = 0 at the starting points (13), the condition Ψ (v) ≤ τ is certainly satisfied at the start of the first outer iteration.
Each outer iteration of the algorithm consists of a feasibility step and some centering steps. At the start of the outer iteration we have a triple (x, y, s) that is strictly feasible for (P ν ) and (D ν ), for some ν ∈]0, 1], and that belongs to the τ -neighborhood of the μ-center of (P ν ) and (D ν ), where μ = νζ 2 . We first perform a feasibility step during which we generate a triple x f , y f , s f which is strictly feasible for the perturbed problems (P ν + ) and (D ν + ), where ν + = (1 − θ)ν with θ ∈]0, 1[, and moreover, close enough to the μ + -center of (P ν + ) and (D ν After the feasibility step, we perform some centering steps to get a strictly feasible triple x + , y + , s + of (P ν + ) and (D ν + ) in the τ -neighborhood of the μ + -center of (P ν + ) and During the centering steps the iterates stay feasible for (P ν + ) and (D ν + ). Hence for the analysis of the centering steps we can use the analysis for FIPMs, presented in Sect. 4. From this analysis we derive that the number of centering steps will not exceed Φ x f , s f , μ + γ κγ , where the parameters γ and κ depend on the kernel function ψ. Hence we are left with the problem of defining a suitable search direction Δ f x, Δ f y, Δ f s for the feasibility step and to determine θ such that after the feasibility step we have Φ x f , s f , μ + ≤ τ f for some suitable value of τ f . The number of outer iterations will be 1 θ log (ζ e,0,ζ e) ε . Thus, the total number of iterations will not exceed

Feasibility Step
For the search direction in the feasibility step, we use the triple Δ f x, Δ f y, Δ f s that is (uniquely) defined by the following system: Then, defining the new iterates according to In the same way one shows that c − A T y f − s f = ν + r 0 c . Hence it remains to find θ such that x f and s f are positive and Φ(x f , s f , μ + ) ≤ τ f . This is the hard part of the analysis of our algorithm, which we leave to the subsection below. We conclude this part with a formal description of the algorithm, in Fig. 2, and a graphical illustration in Fig. 3.
The straight lines in Fig. 3 depict the central paths of the pair (P ν ) and (D ν ) and the pair (P ν + ) and (D ν + ). The τ -neighborhoods of the μand μ + -centers are shown by the dark gray circles. The light gray region specifies the τ f -neighborhood of the μ + -center of (P ν + ) and (D ν + ). The feasibility step is depicted by the first arrow at the right-hand side. The rest of the arrows stand for the centering steps. The iterates are shown by the circlets.

Analysis of the Feasibility Step
The feasibility step starts with some strictly feasible triple (x, y, s) for (P ν ) and (D ν ) and μ = νζ 2 such that Our goal is to find θ such that after the feasibility step, with step size θ , the iterates Ψ To proceed, we define scaled feasibility directions d Using these notations, we may write This shows that Therefore, d Yet we turn to the requirement that Ψ (v f ) ≤ τ f . Using (18), (19) and xs = μv 2 , we write Hence, if θ < θ max then we may write Recall that ρ(s) ≥ 1 for all s ≥ 0 and ρ(s) is monotonically increasing. Furthermore, ψ(t) is monotonically increasing for t ≥ 1. Hence we obtain , the coefficient of n in the above upper bound for Ψ (v f ) does not depend on n. Hence the Lemma follows.
Due to Lemma 5.2 it suffices for our goal to find θ such that Ψ (v) ≤τ wherê τ = O(n). In the sequel, we consider Ψ (v) as a function of θ , denoted as f 1 (θ ). So we have We proceed by deriving a tight upper bound for f 1 (θ ), thereby using similar arguments as in [5]. Since the kernel function ψ(t) is eligible, Ψ (v) is e-convex (cf. [5, Lemma 2.1]), whence we have The first and the second derivatives of f (θ ) are as follows: Since ψ (t) < 0, ∀t > 0, it follows that ψ (t) is monotonically decreasing. From this we deduce that where v min := min(v) and θ small enough, i.e., such that v min − θ d f x > 0. Substitution into (21) gives By integrating both sides of this inequality with respect to θ , while using that f (0) = 0, as follows from (20), we obtain Integrating once more, we get where the last inequality holds because ψ is convex. 3 The first derivative with respect to v min of the right-hand side expression in this inequality is given by Since ψ is (strictly) decreasing, this derivative is negative. Hence it follows that the expression is decreasing in v min . Therefore, when θ and d f x are fixed, the less the v min is, the larger the expression will be. Below we establish how small v min can be when δ(v) is given.
Since the inverse function ρ of −ψ /2 is monotonically decreasing, this is equivalent to Hence the smallest possible value of v min is ρ(δ(v)), and this value is attained in the (exceptional) case where v min is the only coordinate of the vector v that differs from 1. So we may assume that v min = ρ(δ(v)). This implies −ψ (v min )/2 = δ(v) and hence ψ (v min ) ≤ 0, whence v min ≤ 1.
In the sequel, we denote δ(v) simply as δ. Substitution into (22) gives that Hence we certainly have f (θ ) ≤τ if Since f (0) = Ψ (v) ≤ τ , it can be verified that the last inequality holds if Since ρ is decreasing, this is equivalent to Note that if θ approaches zero, then the left-hand side expression converges to ρ(δ) and the right-hand side expression to zero. The left-hand side is decreasing in θ , whereas the right-hand side is increasing. The largest possible θ makes both sides equal. In order to get a tight approximation for this value, we first need to estimate d f x . This is done in the next section, where we follow a similar approach as in [15].

Upper Bound for d f x
One may easily check that the system (16) At this stage we use that the initial iterates are given by (13) and (14), so we have where ζ > 0 is such that for some optimal solutions x * of (P) and (y * , s * ) of (D). Using a similar argument as given in [15,Section 4.3], it can be proven that In [15,Lemma 4.3], it is shown that Using xs = μv 2 and (23) we may write Substitution in (31) yields that Substitution in (30), also using that v min ≥ ρ(δ), we obtain

Condition for θ
Yet we return to the condition (24) on θ : The right-hand side is increasing in d f x . Therefore, due to (32), it suffices if Obviously, this implies that θ n(1+ρ(−δ) 2 ) < ρ(δ) 2 . Therefore, there exists α ∈]0, 1[, independent of n, such that Clearly, using this θ , the expression (33) can be restated as Our objective is to find the largest possible α satisfying this inequality. In order to proceed we need bounds for δ = δ(v) and ρ(δ). The next section deals with some technical lemmas which are useful in deriving these bounds.

Some Technical Lemmas
Recall that ρ is defined as the inverse function of − 1 2 ψ (t) and ρ as the inverse function of ψ(t) for t ≥ 1. We also need the inverse function of ψ(t) for t ∈]0, 1], which we denote as χ . To get tight estimates for these inverse functions, we need the inverse functionsχ andρ of ψ b (t) and −ψ b (t), respectively. Recall that, as it was shown in [5], one has Hence ψ b (t) is monotonically decreasing and ψ b (t) is monotonically increasing. This implies that ψ b (t) and −ψ b (t) have inverse functions and these function are monotonically decreasing. We denote these inverse functions asχ andρ, respectively.
Recall that ρ is the inverse function of ψ(t) for t ≥ 1. The following two results are less trivial than the preceding two lemmas.

Lemma 6.4 One has, for any
Proof The left-hand side inequality in the lemma is due to [5,Theorem 4.9]. The proof of the right-hand side inequality can be obtained by slightly changing the proof of [5,Theorem 4.9] and is therefore omitted.

Iteration Bound of the IIPM Based on ψ q
Since the kernel function ψ q led to the best-known iteration bound for large-update FIPMs based on kernel functions, from now on we restrict ourselves to the case where ψ = ψ q . Thus, we have In the current case the inverse functionsχ andρ are defined as follows: Since χ and − 1 2 ψ are decreasing, − 1 2 ψ χ is increasing. Hence, letting Ψ (v) ≤ τ , Lemma 6.4 implies that By Lemma 6.1 we have χ(τ ) ≥χ τ + 1 2 . Using once more that − 1 2 ψ is decreasing we obtain 2δ ≤ −ψ χ τ + 1 2 .

Complexity Analysis
As we established in (15), the total number of iterations is at most We assume thatτ = O(n). Due to Lemma 5.2 we then also have τ f = O(n), provided that 1/ √ 1 − θ = O(1). Due to (42) the latter condition is satisfied. To simplify the presentation we use τ f = n in the analysis below, but our results can easily be adapted to the case whereτ = O(n). When choosing θ maximal such that (42) holds, we have Hence, substitution of γ = q+1 2q and κ = 1 124 q yields that the total number of iterations does not exceed The expression q 3 n 1 2q is minimal if q = log n/6, and then it is equal to e 4 (log n) 3 /512. This value of q satisfies (41), since log(2e) ≤ 6. Hence we obtain the following iteration bound: It is worth mentioning that if ψ 1 were used, then the iteration bound of our algorithm would be O n 2 log nζ 2 ε , which is a factor √ n/(log n) 3 worse than (44).

Detecting Infeasibility or Unboundedness
As indicated in the introduction, the algorithm will detect infeasibility or/and unboundedness of (P) and (D) if no optimal solutions exist. In that case Lemma 5.1 implies the existence ofν > 0 such that the perturbed pair (P ν ) and (D ν ) satisfy the IPC if and only if ν ∈]ν, 1]. As long as ν + = (1 − θ)ν >ν the algorithm will run as it should, with θ given by (42). However, if ε is small enough, at some stage it will happen that ν >ν ≥ ν + . At this stage the new perturbed pair does not satisfy the IPC. This will reveal itself since at that time we necessarily have θ max <θ . If this happens we may conclude that there is no optimal pair (x * , s * ) satisfying x * + s * ∞ ≤ ζ . We refer to [4] for a more detailed discussion.

Computational Results
In this section, we present the numerical results for a practical variant of the algorithm described in Sect. 5. Theoretically, the barrier parameter μ is updated by a factor (1−θ) with θ given by (42), and the iterates are kept very close to the μ-centers, namely the τ -neighborhood of the μ-centers, with τ = 1 8 . In practice, it is not efficient to do in this way and not even necessary either. We present a variant of the algorithm which uses a predictor-corrector step in the feasibility step. Moreover, for the parameter τ , defined in Sect. 5.1, we allow some larger value than 1 8 , e.g., τ = O(n). We set τ =τ = O(n) withτ defined as in Sect. 5.1. As a consequence, the algorithm does not need centering steps. We chooseτ according to the following heuristics: if n ≤ 500, thenτ = 100n, for 500 ≤ n ≤ 5000, we chooseτ = 10n and for n ≥ 5000, we set τ = 3n. We compare the performance of the algorithm with the well-known LIPSOL package [16].

Starting Point and Stopping Criterion
A critical issue when implementing a primal-dual method is to find a suitable starting point. It seems sensible to obtain a starting point which is well centered and as close to a feasible primal-dual point as possible. The one suggested by theory, i.e., given by (13), being nicely centered, it may be quite far from the feasibility region. Moreover, to find a suitable ζ is another issue.
In our implementation, we use a starting point x 0 , y 0 , s 0 which is proposed by Lustig et al. [17] and inspired by the starting point used by Mehrtora [18]. Since we are interested in a point which is in the τ -neighborhood of the μ 0 -center. If Φ(x 0 ; s 0 ; μ 0 ) > τ, we increase μ 0 until Φ(x0; s0; μ 0 ) ≤ τ .
As in LIPSOL, our algorithm terminates if, for ε = 10 −6 , either or |x T s − (x + ) T s + | < ε occurs. The condition (45) measures the total relative errors in the optimality conditions (1), while the later criterion terminates the program if only a tiny improvement is obtained on the optimality. In fact, it prevents the program from stalling. We include this criterion following Lustig [19].

Feasibility Step Size
As in other efficient numerical experiments, e.g., [16,17], regardless of the theoretical result, we apply different step sizes along the primal step Δx and the dual step (Δy, Δs). This implies that the feasibility improves much faster than when identical step sizes are used. Letting (x, y, s) be the current iterates and (Δx, Δy, Δs) the Newton step, we obtain the maximum step sizes θ p max and θ d max in, respectively, the primal and the dual spaces as follows: The goal is to keep the iterates close to the μ-center, i.e., in itsτ -neighborhood wherê τ is defined in Sect. 5.3. Thus, lettingθ be such that Φ(x +θθ p max Δx, s +θθ d max Δs, μ) ≤τ , the primal and the dual step sizes θ p and θ d are defined as follows: θ p =θθ p max and θ d =θθ d max .

Solving the Linear System
We apply the backslash command of MATLAB ('\') to solve the normal equations whereb is some right-hand side and D, in LO case, has the form Denoting M := AD 2 A T in (46), whenever the matrix M is ill-conditioned, we could obtain some more accurate solution by perturbing M as where I is the identity matrix with size of M.

A Rule for Updating μ
Motivated by the numerical results, and considering the fact that the Mehrotra's PC method has become the most efficient in practice and used in most IPM-based software packages, e.g., [16,[20][21][22], we present the numerical results of the variant of our algorithm which uses Mehrotra's PC direction at the feasibility step. At the feasibility step, we apply the system to obtain the affine-scaling directions (Δ a x , Δ a y , Δ a s ). Then, the maximum step sizes θ p max and θ d max in, respectively, primal and dual spaces are calculated as described in Sect. 9.2. Then, defining We useσ = 0.3, as the default value ofσ . If σ < 1, we calculate the new barrier update parameter μ as follows: Then, if necessary, by increasing μ new by a constant factor, say 1.1, we derive some μ new for which Note that μ new is acceptable only if μ new < μ.
In the case where μ new ≥ μ or σ ≥ 1, we do not perform any μ-update and proceed with the following corrector step which attempts to center toward the current μ-center: Recall that if σ ≥ 1, then we have an increase in the duality gap and hence the use of μ new is no longer sensible. If σ < 1 and μ new < μ, we use the system (47) with μ = 0 as a corrector step.
The feasibility step Δ f x, Δ f y, Δ f s is obtained as follows: Next, we calculate the maximum step sizes θ p max and θ d max in, respectively, the primal and the dual spaces along the step Δ f x, Δ f y, Δ f s , as described in Sect. 9.2. Assuming that μ is the outcome of the predictor step, we obtain a θ such that withτ chosen as described in Sect. 9.4. By setting θ p = θθ d max and θ d = θθ d max , we calculate the new iterates x f , y f , s f as follows:

Results
In this section, we present our numerical results. Motivated by the theoretical result which says that the kernel function ψ q gives rise to the best-known theoretical iteration bound for large-update IIPMs based on kernel functions, we compare the performance of the algorithm described in the previous subsection based on both the logarithmic barrier function and the ψ q -based barrier function. As the theory suggests, we use q = log n 6 in ψ q . Our test was done on a standard PC with Intel Core TM 2 Duo CPU and 3.25 GB of RAM. The code was implemented by version 7.11.0 (R2010b) of MATLAB on Windows XP Professional operating system. The problems chosen for our test are from the NETLIB set. To simplify the study, we choose the problems which are in (P) format, i.e., there is no nonzero lower bound or finite upper bound on the decision variables.
Before solving a problem using the algorithm, we first solve it for singleton variables. This helps to reduce the size of the problem.
Numerical results are tabulated in Table 1. In second and the fourth columns, we listed the total number of iterations of the algorithm based on, respectively, ψ 1 , the kernel function of the logarithmic barrier function, and ψ q . The third and fifth columns stand for the quantity E(x, y, s). The iteration numbers of the LIPSOL package are given in the sixth column of these tables. Finally, in the seventh column, we listed the quantity E(x, y, s) of the LIPSOL package. In each row, the Italic number denotes the smallest of the iteration numbers of the three algorithms, namely our algorithm based on ψ 1 and ψ q , and LIPSOL, and the bold number denotes the smallest of the iteration numbers of ψ 1 -based and ψ q -based algorithms. As it can be noticed from the last row of the table, the overall performance of the algorithm based on ψ 1 is much better than that the variant based on ψ q . However, in some of the problems, the ψ q -based algorithm outperforms the ψ 1 -based algorithm. This amounts to the problems AGG, BANDM, DEGEN2, DEGEN3, SCSD1, SCSD6, SCSD8 and SHARE2B. The columns related to LIPSOL prove that it is still the best; however, our ψ 1 -based algorithm saved one iteration compared with LIPSOL, to solve AGG2 and AGG3, and two iterations when solving STOCFOR1.

Conclusions
In this paper, we analyze infeasible interior-point methods (IIPMs) for LO based on a large neighborhood. Our work is motivated by [4] in which Roos presents a full-Newton IIPM for LO. Since the analysis of our algorithms requires properties of barrier functions based on kernel functions that are used in large-update feasible interior-point methods (FIPMs), we present primal-dual large-update FIPMs for LO based on kernel functions, as well. In Roos' algorithm, the iterates move within small neighborhood of the μ-centers of the perturbed problem pairs. In many IIPMs, the algorithm reduces the infeasibility and the duality gap at the same rate. His algorithm has the advantages that it uses full-Newton steps and hence no calculation of step size is needed, and moreover its theoretical iteration bound is O(n log( (ζ e, 0, ζ e)/ε)) which coincides with the bestknown iteration bound for IIPMs. Nevertheless, it has the deficiency that it is too slow in practice as it reduces the quantity (x, y, s) by a factor (1 − θ) with θ = O( 1 n ) at an iteration.
We attempt to design a large-update version of Roos' algorithm which allows some larger amount of reduction of (x, y, s) at an iteration than 1 − θ with θ = O(1/n). This requires that the parameter θ is larger than O(1/n), even θ = O(1).
Unfortunately, the result of Sect. 5 implies that θ is O(1/(n(log n) 2 )) which yields O(n √ n(log n) 3 log( (ζ e, 0, ζ e)/ε) iteration bound for a variant. Since the theoretical result of the algorithm is disappointing, we rely on the numerical results to establish that our algorithm is a large-update method. A practically efficient version of the algorithm is presented and its numerical result is compared with the well-known LIP-SOL package. Fortunately, the numerical results seem promising as our algorithm has iteration numbers close to those of LIPSOL and, in a few cases, outperforms it. This means that IIPMs suffer from the same irony as FIPMs, i.e., regardless of their nice practical performance, the theoretical complexity of large-update methods is worse. Recall that the best-known iteration bound for large-update IIPMs is O(n √ n log n log( (ζ e, 0, ζ e)/ε) which is due to Salahi et al. [7]. As in LIPSOL [16], different step sizes in the primal and the dual spaces are used in our implementation. This gives rise to a faster achievement in feasibility than when identical step sizes are used. Moreover, inspired by the LIPSOL package, we use a predictor-corrector step in the feasibility step of the algorithm.
Future research may focus on the following areas: -As mentioned before, our algorithm has a factor (log n) 2 worse iteration bound than the best-known iteration bound for large-update IIPMs. One may consider how to modify the analysis such that the iteration bound of our algorithm is improved by a factor (log n) 2 . -As mentioned in Sect. 10, according to the analysis of our algorithm presented in Sect. 5, the barrier-updating parameter θ is O(1/(n(log n) 2 )). This yields the loose iteration bound given by (44). This slender value of θ is obtained because of some difficulties in the analysis of the algorithm which uses the largest value of θ , satisfying (24), to assure that Ψ (v) = O(n). This value of θ is much smaller than the best value we may choose. Assuming n = 60, the largest value of θ satisfying Ψ (v) = n is 0.849114 while the value of θ suggested by theory is 0.141766. A future research may focus on some new analysis of the algorithm which yields some larger value of θ . -Roos' full-Newton step IIPM was extended to semidefinite optimization (SDO) by Mansouri and Roos [23], to symmetric optimization (SO) by Gu et al. [24] and to LCP by Mansouri et al. [25]. An extension of large-update FIPMs based on kernel functions to SDO was presented by El Ghami [10]. One may consider how our algorithm behaves in theory and practice when it is extended to case of SDO, SO and LCP.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.