Polynomial worst-case iteration complexity of quasi-Newton primal-dual interior point algorithms for linear programming

Quasi-Newton methods are well known techniques for large-scale numerical optimization. They use an approximation of the Hessian in optimization problems or the Jacobian in system of nonlinear equations. In the Interior Point context, quasi-Newton algorithms compute low-rank updates of the matrix associated with the Newton systems, instead of computing it from scratch at every iteration. In this work, we show that a simplified quasi-Newton primal-dual interior point algorithm for linear programming enjoys polynomial worst-case iteration complexity. Feasible and infeasible cases of the algorithm are considered and the most common neighborhoods of the central path are analyzed. To the best of our knowledge, this is the first attempt to deliver polynomial worst-case iteration complexity bounds for these methods. Unsurprisingly, the worst-case complexity results obtained when quasi-Newton directions are used are worse than their counterparts when Newton directions are employed. However, quasi-Newton updates are very attractive for large-scale optimization problems where the cost of factorizing the matrices is much higher than the cost of solving linear systems.


Introduction
Let us consider the following general linear programming problem where x, c ∈ R n , b ∈ R m and A ∈ R m×n . We assume that (1) is feasible and the rows of A are linearly independent. Define function F : R 2n+m → R 2n+m by where X, Z ∈ R n×n are diagonal matrices defined by X = diag(x) and Z = diag(z), respectively, and e is the vector of ones of appropriate size. First order necessary optimality conditions for (1) state that, if x * ≥ 0 is a minimizer, then there exist z * ∈ R n , z * ≥ 0, and λ * ∈ R m such that F (x * , λ * , z * ) = 0 holds. Interior point methods (IPMs) try to follow the so-called central-path of problem (1), defined by the solution of the perturbed KKT system F (x, λ, z ) = 0 0 µe T , as µ → 0. Instead of solving such a system exactly, primal-dual IPMs apply one iteration of Newton method for a given value of µ k at iteration k. In order to calculate this step, the Jacobian of F is needed With an iterate (x k , λ k , z k ) at step k, the classical Newton direction is calculated by solving the following system

