Convergence of a Class of Inexact Interior-Point Algorithms for Linear Programs.
We present a convergence analysis for a class of inexact infeasible-interior-point methods for solving linear programs. The main feature of inexact methods is that the linear systems defining the search direction at each interior-point iteration need not be solved to high accuracy. More precisely, we allow that these linear systems are only solved to a moderate accuracy in the residual, but no assumptions are made on the accuracy of the search direction in the search space. In particular, our analysis does not require that feasibility is maintained even if the initial iterate happened to be a feasible solution of the linear program. We also present a few numerical examples to illustrate the effect of using inexact search directions on the number of interior-point iterations.
Introduction
Since the publication [6] of Karmarkar's original interior-point algorithm for linear programs, numerous
variants of the method have been developed; see, e.g., Kojima, Mizuno, and Yoshise [8],
Megiddo [13], Monteiro and Adler [17], and Tanabe [21]. Especially interesting is the algorithm by
Kojima, Megiddo, and Mizuno [7], since it is also practically efficient; see, e.g., [9].
Numerical Analysis Manuscript No. 96-16, Bell Laboratories, Murray Hill, New Jersey, September 1996. Available
on WWW at http://cm.bell-labs.com/cs/doc/96.
1.1 Inexact Interior-Point Methods
The motivation for the study in this paper is to provide a theoretical justification for recent modifications
of the algorithm [7] in which the search directions are computed only inexactly at each
interior-point iteration. One feature of such inexact interior-point methods is that they allow the
use of iterative techniques to approximately solve the linear systems defining the search directions.
For example, in [3], we employ a preconditioned QMR algorithm for the iterative solution of these
linear systems.
First, we describe a modification of the algorithm in [7] that is adapted to inexact computations
of the search directions, and then we show its global convergence. Our goal is to design an algorithm
such that the number of iterations to solve a given linear program under "reasonably" inexact
computations is not much higher than for the algorithm in [9]. We are aware that the algorithm
in [9] has been improved in the meantime; see, e.g., [10] and many other more recent publications.
However, we are not aware of any more efficient algorithm that is readily extended to inexact
computations.
In some applications, an initial strictly dual-feasible solution of the linear program is known. For
such problems, in [18], an implementation of a primal-infeasible dual-feasible interior-point method
is presented and results of numerical experiments are reported. Our approach originates from [3],
and it is distinctly different from [18]. In particular, the method of our paper also converges even
if there is no strictly dual-feasible solution.
1.2 Outline
The outline of this article is as follows. In Section 2, we formulate the linear program in standard
form and state our assumptions. In Section 3, we describe the class of inexact interior-point
methods studied in this paper. In Section 4, we discuss issues in defining appropriate inexact
search directions. In particular, we explain why we assume a certain accuracy requirement, and
how we can achieve it. In Section 5, we establish global convergence. In Section 6, we make some
concluding remarks. Finally, in an Appendix, we define a certain condition number of a linear
program, and explain how it relates to our analysis in Section 5.
1.3 Notation
Throughout the paper, we use the following notation. Lower-case letters are used for vectors, and
upper-case letters are used for matrices. For a vector denoted by a lower-case letter, the upper-case
version means the diagonal matrix whose diagonal elements are the components of the vector. For
example,
is the diagonal matrix associated with the vector . By e we denote
the vector . The vector norm
x T x is always the Euclidean norm,
and denotes the corresponding Euclidean matrix norm. The notation
means that all components of the vector x are positive (respectively,
nonnegative). Instead of x - 0, we also write x 2 IR n
. The same symbol, 0, will be used to denote
the number zero, the zero vector, and the zero matrix. When 0 denotes the zero vector or the
zero matrix, then its length or dimension, respectively, will always be clear from the context. The
capital letter I is used to denote the (square) identity matrix; its size will always be apparent from
the context as well.
2. The Problem
2.1 Formulation and Optimality Conditions
We consider linear programs in the standard form
minimize x2IR n
c T x subject to
(1)
We assume that A is a given real m \Theta n matrix, and b and c are given real vectors of length m and
n, respectively. The dual problem to (1) can be written as
maximize
subject to A T y
(2)
The Karush-Kuhn-Tucker optimality conditions for the primal-dual problems (1) and (2) can be
stated as follows:
Equations (3) and (4) are referred to as primal and dual feasibility, respectively, and the nonlinear
equation (5) is called the complementarity condition.
2.2 Assumptions
In the sequel, we always assume that A has full row rank. Theoretically, this assumption could be
enforced by computing a singular value decomposition of A, but given a sparse constraint matrix
A, this is not done in practice. Instead, a given linear program is typically preprocessed before
an interior-point method for its solution is applied; see, e.g., [1, 5]. For example, simple heuristics
are used to detect fixed variables, empty rows and columns, redundancies, and certain implicit
dependencies. These heuristics are often sufficient to reduce the linear program not only in size
but also to a form that the matrix A has linearly independent rows, even if the constraint matrix
of (1) in its initial form has not.
The convergence analysis of our inexact interior-point method depends on the smallest singular
value, oe min (A), of A. The smaller oe min (A), the higher the accuracy that is required at each step
of the method. We now briefly discuss the size of oe min (A), assuming that all rows of A are scaled
to Euclidean norm 1. We then have 1 - oe
is the largest singular
value of A. When oe min (A) is very small, the linear program (1) is not well-conditioned. However,
assuming that oe min (A) is large (say greater than 0:1) does not imply that the problem (1) is well-
conditioned. In fact, also for 0 - oe min (A) - 1, the linear program may be degenerate, and then the
pair of primal-dual problems has infinite condition see the Appendix. In particular, prior
to solving a given linear program (1), the problem could always be transformed into an equivalent
linear program with a constraint matrix A 0 such that oe min since A is assumed to
have full rank, the matrix A T admits a QR factorization of the form
where Q is an orthogonal n \Theta n matrix with Q 1 denoting its first m columns, and R 1 is a nonsingular
upper triangular matrix. By means of (7), the equality constraints in (1) can be rewritten
as follows:
The reformulation Q T
1 b of the equality constraints is obviously well-conditioned with
oe min (Q T
theoretically, one might first compute the QR factorization (7) as well as
R \GammaT
(possibly with higher accuracy than the standard machine precision), and then solve the
transformed linear program. Unfortunately, the reformulation (8) typically destroys the sparsity of
the constraint matrix, as Q T
may have considerably more nonzero entries than A, and therefore,
this reformulation is not used in practice.
We do, however, assume in the sequel that a positive lower bound,
for the smallest singular value of the original constraint matrix A is available.
3. A Class of Inexact Interior-Point Algorithms
The class of algorithms considered in this paper generates a sequence of iterates that simultaneously
converge to feasibility and optimality. To guarantee convergence, we require that all iterates lie in
a certain set N .
3.1 An Infinity-Norm-Neighborhood of the Perturbed Path
The set N depends on some initial point
and a real number - ? 0 be given such that x
For a real number fi ? 1, a neighborhood N defined by the set of equalities
(11) and the two sets of inequalities (12) and (13) below. A point
belongs to N (fi) if there exist
vectors
such that
and
e
A closely-related neighborhood was proposed in [7]. The new element in the above definition is
the perturbation vectors ~ b; ~ c in (11) and (12); in [7], the vectors ~ b and ~
c are zero.
Unless some good initial estimate the solution of problems (1) and (2) is available,
we will think of - as a large number, - AE 0, and of x large vectors. The set N is a
"wide" infinity-norm-neighborhood of the "perturbed path"
It is well known that there exists a unique solution (x; of (14)-(16) for each t 2 (0; 1], provided
(3), (4), and (6) have a solution and A has full rank; see, e.g., [8, 13]. Above, "infinity-norm" refers
to the neighborhood with respect to the nonlinear equation
-e, while the deviation from
the linear equalities is measured in the Euclidean norm. The requirement on the size of ~ b in (12)
could be weakened slightly by requiring kR \GammaT~ bk -=2 where R 1 is as in (7) and (8). However,
our goal is to formulate a method that does not require the QR factorization (7). To emphasize
this principle, we formulate the neighborhood in terms of the lower bound (9), oe, on the smallest
singular value of A.
3.2 The Newton Step
The class of interior-point algorithms considered in this paper can be viewed as variants of damped
Newton's method applied to the system of equations
c;
-e:
are nonnegative numbers. They are adjusted at each iteration of the interior-point
method such that they converge to zero in the course of the algorithm. The step lengths in
the damped Newton method are suitably chosen to guarantee that (x;
Applying Newton's method to (17)-(19) yields the system of linear equations
A T \Deltay
for the Newton correction (\Deltax; \Deltay; \Deltas). We now denote by
the right-hand sides in (20)-(22). Then, after a reordering, the linear system (20)-(22) can be
written in the form 2X S 0
I 0 A T
\Deltas
\Deltax
\Deltay5 =4
with a 3 \Theta 3-block matrix. In our algorithm, the systems (24) are not solved exactly but only to a
given accuracy in the primal and dual residual,
As we will explain in Section 4, it is meaningful and easily possible to solve the last equation in
exactly (up to machine precision), even when allowing possibly larger errors " x and " s for (20)
and (21). In Section 4, we will also define appropriate values for " x and " s .
3.3 An Infeasible-Interior-Point Algorithm
Let be the iterate at step k of the interior-point method, and
let (\Deltax; \Deltay; \Deltas) be the solution of the linear system (24). Then, an improved iterate
at the next step, k + 1, is obtained via
where the step lengths ff x and ff s are chosen such that
is a parameter that can be chosen very close to 1. The sole purpose
of this parameter is to ensure that the possibility of ff
The class of algorithms considered in this paper is designed to support the iterative solution of
the linear systems (24) (as in [3]). The concept of an algorithm using an iterative solver is based
on the following principles.
predictor-corrector strategy is applied since solving for the corrector step is approximately
as expensive as solving for the predictor step. (This is in contrast to the standard methods
based on direct solvers that use the same sparse factorization of some reduced version of the
coefficient matrix of (24) for the predictor and the corrector step.)
ffl Relatively short steps (i.e., not the usual 99.995% of the way to the boundary) are taken. This
is done for two reasons. First, if the search direction is computed to low relative accuracy, a
long step is not meaningful, and secondly, very long steps might worsen the condition number
of the linear systems.
The stopping criterion in the interior-point algorithm is based on the error measure
ck
Algorithm 1 (Inexact infeasible-interior-point algorithm)
INPUT: Parameters and an initial point
For
If
by solving the linear system (24) (with right-hand side (23) and
to an error satisfying (25) with "
Choose ff k
x \Deltax;
s \Deltas:
3 to
-n
Next, we specify the choice of the parameters - k
s and the step lengths ff k
s in
the above algorithm. This choice guarantees that all iterates remain in the neighborhood N (fi);
see Theorem 1.
3.4 Specifications of the Algorithm
The quantities ff k
s , and - k
3 depend on three additional input parameters - 1 , - 2 , and - 3 of the
Step
We set
and - k
For later use, we note that
The number - k
3 depends on an additional quantity fl k . First, the number fl k is chosen such that
min
and then we set - k
. The possibility of a small choice of fl k for small ffi(x k ; y is to allow
for superlinear convergence.
Step
The accuracy in (24) for solving A(x is chosen as
x := -oen
for the primal residual, and as
s := -n
for the dual residual. Here,
x
s
where ff ant
x and ff ant
s are anticipated values for the actual step lengths ff k
x and ff k
s , respectively.
These values need to be chosen such that
In general, the values of ff ant
x and ff ant
s will depend on k. Given any approximation or guess
(\Deltax; \Deltay; \Deltas) for the solution of (24), values for ff ant
x and ff ant
s can be computed easily, for example,
by using the step-length strategy (35) and (36) below. With these values of ff ant
x and ff ant
s , one then
s in (29) and (30), and subsequently, one checks if the current approximation
(\Deltax; \Deltay; \Deltas) satisfies the accuracy requirements (29) and (30). If yes, (\Deltax; \Deltay; \Deltas) is accepted
as a sufficiently accurate solution of (24). If not, one continues the iterative method used for the
inexact solution of (24), until the stopping criteria (29) and (30) (with appropriately updated values
of ff ant
x and ff ant
s are met.
Step
To determine the step lengths ff k
s at iteration step k, we first compute -
ff k as the maximum
such that the point
for all 0 - ff -
ff. As we will see below, the definition of " k
s in (29) and (30) guarantees
that (12) holds as well, and thus,
Furthermore, we
will show that -
ff k is always positive and bounded away from zero. Then ff k
are chosen
as to satisfy i
x \Deltax
s \Deltas
as well as
(fi). The rules (33) and (34) were proposed in
[7], and there, a cheap way of computing -
ff k is also discussed.
ffl In [3], an adaptive choice of the step lengths ff k
s was proposed using the following strategy
that is based on two further parameters - 4 - 1 and - 5 2 (0; 1) with suggested default values
First, one determines
and chooses ae k such that
Then one sets
x and ff k
The motivation for this strategy is that fl k being small implies a large reduction in comple-
mentarity, and -
x and -
s being close to 1 implies a large reduction of infeasibility by the
full Newton step. If both is the case, a large ae k seems profitable even though it may bring
the iterate close to the boundary of the feasible set. Therefore, ff k
s may be a practically
efficient choice. When - 3 ! 1 is close to 1, the constraint (34) on the step lengths is fairly
weak, and the above choice of ff k x ; ff k s may satisfy (34) for most iterations k. If (34) is not
satisfied or if the result does not lie in N (fi), then step lengths that do satisfy these criteria
may be computed as in [7].
The main computational effort in Algorithm 1 lies in the solution of the systems (24). Before
solving these systems, we first employ the stable reduction described in Section 4.1 below to reduce
the 3 \Theta 3-block system to a smaller 2 \Theta 2-block system. The smaller systems can be solved by either
a direct method (see, e.g., [2, 22]), or by an iterative method, such as the one proposed in [3].
3.5 The Main Result
The main result, which is stated as Theorem 1 in Section 5 below, is that all iterates generated
by Algorithm 1 remain in the neighborhood N (fi). Furthermore, the iterates either converge to
the set of optimal solutions if this set is nonempty, or it is detected that no optimal solution exists
whose norm is less than a given bound.
The result that the iterates remain in N (fi) is a consequence of the following two features of
Algorithm 1. (i) The first two equality constraints in (11) are trivially maintained by a suitable
update of ~ b and ~ c, while the third equality constraint in (11) is trivially satisfied due the choice of
3 in Algorithm 1. (ii) The norm constraints (12) on ~ b and ~ c are guaranteed by a proper choice of
" x and " s , and the inequality constraints (13) on Xs and on the t i 's are maintained by selecting
suitable parameters
4. Inexact Search Directions
In this section, we discuss how the accuracy requirement in Step 3) of Algorithm 1 can be satisfied,
and why we have chosen this type of accuracy requirement.
4.1 A Stable Reduction of the 3 \Theta 3-Block System
To show how a search direction satisfying an accuracy requirement of the form (25) can be obtained
in practice, we first review the stable reduction (introduced in [3]) of the 3 \Theta 3-block system (24)
to a linear system with the 2 \Theta 2-block matrix
The basis of the stable reduction is a partition 1 of the vectors x and s into two parts x 1 , x 2 and
respectively, such that x 1 Moreover, for the ease of notation, we assume
that x 1 and s 1 are the leading entries of x and s, i.e.,
and
In this subsection, we do not refer to specific components of the vectors x or s; instead x 1 and s 2
always denote vectors containing all "large" components of x and s, respectively. We partition the
constraint matrix A conforming with (37), and write
Similarly, (37) induces the partition
and
of the residual vectors p and q defined in (23). As usual, we denote by X 1 , X 2 , S 1 , and S 2 the
diagonal matrices associated with the vectors x 1 , x 2 , s 1 , and s 2 , respectively.
Instead of the original system (24), we propose to solve the equivalent reduced 2 \Theta 2-block
system -
\Deltay
~ r
~
where
~ r :=
and ~
Furthermore, the new variable u is connected with x 1 and s 2 via the relation
\GammaX
Using (37)-(39), (41), and (42), together with
one readily verifies that the reduced 2 \Theta 2-block system (40) is indeed equivalent to the original
system (24).
We consider interior-point methods that solve the reduced system (40) (iteratively) to a low
relative accuracy, and thus obtain only an approximate solution. Given such an approximate
solution u, \Deltay of (40), the quantity \Deltax 1 is computed by using the relation (42). Note that
We stress that the choice of this partition may be changed, for example, to further reduce the primal residual
at the cost of the dual residual. In preliminary numerical experiments, it was our observation that a partitioning
according to x1 - s1=10 and x2 ! s2=10 might be slightly more efficient.
computing \Deltas 2 from (42) would involve the inverse of the "small" matrix X 2 . This can be avoided
by computing \Deltas 2 via
and then obtaining \Deltax 2 and \Deltas 1 from (43). In the above computations, only the inverses of the
"large" matrices X 1 and S 2 are used, but not X \Gamma1
2 or S \Gamma1
1 . Furthermore, we remark that the
partition of x into x 1 and x 2 has nothing to do with guessing active indices. The partition is
merely done to improve the stability for solving the current linear system, it does not involve any
additional computational cost, and there is no loss if the partition changes at each iteration.
4.2 Motivation of Inexactness Requirements
Next, we further motivate our requirement in Step 3) of Algorithm 1 on the accuracy of the
computation of the search directions in (25). We implicitly assume that we could compute the search
direction to arbitrarily high precision, but that we are not willing to invest more computational
effort for more accuracy than necessary.
To guarantee convergence of the interior-point method even for inexact search directions, we
need a careful analysis of the errors introduced by the inexactness of the search direction.
Assume for the moment that we set - in (24). If (\Deltax; \Deltas; \Deltay) solve the linear
systems (24) exactly, then, by the linearity of r, the primal residual r + at
Thus, the norm of the primal residual r is reduced by the factor 1 \Gamma ff x in the case of exact solution
of (24). If we allow for some tolerance " in the approximate solution of (24) by requiring that \Deltax
is merely computed to a certain relative accuracy in the residual,
instead of straightforward computations, we have
inexact
The result r
inexact obtained from inexact computations is almost as good as the result obtained
from exact computations.
Above inequality also shows that the size of ff x is important for the reduction of r. It is obvious
that the maximum feasible step lengths -
ff s are very sensitive functions of the accuracy of the
third residual q in (19), especially if the current iterate is close to the boundary of the positive
orthant. We therefore assume, that the stable reduction (43) above is used so that q
holds true.
Of course, analogous considerations as above for r and r + also hold for p and p + . Finally, we
remark that criterion (44) is not satisfactory in this form for the following reason. An iterative
approach for solving the linear systems requires that we allow for a low relative accuracy even if the
initial residual happened to be zero, i.e., if Ax criterion (44) does not allow any increase
in the residual. We therefore do not use criterion (44) but a more appropriate modification of (44),
which is discussed next.
4.3 Analysis of the Perturbed Feasible Sets
We consider the perturbed constraints,
that are part of the relations defining N (fi). In the following proposition, we give a condition on
the perturbation vectors ~ b and ~ c that still guarantees the feasibility of (46)-(48). This result is
crucial for the convergence analysis in the next section. Furthermore, we remark that the lower
bound (9) on the smallest singular value of A is used only in the proof of this result.
Proposition 1 Assume that the problems (1) and (2) have an optimal solution. If x
then the constraints (46)-(48) have a feasible solution for all perturbation vectors ~ b; ~
c with
and for all In particular, there exists a feasible solution (~x; ~
s) of (46)-(48) with
~
x -2 -t 1 e and ~
e:
Moreover, if t
can be chosen to satisfy
Proof: Let the point (-x; -
s) be a feasible solution of (1) and (2), i.e., it satisfies the constraints
(3), (4), and (6). We need to construct a feasible solution (~x; ~
s) of (46)-(48). To this
end, we first set h := A +~ b, where A + denotes the pseudoinverse of A (see, e.g., [20]). The vector h
clearly satisfies
Moreover, since and by the first inequality in (49), we have
We now define (~x; ~
~
Using (3), (4), (52), and the definition of - b; - c in (10), one readily verifies that the point (~x; ~
satisfies (46) and (47).
It remains to show (50) and (51). Note that (53) implies h - \Gamma- e=2. Together with - x - 0,
it follows from (54) that
which is just the first inequality in (50). Note that, by the second inequality in (49), we have
and thus ~
Together with (55), this implies the second inequality in (50). Finally, (51)
follows from (54) (with t
Here, we would like to stress our observation of Section 2, that theoretically, one can transform
problem (1) to an equivalent problem such that prior to applying the interior-point method,
and practically, when this transformation is not appropriate, we may still anticipate that 1=oe
is moderate even if (1) happened to be poorly conditioned. In particular, degeneracy or near-
degeneracy do not depend on oe.
4.4 Required Accuracy
In the course of Algorithm 1, the vectors - are reduced to t k- b
and t k
respectively. Furthermore, by the inexactness of the search direction, "new" perturbation
vectors ~ b k and ~ c k are introduced. More precisely, we have
where ~ b k and ~ c k are defined as follows:
We now show that the vectors (56) remain bounded, provided that in Step 3) of Algorithm 1 the
search directions are computed to the accuracies "
s defined in (29) and (30),
respectively.
Proposition 2 For all iterates computed by Algorithm 1, the associated perturbations vectors ~ b k
and ~ c k defined in (56) satisfy
Proof: The proof is by induction on k. For
A T y thus (57) is trivially satisfied for
Now assume that (57) holds for some k - 0. We show that then (57) is also satisfied for k + 1.
We will only verify the first relation in (57); the second relation is treated analogously. Recall that
where we omitted indices in
1 . Using (56) (with k replaced by
(58), we obtain
By the first relations in (56) and (57), respectively, we have
Taking norms in (59) and using (29) and (60), it follows that
Next, recall from (27) and (32) that - 1 - t 1 and ff - ff ant
x . Thus, together with the definition of
1 in (31), it follows that
Using the definition of " x in (29), we further obtain
Inserting (62) into (61) and using t k+1
This is just the desired first relation in (57) (with k replaced by k 1).
Remark: When - k
the upper bound (62) for " k
x reduces to
x
To relate this to some statement about relative accuracy, we point out that the norm of the right-hand
side r in (23) is bounded by 2t k
1). Note that - and k - bk are known and k - bk=-
can easily be bounded for a suitable choice of x 0 . If k - bk is zero or very small compared to -, the
concept (44) of solving for a certain relative accuracy is more restrictive than (29), and if k - bk is
large, say k
roughly means that
x
This is in the spirit of (44), except for the additional factor t k
1. (Likewise for p.) Conditions
implying that for very small t k
2 a high relative accuracy is required, coincides with
the desire of a higher relative accuracy in the final iterations of the interior-point method to obtain
superlinear convergence.
5. Convergence Analysis
In this section, we present a proof of convergence of the inexact interior-point Algorithm 1.
Proposition 3 If the primal and dual linear programs (1) and (2) have an optimal solution, then
the x- and s-components of the neighborhood N (fi) are bounded, i.e., there exists a constant M 1
such that
Furthermore, if (1) and (2) have a solution (-x; -
s) with
ks
for all (x;
Proof: This proposition was used in slightly different form in [14, 15]. Since the form stated here
is more general than in [14, 15], we give a complete proof.
Let be an
optimal solution of (1) and (2). By Proposition 1, there exists a point (~x; ~
s) such that
c; ~
We now set
and hence, using (65),
Straightforward algebra leads to
x
-:
This concludes the proof of (63).
Now assume that (1) and (2) have a solution (-x; -
!. We can then
choose x opt and s opt as minimum-norm solutions in the inequality (66). Also, note that, by (51),
~
. The proof of (64) then follows from the inequality (66).
Proposition 4 For all fi ? 1, there is a - ? 0 such that
e for all (x;
Proof: The proof is straightforward by using the definition of N (fi).
Proposition there is an M 3 ! 1 such
that for any point generated by Algorithm 1 with
and ks , the bound k\Deltaxk; k\Deltask - M 3 holds. Here, \Deltax and \Deltas are determined
from (24) and (23) for a given accuracy specified by (25).
Proof: Note that kA T y ck is bounded for all iterates. From boundedness of ks k k thus
follows boundedness of k. Note that " x and " s are bounded as well, and \Deltax, \Deltas are the exact
solution of a system (24) with r and q perturbed by a (bounded) quantity of size " x and " s . Thus,
we may assume that the right-hand sides (23) are bounded, and from Proposition 4 and our full-rank
assumption on A follows that the inverses of the matrices in (24) are uniformly bounded. (The
inverse exists for any x ? 0 and s ? 0, and is a continuous function of x and s!) This concludes
the proof.
Next, we state and prove our main result.
Theorem 1 The iterates generated by Algorithm 1 are contained in N (fi). If the iterates generated
by Algorithm 1 are unbounded, then the linear program (1) has no solution. If all iterates generated
by Algorithm 1 are bounded, then the stopping criterion in Step 1) of Algorithm 1 will be satisfied
after a finite number of iterations.
Proof: This proof is based on ideas in Kojima et. al. [7].
By Proposition 2, the vectors ~ b and ~ c associated to x k and s k satisfy (57), and hence (12).
Relation (13) is guaranteed by the update in Step 4) of Algorithm 1, and hence, the iterates remain
in N (fi).
If the iterates generated by Algorithm 1 are unbounded, then, in view of Proposition 3, the
linear program (1) has no solution.
We now assume that the sequence of iterates is bounded, say kx ks
We show that then Algorithm 1 finds an approximate optimal solution after a finite
number of iterations.
Suppose that t k
by (13), and by the
bounds (12), also ck ! 0. By (11), we also have
Therefore, in view of (26), the stopping criterion of Algorithm 1 is satisfied after a finite
number of steps.
Conversely, as long as Algorithm 1 is running, we have Hence there is some
(for which Algorithm 1 does not halt). By
Proposition 4, there is some - 3 ? 0 such that x k - 3 e and s k - 3 e.
Also, by (28), fl k is bounded away from zero, without loss of generality,
Further, Proposition 5 is applicable, and there is some M 3 ! 1 such that k\Deltax k k; k\Deltas k k - M 3 for
all k. (Without loss of generality, we assume M 3 - 1.) Finally, in view of (28), we may assume
that independent of k.
To conclude the proof, we show that -
ff k is bounded away from zero, i.e., that there is some
Note that both, ff k
x and ff k
s can be chosen equal to -
ff k . Thus (67) guarantees that a step length as
required by (33) can indeed be found, and (34) ensures that t k
converges to zero.
Hence, showing (67) indeed completes the proof of convergence.
The computation of -
ff k depends on (13) and (33). By the definition of " x and " s , all the
equalities defining N (fi), as well as the bounds on ~ b and ~ c are maintained for all ff 2 [0; -
consider some iteration k, and verify that all bounds are maintained.
Bound on Xs=(t 3 -)
For completeness, we write the proof from [7] in our notation.
ff\Deltas. Clearly, since x k
(independently of k), there is some ff 1 ? 0 such that x(ff); s(ff) ? 0 for 0 - ff - ff 1 . We consider
This number is positive if, and only if, the first relation of (13) is satisfied. From (45) it follows
that
fin
with
We also have
Thus we can bound (68) from below by
\Deltax T \Deltas
Note that j\Deltax i \Deltas
3 and j\Deltax T \Deltas=nj - M 2
3 . It follows that for all
the relation (68) holds provided that it holds for
Bounds linking t 1 and t 3
First we show that there is some ff 3 ? 0 (independent of k) such that t 1
assuming that t 1 are
defined in the obvious way,
n-
n-
Recall that, in Algorithm 1, either
g. We refer to the first instance as "Case a)" and to the latter as "Case b)". In
Case a), it follows that
calculations yield that t 1
(j
\Deltax T \Deltas
n-
In Case b), it follows from fit 1 -
for some - 6 ? 0 independent of k. Since fit further conclude that t 1
3 . The minimum of these bounds yields ff 3 .
Next we show that t 3 This expression is
written as
\Deltax T \Deltas
n-
In Case a), the quantity - 3 bounded away from zero, and in Case b), fit
away from zero, so that the desired result follows again in a straightforward manner.
Finally, the bounds relating t 2 and t 3 are analogous to the ones between t 1 and t 3 .
We close by noting that the rules (33), (34) are identical to [7], so that the analysis there fully
applies to our algorithm, and guarantees, that the step lengths are bounded away from zero as long
as the algorithm does not stop.
6. Concluding Remarks
In this paper, we proved global convergence of the infeasible-interior-point method stated as Algorithm
1. This method allows that the linear systems at each iteration are solved only to moderate
relative accuracy in the residual. Here, the accuracy requirement depends on an estimate oe for the
smallest singular value oe min (A) of the constraint matrix A. The numbers oe or oe min (A) are not used
elsewhere in the algorithm.
Algorithm 1 is very similar to the algorithm in the paper [7], which, however does not provide
an analysis for the case where the search directions are not computed exactly. For very large
linear programs, direct methods for solving the linear systems are often prohibitively expensive,
and iterative methods must be used. The iterative solutions, of course, contain a non-negligible
error. In [3], an efficient iterative method, namely a variant of the QMR algorithm [4] tailored to
indefinite symmetric systems, is discussed along with different strategies for preconditioning the
linear systems arising in Algorithm 1. The assumptions made for Algorithm 1 originated from and
were motivated by the method in [3].
We chose the method in [7] as a reference, since it can be adapted to inexact computations
without substantial loss in its performance, and since it is known to be efficient also for numerical
implementations; see, e.g., [9]. Other efficient interior-point methods, for example the methods
based on a self-dual formulation of the linear program, do not allow such a simple extension to
inexact computations. Apart from the specifications of the inexactness, the only modification of
the algorithm in [7] that was made for Algorithm 1 is a slightly different choice of the parameters - 1
and - 2 . These parameters are set to zero in [7]. They are also set to zero in Step 2) of Algorithm 1,
except for the case that the corresponding residual r or p is already very small. In the latter case
we do not expect much loss by not further reducing the residual at the current iteration.
In [16], another inexact infeasible-interior-point method is analyzed. The main difference lies in
the choice of the parameters - 1 and - 2 . The method in [16] uses positive values for - 1 and - 2 that
must satisfy certain conditions, and that may not be close to zero. Moreover, the algorithm in [16]
uses the value oe, or more precisely a certain norm that can be bound in terms of oe, not only for the
accuracy of the search direction, but also for determining a certain parameter ' and the step length.
Therefore, the method in [16] does not lend itself to an actual practical numerical implementation.
However, the theoretical results for the method in [16], are somewhat stronger; linear convergence
could be established, and for a special case, also a proof of polynomiality is given.


Appendix

Forms of Ill-Conditioning
When looking at standard optimization literature, we found that the definition of a condition
number of a linear program is often not explained there. Some books do talk about condition
but only for linear systems of equations. In [11, 12], a condition number of linear programs
is defined that only takes into account perturbations in the right-hand sides of the constraints. For
convex optimization problems, Renegar [19] introduced a condition-number-like distance to the
nearest ill-posed (degenerate, primal or dual having empty interior) problem. In this Appendix, we
include a precise definition of a condition number of a linear program, and explain its derivation.
We further relate the quantity oe used in our analysis to this condition number. We note that the
described condition number is different from the ones in [11, 12, 19].
The following definition is essentially taken from [20].
Let D be some open set in IR N . The condition number of a function OE
at some point x 2 D is defined as
lim sup
\Deltax!0
If both the condition number is not defined. If only the condition is
defined as +1.
Observe that this condition number still depends on the choice of a suitable norm. It is obvious,
that (in the limit, for small changes in x) this condition number is the lowest upper bound for the
quotient ' relative change of x
resulting relative change of OE(x)
For linear programs, we may define a map OE that maps the data (A; b; c) to the set of primal-dual
optimal solutions. Let D be the open set of data such that linear programs given by (A; b; c)
have a unique solution. Then OE is a function (and not a set-valued map), and the
above definition of a condition number can be applied to the linear program given by (A; b; c) . If
the set of primal-dual optimal solutions of a given linear program contains more than one point, we
define the condition number as +1. (This is a "continuous" completion of the condition number
One can also define a condition number for infeasible programs, by defining as
"solution" a certificate that the primal or dual program is infeasible. We omit this case.
We remark that for systems of linear equations, the condition number defined in the above
manner coincides with a condition number defined as the inverse of the "distance" to the nearest ill-posed
system (i.e., to the nearest singular matrix). Recall that, for convex problems, Renegar [19],
introduced a distance to the nearest ill-posed problem; it is easy to see that the inverse of this
distance does not coincide with the above condition number. Both, the above definition as well
as the definition based on the nearest ill-posed problem have properties that one would like to
associate with the condition number of the optimization problem. From the point of view of
sensitivity analysis, the above definition seems to be more appropriate, and we therefore use this
definition.
Thus, the condition number of a linear program (A; b; c) with a unique solution [x;
can be defined as
Next, we assume that (A; b; c) 2 D and briefly discuss the norm to be used.
First, we observe that for any two given positive diagonal matrices   n and   m of order n and m,
respectively, the data (A; b; c) and (  mA  n ;   m b;   n c) are equivalent in the sense that (x; y; s) is a
solution of the first problem if, and only if,
is a solution to the second problem.
Likewise, we see that b and c can be multiplied by some positive scalars - 1 and - 2 resulting in
equivalent programs by multiplying x and . The above transformations do not
encounter any cancellation error; they have condition number one. Hence, any such transformation
can be carried out almost exactly prior to solving a given linear program.
ffl Example: Consider the linear program
subject to [
with some 1. For simplicity, we ignore the dual variables. The solution is
. The 1-norm of the data 1. The above
condition number with respect to the 1-norm is at least
lim
\Deltaffl!0
This is seen when restricting oneself to equal changes \Deltaffl to both ffl-entries of the data. Note
that the last "column" of A is of norm ffl - 1. If the system is rescaled as above with
, then the rescaled version,
minimize
subject to [
has optimal solution -
perfectly well-conditioned.
It seems appropriate to eliminate the effect of such scalings from the definition of a condition of a
linear program. We may, for example, assume that the data is scaled such that b, c, and all rows
and columns of A have norm approximately 1.
As mentioned before, prior to solving the linear program, one might even compute a QR factorization
(7) of the matrix A T in a stable fashion and possibly at a higher accuracy than the standard
machine precision. If A partitioned in
a square nonsingular upper triangular matrix R 1 , then,
Given Q 1 and R \GammaT
b, the right hand side is a stable representation of the equality constraints.
After computing such a stable representation, the more difficult computations of solving the linear
program may then be performed by using Q T
rather than A. We stress again, because of the
anticipated loss of sparsity, the QR factorization is not computed in practice. We would like to
know, however, what kind of ill-conditioning of the linear program can occur when in addition to
eliminating the effects of diagonal scaling, we also assume a stable choice
1 . For this purpose
we will rewrite the linear program as follows.
Given a primal program
c T x subject to
its dual can be written as either
(D
subject to A T y - c;
or as
(D
subject to A T y
Given
1 , the matrix Q 2 is not unique in general. The dependence of Q 2 on Q 1 , however, is
well conditioned with respect to the Euclidean norm k \Delta k, in the sense that
lim sup
ae
oe
Since y is not restricted in its sign, the constraint "A T is equivalent to "Q T
Under this constraint, the dual objective value is given by b T
Thus, up to an
additive constant b
problem (D 2 ) is equivalent to
(D 0
subject to Q T
c and s - 0;
1 b. In particular, the solution s does not depend on b T R \Gamma1
c. The primal-dual
algorithm of this paper does not depend on the variable y, except for estimating the size of
the dual residual. We therefore concentrate on the changes in x and s due to small perturbations
of the data.
The set of vectors s that are feasible for (D 2 ) and for (D 0
2 ) is of course the same, and also the
angles between the rows of [ A T I ] and Q T [ A T I ] is the same. In particular,
oe min ([ A T I
We consider the condition number of (D 0
rather than the one of (D 2 ). Various sensitivities of
the solutions x and s with respect to the data of (P ) and (D 0
2 ) can be interpreted geometrically as
follows.
For the primal problem, we can identify three possible forms of ill-conditioning that we characterize
by the condition numbers C 1 , C 2 , and C 3 .
1. Assuming that the columns of A are scaled to Euclidean norm 1, we have 1 - oe
where oe max is the largest singular value of A. The condition number C 1 of computing any
point satisfying the linear equality constraints depends on oe min (A), the smallest singular value
of A. The number C 1 can be interpreted as the inverse of the smallest angle between two
linear manifolds generated by two disjoint sets of rows of the matrix A; see the angle d 1 in

Figure

1. We point out that oe min (A)=oe max (A) is not invariant under scaling of rows and
columns of A. The assumption that rows and columns have approximately constant norm is
to reduce the effect of such scaling on oe min (A). By assuming
and an ill-conditioning where C 1 is large does not occur.
2. The condition number C 2 refers to the problem of computing a certain vertex of the feasible
set, and may be interpreted as the inverse of the angle with which the affine manifold
intersects the manifold spanned by the active inequality constraints; see d 2 in Figure 1. This
condition number depends on the particular vertex, and can approximately be expressed as
the inverse of the smallest singular value of the matrix [ A T e
selection of active indices f g.
3. The condition number C 3 refers to the problem of finding which vertex is optimal in the
primal space. It depends on the angle that the objective vector c forms with the affine
manifold or with some submanifold f x
d 3 in

Figure

1. If oe min (A) is large, the condition number C 3 is linked to the inverse of the
smallest singular value of [ A T e
Program (D 0
Program (D 0
2 ) has the same structure as (P ). Tts constraint matrix Q T
, however, has perfect
condition, and thus the equivalent to C 1 vanishes. If oe min (A) is small, the dependence of Q 2 on A,
however, is not well conditioned. In fact, the equivalent to C 1 is still present in problems (D 1 ) and
d

Figure

1: Illustration of the primal problem (P )
d
s=consts
~d
Figure

2: Illustration of the dual problem (D 0
(D 2 ). The equivalent to C 2 will be denoted by C 4 , and depends on the inverse of the angles d 4 in

Figure

2, and the equivalent to C 3 will be denoted by C 5 . The corresponding angle is denoted by
d 5 in

Figure

2.



--R

Presolving in linear programming.
Solving symmetric indefinite systems in an interior-point method for linear programming

QMR: a quasi-minimal residual method for non-Hermitian linear systems
A computational view of interior-point methods for linear programming
A new polynomial-time algorithm for linear programming
A primal-dual infeasible-interior-point algorithm for linear programming
A primal-dual interior point algorithm for linear programming
Computational experience with a primal-dual interior point method for linear programming

A condition number for linear inequalities and linear programs.
Lipschitz continuity of solutions of linear inequalities
Pathways to the optimal set in linear programming.
Polynomiality of infeasible interior point algorithms for linear programming
An infeasible-interior-point algorithm using projections onto a convex set

Interior path following primal-dual algorithms
A truncated primal-infeasible dual-feasible interior point network flow method
Incorporating condition numbers into the complexity theory of linear programming.
Introduction to Numerical Analysis
Centered Newton method for mathematical programming.
Symmetric indefinite systems for interior point methods.
--TR

--CTR
Hans D. Mittelmann, Interior Point Methods for Second-Order Cone Programming and OR Applications, Computational Optimization and Applications, v.28 n.3, p.255-285, September 2004
Stefania Bellavia , Sandra Pieraccini, Convergence Analysis of an Inexact Infeasible Interior Point Method for Semidefinite Programming, Computational Optimization and Applications, v.29 n.3, p.289-313, December 2004
infeasible-interior-point method;convergence;linear program;residual;inexact search direction;linear system