Background
Given a function G : R N → R N , suppose that we want to solve the nonlinear system G(x) = 0. Secant methods iteratively construct a linear model M k (x) of G which interpolates the last two computed iterates of the method. At each iteration, they need to compute an approximation to the Jacobian of G, which has to satisfy the secant equation where s k−1 =x k −x k−1 and y k−1 = G(x k ) − G(x k−1 ). There are infinitely many solutions to the secant equation for N ≥ 2 and different approaches generate different secant methods [16]. Among them, the Broyden "bad" approach uses the already computed approximation to the inverse of G atx k−1 , called H k−1 , to compute the current approximation H k as where and ρ k−1 = y T k−1 y k−1 . The Broyden "bad" update is a rank-1 update where H k is the matrix closest to H k−1 in the Frobenius norm which satisfies the secant equation. After ℓ updates of an approximation H k−ℓ , current approximation H k is given by For the specific case of this work, where G is given by F defined in (2), we have that N = 2n + m,x = (x, λ, z ) and the vectors s k−1 and y k−1 from the secant equation assume a more specific description whereᾱ k−1 ∈ (0, 1] is the step-size taken at iteration k − 1 towards the solution of (4). In [10], the authors described an interior point method based on low rank quasi-Newton approximations to the Jacobian of F . The Broyden updates were tested, and the computational experience revealed that the most efficient one was the Broyden "bad" update.
Since we are interested in finding an approximate solution of the linear system given by the Newton method (4), in the Broyden "bad" approach, given ℓ ≥ 0, the following direction is computed If H k−ℓ = J −1 k−ℓ and ℓ = 0, system (8) turns out to be exactly (4). Therefore, in the same way as discussed in [10], we assume that the initial approximation H k−ℓ is given by the perfect approximation J −1 k−ℓ . When ℓ > 0, the quasi-Newton procedure strongly uses the fact that the factorization of J k−ℓ (or a good preconditioner) has already been computed. In addition, Lemma 1, taken from [10], shows that, with this choice of initial approximation, the Broyden "bad" update has an interesting alternative interpretation. Lemma 1. Assume that k, ℓ ≥ 0 and H k is the approximation of J −1 k constructed by ℓ updates (6) using initial approximation H k−ℓ = J −1 k−ℓ . Given v ∈ R 2n+m , the computation of r = H k v is equivalent to the solution of Lemma 1 is the basis of the analysis developed in this work. It states that we can study quasi-Newton steps using the Jacobian of the Newton step. The only difference is the right-hand side. Using this property of the Broyden "bad" update we are able to extend the well known complexity results described in [21]. The difficulty in the analysis will be mostly caused by the extra term, added to the usual right-hand side of (4). It is important to note that Lemma 1 does not assume that the iterates are feasible, hence it is useful in both feasible and infeasible cases. Although by (6), the sparsity structure of the third row of B k (the inverse of H k ) is lost when ℓ ≥ 1, we can see that the structural sparsity of J k−ℓ can still be used to solve the linear systems.
Let us define a skeleton primal-dual quasi-Newton interior point algorithm. It is given by Algorithm 1 and generates a sequence of alternating Newton and quasi-Newton steps. Clearly, by the nature of update (5), the first step needs to be a Newton step.
Algorithm 1 Conceptual Quasi-Newton Interior Point algorithm. Input: Let us analyze what happens when a sequence of two steps is performed: at iteration k the Newton step is made (with stepsizeᾱ k ) and then at iteration k + 1 the quasi-Newton step is taken (with stepsizeᾱ). For Newton step at iteration k, we observe that Later, in the proofs of several technical results, we will need to analyze the error produced when the quasi-Newton direction (∆x k+1 , ∆λ k+1 , ∆z k+1 ) is multiplied by J k+1 : Applying Lemma 1 for iteration k + 1 with ℓ = 1 and then observing that (∆x k , ∆λ k , ∆z k ) solves the Newton system (4) and using (9), the third block equation in Lemma 1 gives Hence, using (7) and (10) where ∆Z k , ∆X k , ∆Z k+1 and ∆X k+1 are given by diag(∆z k ), diag(∆x k ), diag(∆z k+1 ) and diag(∆x k+1 ), respectively. Next we compute the new complementarity products obtained after a sequence of two steps, apply (11), and add and subtract the termᾱ 2 k ∆X k ∆Z k e to derive By multiplying both sides of equation (12) with e T we get the complementarity product at iteration k + 2: It is worth noting that the final expression in (12) involves a composite direction (ᾱ k ∆x k +ᾱ∆x k+1 ,ᾱ k ∆z k +ᾱ∆z k+1 ) which corresponds to an aggregate of two consecutive steps: in Newton direction at iteration k and in quasi-Newton direction at iteration k + 1. Much of the effort of the analysis presented in this paper is focused on this composite direction. Let us mention that we will also use the component-wise versions of equations (9), (11) and (12). For example, in case of (12) this gives Observe that equations (9)- (14) are valid for both feasible and infeasible algorithm. However, the analysis in Sections 3 and 4 will distinguish between these two cases because in the feasible one we are able to take advantage of the orthogonality of primal and dual directions and exploit it to deliver better final worst-case complexity results.
Before we conclude the brief background section and take the reader through a detailed analysis of different versions of the primal-dual quasi-Newton interior point algorithm, let us observe that equation (12) involves an important term γ 1 . By Lemma 1, γ 1 can be seen as the scalar coefficient of the projection of vector v onto the subspace generated by vector y k : Using the non-expansive property of projections we conclude that In the next lemma, a lower bound for y k is derived. It states that the denominator of (15) can be bounded away from zero if a sufficient decrease of µ = x T z/n is ensured and non-null step-sizes are taken. The bound for v involves the right-hand side in (8) and therefore depends on the feasibility of iterates and on the choice of the centering parameter σ. Lemma 2. Let k + 1 be a quasi-Newton iteration of Algorithm 1 and y k be the quasi-Newton vector defined by (7) to construct H k+1 by the Broyden "bad" update (6).
Proof. Ifᾱ k = 0 or ρ k = 0 the result trivially holds, so we can assume thatᾱ k , ρ k ∈ (0, 1]. Suppose, by contradiction, that y k < ρ kᾱk µ k /2. Therefore, by definition of y k , where the last result was obtained by adding up all the n previous inequalities and dividing by n. By hypothesis we have that µ k+1 ≤ (1 − ρ kᾱk )µ k and, therefore, which implies ρ k /2 > ρ k and is a clear absurd. Thus, we conclude that y k ≥ ρ kᾱk µ k /2.

Worst-case complexity in the feasible case
For all the results in this section, we suppose that (x 0 , λ 0 , z 0 ) is primal and dual feasible, given by Assumption 1.
Our analysis follows closely the theory in [21]. Under Assumption 1 the primal and dual directions are orthogonal to each other [21,Lemma 5.1]. We show that the same holds for quasi-Newton directions.
Proof. Using Lemma 1 with r = Hv defined by their respective terms in (8) and by the primal and dual feasibility of (x 0 , λ 0 , z 0 ), we observe that the first two block rows of system (8) (at iteration k + 1) given by are the same as in the system solved by the usual Newton step. Therefore, Corollary 1. If k is a Newton iteration and k + 1 is a quasi-Newton iteration, then ᾱ k ∆x k +ᾱ∆x k+1 T ᾱ k ∆z k +ᾱ∆z k+1 = ∆x k T ∆z k = ∆x k T ∆z k+1 = ∆z k T ∆x k+1 = 0.
Proof. Since, by equation (7), (∆x k , ∆λ k , ∆z k ) was computed by the Newton step at iteration k we have that ∆x k T ∆z k = 0 by the same arguments of Lemma 3. By Lemma 1 the first two equations do not change between iterations k and k + 1, and we have the desired results.
The second important result is that the quasi-Newton step can decrease the barrier parameter µ in exactly the same way as the Newton step (see [21,Lemma 5.1]). Recall that by definition µ = (x T z)/n.
Proof. By equation (13) and using Corollary 1, we obtain By dividing both sides of the last equation by n we get the desired result.
Then, using Lemma 4, after Newton step at iteration k and quasi-Newton step at iteration k + 1, we get It is worth noting that (in the feasible case) the term γ 1 originating from Lemma 1 does not have any influence on the value of µ(ᾱ).

The N 2 neighborhood
In this section we will consider a short-step interior point method and employ the notion of N 2 neighborhood of the central path where F is the set of primal and dual feasible points such that x, z > 0, see Assumption 1. For all considerations in this subsection we add the following assumption.
To further exploit (19) we need bounds on two terms which appear in it: 1 +ᾱγ 1 and the second-order error contributed by the composite direction ᾱ k ∆X k +ᾱ∆X k+1 ᾱ k ∆Z k +ᾱ∆Z k+1 e . The following technical result delivers a bound for |γ 1 |. (Observe that γ 1 is evaluated only when the quasi-Newton iteration is performed.) Lemma 5. Let k + 1 be a quasi-Newton iteration of Algorithm 1. Suppose that v in Lemma 1 is given by the right-hand side of (4) and Assumption 2 holds. Ifᾱ k ∈ (0, 1] and σ k ∈ [0, 1), then where γ 1 is defined in Lemma 1.
Proof. We use the assumptions of the lemma, the fact that the iterates are primal and dual feasible, the property e T (µ k+1 e − X k+1 Z k+1 e) = 0 and the relation By defining ρ k = 1 − σ k we can see that µ k+1 = (1 − (1 − σ k )ᾱ k )µ k = (1 − ρ kᾱk )µ k which ensures the sufficient decrease condition of Lemma 2. Sinceᾱ k > 0 and σ k < 1, by the assumptions of the lemma, we have that y k > 0. Then, by simple substitution of the previous equation and Lemma 2 in (15) we have Next we turn our attention to the error in the composite direction ᾱ k ∆X k +ᾱ∆X k+1 ᾱ k ∆Z k +ᾱ∆Z k+1 e and start from a technical result.
Proof. We use equations (4), (10) and (9) and some simple manipulations to obtain After adding and subtracting the term (ᾱ +ᾱ k (1 −ᾱ))µ k e we further rearrange the previous equation Then using µ k+1 = (1 −ᾱ k (1 − σ k ))µ k (which clearly holds for a step in Newton direction) and Lemma 4 which delivers a similar result for a step in quasi-Newton direction, we get: By substituting (23) in (22), the desired result is obtained.

Lemma 7. Let k + 1 be a quasi-Newton iteration of Algorithm 1 and suppose that Assumption 2 holds. Then
Proof. Let us first define D k = (X k ) 1/2 (Z k ) −1/2 and the scaled vectors Using Lemma 3, Corollary 1 and observing that all the involved matrices are diagonal, we get u T k v k = 0. Hence vectors u k and v k satisfy the assumptions of Lemma 5.3 in [21]. With U k = diag(u k ) and V k = diag(v k ), we write After multiplying both sides of (20) by (X k Z k ) −1/2 and replacing it in (24) we obtain where the last inequality comes from the fact that ( Lemma 3 again and the definition of µ k to observe that and further rearrange (25): The second norm on the right-hand side of (27) is given by (17). Using the definition of N 2 (θ k ) neighborhood and Lemma 5.4 from [21], for the first norm we get Finally, by substituting (17) and (28) in (27) we obtain the required result.
We are ready to state the main result of this subsection. In Theorem 1 we show that it is possible to choose sufficiently small values for the step-sizesᾱ k andᾱ such that (x k+2 , λ k+2 , z k+2 ) ∈ N 2 (θ k ). Therefore, we ensure that the quasi-Newton step taken after a Newton step remains in the N 2 neighborhood. This implies that all the iterates generated by the algorithm belong to N 2 . The upper bounds for the step-sizes as stated in the theorem are then used to determine the worst-case complexity of Algorithm 1 operating in N 2 . First, recall [21,Theorem 5.6] that it is possible to choose parameters θ k , θ k+1 ∈ (0, 1) and σ k , σ k+1 ∈ (0, 1) so that Theorem 1. Let k + 1 be a quasi-Newton iteration of Algorithm 1. Suppose that Assumptions 1 and 2 hold and that θ k+1 = θ k and σ k+1 = σ k . If the step-sizes in Newton and quasi-Newton iterationsᾱ k andᾱ satisfȳ for σ k ∈ [0, 1) and θ k ∈ (0, 16/25), then Proof. We will first show that for allᾱ k andᾱ satisfying the conditions of the theorem holds. This condition ensures that the Newton step at iteration k + 2 will be successfully performed and the algorithm converges. By inequality (19), condition (31) is satisfied if We will derive bounds to each term on the left-hand side of this inequality in order to find an expression in a form K 1ᾱ 2 − K 2ᾱ , K 1 , K 2 > 0, which will be nonpositive for small values ofᾱ. For the first term, we use the fact that θ k+1 = θ k , σ k+1 = σ k andᾱ k ∈ (0, 1). In addition, we use Lemma 4 to expand µ(ᾱ) and the fact that µ k+1 was calculated in the Newton step k to obtain For the second and third terms, we first apply (29) in Lemma 5 to simplify the bound of γ 1 : Therefore, using (29) again, we derive the following bound to the second term where in the last inequality we used the conditionᾱ k ≤ᾱ ≤ 1 from (30). For the third term in (32), we use the bound obtained in Lemma 7 and analyze each part of it independently. First, sinceᾱ k ≤ᾱ and σ k+1 = σ k , we observe that Using bound (34), assumptionᾱ k ≤ (1 − σ k )/(4σ k ) in (30), and (29) again, we also obtain By combining (36) and (37) in the statement of Lemma 7 and applying (29) once more, we derive a bound to the third term of (32) By (33), (35) and (38), the following bound on expression in (32) is obtained which is negative only ifᾱ ≤ σ k (1 − σ k )/(10(1 − σ k ) + 4), as requested by (30). This bound onᾱ implies a bound onᾱ k , due to conditionᾱ k ≤ᾱ. Therefore, we arrive in the step-size conditions (30) of the theorem. Using (31) we also have that We now show that the new iterate belongs to F . By Assumption 1 and Lemma 1 we know that all iterates remain primal and dual feasible. It remains to show that x k+2 and z k+2 are strictly positive. We follow the same arguments as [17] and adapt them to our case. Suppose by contradiction that x k+2 i ≤ 0 or z k+2 i ≤ 0 hold for some i. By (40), we have that x k+2 i < 0 and z k+2 i < 0 and that implies Since (x k , λ k , z k ) ∈ N 2 (θ k ) andᾱ ≤ 1/4, by (30), we conclude that (1 − θ k )µ k < (9/16)θ k µ k and, therefore, θ k > 16/25, which contradicts the choice of θ k . Hence, (x k+2 , λ k+2 , z k+2 ) belongs to N 2 (θ k ) and the Newton step at iteration k + 2 also falls in the N 2 (θ k ) neighborhood.
For the polynomial convergence of Algorithm 1, we define σ k = σ k+1 = 1 − 0.4/ √ n and θ k = θ k+1 = θ k+2 = 0.4, which satisfy condition (29) [21] and maintain all the previous results. By Lemma 4, µ(ᾱ) ≤ µ k+1 ≤ µ k . So it is enough to look just at the Newton steps, which are easier to analyze. Using (30) it is not hard to see that Therefore, from which the convergence with the worst-case iteration complexity of O(n) can be derived.

The N s neighborhood
Colombo and Gondzio [3] used the symmetric neighborhood N s (γ), defined by for γ ∈ (0, 1), which is related with the N −∞ (γ) neighborhood used in long-step primal-dual interior point algorithms. The idea of the symmetric neighborhood is to add an upper bound on the complementarity pairs, so that their products do not become too large with respect to the average. The authors showed that the worst-case iteration complexity for linear feasible primal-dual interior point methods remains O(n) and the new neighborhood has a better practical interpretation. As HOPDM [7] implements the N s neighborhood and it was used in the numerical experiments in [10] for quasi-Newton IPM, it is natural to ask about the iteration complexity of Algorithm 1 operating in the N s neighborhood. The analysis presented below will follow closely that from Subsection 3.1. We start from an assumption, but it is worth observing that, from [3,21], this assumption holds if the step-size in the Newton direction is sufficiently small:ᾱ k ∈ 0, min 2 3/2 1−γ 1+γ σ k n , 2 3/2 γ 1−γ 1+γ σ k n . Assumption 3. Let γ ∈ (0, 1) and (x k , λ k , z k ) ∈ N s (γ). Let the iterate after a stepᾱ k in Newton direction also satisfy (x k+1 , λ k+1 , z k+1 ) ∈ N s (γ).
Our main goal is to show that the next iterate obtained after a step in the quasi-Newton direction also belongs to N s (γ) if suitable step-sizesᾱ andᾱ k are chosen. To demonstrate this, we will consider lower and upper bounds on the complementarity products in the N s (γ) neighborhood using two possible values of ζ ∈ {γ, 1 γ }. We start the analysis from expanding the complementarity product at the quasi-Newton iteration whereᾱ ∈ (0, 1) is the step-size associated with the quasi-Newton direction. To guarantee that the new iterate belongs to N s (γ), for ζ = γ the expression in (41) should be non-negative and for ζ = 1/γ it should be non-positive, for i = 1, . . . , n. Using (14) and Lemma 4, we rewrite the expression in (41) To deliver the main result of this section we will need a bound for the quasi-Newton term γ 1 defined in Lemma 1 when the algorithm operates in the N s (γ) neighborhood. Lemma 8. Let k + 1 be a quasi-Newton iteration of Algorithm 1 operating in the N s (γ) neighborhood. Suppose that v in Lemma 1 is given by the right-hand side of (4) and Assumption 3 holds. Ifᾱ k ∈ (0, 1] and σ k ∈ [0, 1), then Proof. Using the conditions of the lemma, the definition of µ and the fact that σ k+1 ∈ [0, 1], we have that Sinceᾱ k ∈ (0, 1] and σ k ∈ [0, 1), by defining ρ k = 1 − σ k > 0 we can again use µ k+1 = (1 −ᾱ k (1 − σ k ))µ k and Lemma 2 to ensure that y k > 0. Therefore, by (15), Lemma 2 and the previous result, The next lemma delivers a bound for term ᾱ k ∆x k +ᾱ∆x k+1 i ᾱ k ∆z k +ᾱ∆z k+1 i under Assumption 3. Most of the calculations have already been made in the proof of Lemma 7. Let us mention that Lemmas 7 and 9 can be viewed as the quasi-Newton versions of [21, Lemma 5.10].

Lemma 9. Let k + 1 be a quasi-Newton iteration of Algorithm 1 and suppose that Assumption 3 holds. Then,
Proof. We observe that many arguments used at the beginning of the proof of Lemma 7 remain valid for the N s neighborhood. (Only the bound x k i z k i ≥ (1 − θ k )µ k needs to be replaced with x k i z k i ≥ γµ k .) Then inequalities (25) and (27) are replaced with the following where in the last equality we have used equation (26). We already have the expression for (µ(ᾱ)−µ k )e 2 (see (17)), but we need a bound for the first norm in (43). We observe that Therefore, by (44) and the bound on ∆X k ∆Z k e obtained in [21, Lemma 5.10] (which also holds for N s [3]) since n ≥ 1. Using (17) and (45) in (43), we obtain the desired result.
Using (42) and the bounds obtained so far, we now show that, for sufficiently small step-sizesᾱ k andᾱ, in the Newton and quasi-Newton iterations, respectively, the point (x k+2 , λ k+2 , z k+2 ) also belongs to N s (γ). The upper bounds for the step-sizes delivered by the theorem below are then used to determine the O(n 3 ) iteration worst-case complexity of Algorithm 1 operating in N s (γ).
Proof. By construction we guarantee thatᾱ ≥ 2ᾱ k as needed in (46). We start the proof by setting ξ = γ in (42) and showing that x k+1 +ᾱ∆x k+1 i z k+1 +ᾱ∆z k+1 i − γµ(ᾱ) ≥ 0 for sufficiently small step-sizes. By (42), using Assumption 3, Lemma 9 and [21, Lemma 5.10], we obtain We will rearrange this expression and represent it in a form K 1ᾱ − K 2ᾱ 2 with K 1 , K 2 > 0. Next, we will show that (47) is non-negative for sufficiently small values ofᾱ k andᾱ. To deliver the desired results every term inside the curly brackets in (47) will be bounded.
It remains to show that x k+2 and z k+2 are strictly positive. Similarly to Theorem 1, if, by contradiction, x k+2 i ≤ 0 or z k+2 i ≤ 0 for some i, then we must have x k+2 i < 0 and z k+2 i < 0. On one hand, we already know that x k i z k i ≥ γµ k , by Assumption 3. Hence, by Lemma 9, (51), (52) and similar arguments to those used in this proof 2l .
To deliver the worst-case polynomial complexity, we follow the same arguments as those presented after the proof of Theorem 1 and consider only Newton iterations which are simpler to analyze. By (46), we have that  Because of the close relation of the symmetric neighborhood and the N −∞ (γ) neighborhood, one should expect similar complexity results to that of Theorem 2. However, if the N s (γ) neighborhood is not applied in Lemma 8, then a naive approach is to bound n i=1 x k+1 i z k+1 i 2 by µ 2 k+1 n 2 . Such approach would increase the worst-case iteration complexity to O(n 4 ), as the degree would be increased in equation (51). A different approach to reduce the worst-case polynomial degree when working with the N −∞ (γ) neighborhood is subject to future developments.

Worst-case complexity in the infeasible case
In this section, we analyze worst-case iteration complexity for linear programming problems without assuming feasibility of the starting point. The analysis uses some ideas developed in Section 3 and follows [21,Chapter 6], by adding extra requirements on the computation of step-sizeᾱ k . The proofs are modified to consider the symmetric neighborhood and admit the quasi-Newton steps.
We start by defining r k b and r k c to be the primal and dual infeasibility vectors at point (x k , λ k , z k ) which appear in the right-hand side of equation (4): Then, given an initial guess (x 0 , λ 0 , z 0 ) and parameters β ≥ 1 and γ ∈ (0, 1), the infeasible version of the symmetric neighborhood is defined as follows In order to obtain complexity results,ᾱ k needs to be computed in such a way that the new iterate remains in N s (γ, β) and a sufficient decrease condition for µ k is satisfied. More precisely, given α dec ∈ (0, 1),ᾱ k is the largest value in [0, 1] (or a fixed fraction of it) such that The goal of this section is to show that, assuming that one quasi-Newton step is performed from a point in N s (γ, β) by Algorithm 1, then conditions (56) are satisfied by the new point. When only Newton steps are taken, it is shown in [21, Lemma 6.7] that there is an interval [0,α] such that (56) holds for the N −∞ (γ, β) neighborhood. This is the case if a special starting point is used andα ≥δ/n 2 , whereδ is a constant independent of n. In Lemma 10, we extend those results in order to ensure that the iterates belong to N s (γ, β) when only Newton steps are taken.
First, let us recall some results from [21] that will be used frequently. We define the scalar ν k = which allows us to write vectors r k b and r k c as r k b = ν k r 0 b and r k c = ν k r 0 c , respectively. The parameter σ k used in (8) satisfies 0 < σ min ≤ σ k ≤ σ max < 1 for all k. We also assume that a special initial point is used in Algorithm 1, given by where ξ is such that (x * , z * ) ∞ ≤ ξ, for some primal-dual solution (x * , λ * , z * ). Let k be a Newton iteration, (x k , λ k , z k ) ∈ N s (γ, β) and D k = (X k ) 1/2 (Z k ) −1/2 . Also, let ω = 9β/γ 1/2 be a constant independent of n. When (57) is used as the starting point, the following bounds hold: The proofs for these bounds can be found in Lemmas 6.4 and 6.6 of [21], respectively. We start the analysis from looking at Newton step in N s (γ, β) neighborhood.

Lemma 10.
If (x k , λ k , z k ) ∈ N s (γ, β) and the Newton step is taken at iteration k, then there exists a constantδ independent of n and a valueα ≥δ/n 2 , independent of k, such that for allᾱ k ∈ [0,α] Proof. We only need to address the right inequality of (62), since all the other results were detailed in [21,Lemma 6.7]. Since k is a Newton iteration, by (60), we obtain As (∆x k , ∆λ k , ∆z k ) solves (4) and (x k , λ k , z k ) ∈ N s (γ, β), using a component-wise version of (9) we get Using similar arguments and equation (9) again, we also obtain Using (63), (64) and the fact that n ≥ 1, we obtain The right-hand side of this inequality is non-positive ifᾱ k ≤ σ min (1 − γ) ω 2 n 2 (1 + γ) . By definingδ as the minimum of [21,Lemma 6.7], we obtain the desired result.
Based on the previous paragraphs, we set some assumptions that will be used in the remaining of this section. As usual, at iteration k the Newton step is taken and at iteration k + 1 the quasi-Newton step is made. We assume that (x 0 , λ 0 , z 0 ) satisfies (57), both (x k , λ k , z k ) andᾱ k satisfy (56), and (x k+1 , λ k+1 , z k+1 ) ∈ N s (γ, β). In our analysis, we will use extensively various well known results for Newton steps, such as properties (58)-(62). We observe that, while in the standard Newton approach one has to compute bounds for ∆x k T ∆z k , in the quasi-Newton approach (see for example (13)), it is necessary to additionally get bounds for γ 1 and (ᾱ k ∆x k +ᾱ∆x k+1 ) T (ᾱ k ∆z k +ᾱ∆z k+1 ). In the next lemma, we give a bound for γ 1 without assuming feasibility of iterates.
Lemma 11. Let k + 1 be a quasi-Newton iteration of Algorithm 1. Suppose that v in Lemma 1 is given by the right-hand side of (4) at iteration k + 1 andᾱ k ∈ (0, 1]. Then, there exists a constant C 3 ≥ 1, independent of n and k such that Proof. Using scalar ν k+1 defined in the beginning of this section, vector v at iteration k + 1 is given by We observe that, since e T (µ k+1 e − X k+1 Z k+1 e) = 0, by (56) we obtain Taking the 2-norm in (66) and using (59) and (67), we get By defining ρ k = α dec > 0, we can see that condition (56) results in the sufficient decrease condition of Lemma 2. Sinceᾱ k ∈ (0, 1], this ensures that y k > 0. Using Lemma 2, (56) and (68) in equation (15), we conclude that . In the last inequality we used (56) to yield µ k+1 /µ k ≤ 1.
Our goal now is to bound the term (ᾱ k ∆x k +ᾱ∆x k+1 ) T (ᾱ k ∆z k +ᾱ∆z k+1 ) in (13). We follow the same approach as [21] but compute, instead, bounds of (D k ) −1 (ᾱ k ∆x k +ᾱ∆x k+1 ) and D k (ᾱ k ∆z k +ᾱ∆z k+1 ), where D k = (X k ) 1/2 (Z k ) −1/2 . It is important to note that matrix D k is related to the Newton iteration, hence we can use Lemma 1 and the enjoyable properties of the true Jacobian.
First, we show what happens when matrix A multiplies the combined directionᾱ k ∆x k +ᾱ∆x k+1 : where x * is the primal solution of (1) used to define constant ξ in (57). By definingν k = −(ᾱ kᾱ −ᾱ k −ᾱ)ν k = (1 − (1 −ᾱ k )(1 −ᾱ))ν k , we can see that the pointx =ᾱ k ∆x k +ᾱ∆x k+1 +ν k (x 0 − x * ) is such that Ax = 0. Similar arguments can be used to show thatẑ =ᾱ k ∆z k +ᾱ∆z k+1 +ν k (z 0 − z * ) andλ =ᾱ k ∆λ k +ᾱ∆λ k+1 +ν k (λ 0 − λ * ) satisfy A Tλ +ẑ = 0. Therefore, it is not hard to see that We now multiply the third row of the coefficient matrix in (4) by xλẑ to obtain Multiplying this equation on both sides by (X k Z k ) −1/2 , we get Using (69) and (70), we conclude that Using this relation we get and observe that the same bound holds for D kẑ . In the next two lemmas, we compute the bounds for the terms in the right-hand side of (71).
Lemma 12. Let k + 1 be a quasi-Newton iteration of Algorithm 1. Let ω and C 3 be the constants (independent of k and n) defined in equation (60) and Lemma 11,respectively. Then First, we observe that the diagonal matrices ∆X k and ∆Z k were computed in the Newton iteration and use (60) to obtain k . Since both (D k ) −1 and ∆X k are diagonal matrices, using the property of the induced 2-norm of matrices, we get Hence, ∆X k ∆Z k e ≤ ω 2 n 2 µ k .
Now, we use equation (21) (which does not depend on the feasibility of the iterate or the type of the neighborhood) to expand the desired expression in the statement of this Lemma. Additionally, we use equation (56) to bound µ k+1 and X k Z k e as well as Lemma 11 and equation (72) to derive: Lemma 13. Let k + 1 be a quasi-Newton iteration of Algorithm 1. Let ω and C 3 be the constants (independent of k and n) defined before. Then, We will only consider (D k ) −1 (ᾱ k ∆x k +ᾱ∆x k+1 ) , since getting a bound for D k (ᾱ k ∆z k +ᾱ∆z k+1 ) follows the same arguments. We use the definition ofx (see (69)), add and subtractν k (D k ) −1 (x 0 − x * ) inside the norm, and then use the triangle inequality and (71) to obtain where another termν k D k (z 0 − z * ) was added in the last inequality. We already have bounds for the first term in the right-hand side of (73), by Lemma 12. From [21, Lemma 6.6] we know that and, similarly, D k (z 0 − z * ) ≤ ξ (X k Z k ) −1/2 x k 1 . We recall thatν k = (1 − (1 −ᾱ k )(1 −ᾱ))ν k and apply all these inequalities together with (58) to obtain which provides bounds for the last two terms in (73). Therefore, using Lemma 12 and the above inequality we write kᾱ γ −1/2 n 5/2 µ 1/2 k and the lemma is proved.
If we restrict the choices ofᾱ k andᾱ, then the bounds obtained in Lemma 13 can be significantly simplified as shown in the corollary below.
3ᾱ k , then there exists a constant C 6 , independent of n, such that Proof. Using the bounds onᾱ andᾱ k assumed in the lemma we obtain ( and the conclusion follows by defining C 6 = (σ max + γ −1 + 8β)(1 + 2C −1 We are now ready to prove the polynomial worst-case iteration complexity of Algorithm 1 in the infeasible case. Lemma 10 dealt with the Newton step of the method. Theorems 3 and 4 provide the results for the quasi-Newton step. Theorem 3. Suppose that k + 1 is a quasi-Newton step of Algorithm 1 and all the hypotheses of this section hold. If the following condition holds α dec + σ max ≤ 1 − σ min (74) then, there exists a constant C 5 , independent of n, such that, if α k ∈ 0, min 1, (1 + C 3 )ω 2 −1 , C 5 n 5 andᾱ ∈ (C 3 C 5 ) −1 n 5ᾱ2 then the iterate after the quasi-Newton step satisfies [x k +ᾱ k ∆x k ] i [z k +ᾱ k ∆z k ] i ≥ γµ(ᾱ) [x k +ᾱ k ∆x k ] i [z k +ᾱ k ∆z k ] i ≤ 1 γ µ(ᾱ).
For sufficiently large n, we clearly have min 1, (1 + C 3 )ω 2 −1 , C 5 n 5 = C 5 n 5 <δ n 2 , whereδ comes from Lemma 10 and this guarantees that the hypotheses of this section hold. By (75) we have that bothᾱ k andᾱ are greater than or equal to for all k = 0, 1, . . . , which yields O(n 5 ) worst-case iteration complexity for the infeasible case of Algorithm 1. Unlike in the feasible case, because of the close relation of the symmetric neighborhood and the N −∞ (γ, β) neighborhood, we expect that a similar complexity result to that of Theorem 3 should hold for the algorithm which operates in N −∞ (γ, β). However, to save space we do not derive it here.

Conclusions and final observations
This work provided theoretical tools to analyze the worst-case iteration complexity of quasi-Newton primal-dual interior point algorithms. A simplified algorithm was considered, which consisted of alternating Newton and quasi-Newton steps. The quasi-Newton approach was based on the Broyden "bad" low-rank update of the inverse of the unreduced system. Feasible and infeasible algorithms and well established neighborhoods of the central path have been considered.
The results showed that in all cases, the degree of the polynomial in the worst-case result has increased. This behavior has already been observed by [10], where the number of overall iterations increased, but the number of factorizations of the true unreduced system actually decreased.
Complexity results might also be obtained in the feasible case for convex quadratic programming problems and the N 2 (θ) neighborhood, following [8], for example. However, since Lemma 3 does not hold anymore and thus γ 1 is not eliminated from µ(ᾱ), a further study needs to be carried out.
Another interesting and complicated question is the case ℓ > 1, where we allow more than one quasi-Newton step after a Newton step. The authors in [10] observed a trade-off between the number of consecutive quasi-Newton iterations and the overall speed of convergence of the algorithm. Numerical results showed that ℓ = 5 is a reasonable choice. We expect that similar worst-case complexity bounds should be obtained for the case ℓ > 1, with possibly higher degrees of polynomials, but this still remains an open question.
The worst-case complexity results obtained when the N s neighborhood was considered in both feasible and infeasible cases seem rather pessimistic. The high polynomial degrees originate from Lemma 9 for the feasible and from Lemma 13 for the infeasible cases. Finding new ways of reducing the degrees of such expressions would definitely result in better complexity results, since all the other terms are at most of order n.