convergence analysis of inexact infeasible-interior ... - Semantic Scholar
CONVERGENCE ANALYSIS OF INEXACT INFEASIBLE-INTERIOR-POINT-ALGORITHMS FOR SOLVING LINEAR PROGRAMMING PROBLEMS JANOS KORZAK
Abstract. In this paper we present a convergence analysis for some inexact (polynomial) variants of the infeasible-interior-point-algorithm of Kojima, Megiddo and Mizuno. For this analysis we assume that the iterates are bounded. The new variants allow the use of search directions that are calculated from the de ning linear system with only moderate accuracy, e.g. via the use of Krylov subspace methods like CG or QMR. Furthermore, some numerical results for the proposed methods are given. Key words. Linear programming, infeasible-interior-point method, inexact search direction AMS subject classi cations. 90C05, 65K05, 90C06
1. Introduction. Considering the fact that the primal-dual algorithm of Kojima, Megiddo and Mizuno [3] (henceforth called the KMM-Algorithm) does not use a predictor-corrector approach, it is surprising that it is an ecient algorithm for solving linear programming problems in practice. Of course, this eciency does not only follow from the use of the Newton search direction, which is clearly inferior to the predictor-corrector direction proposed by Mehrotra [7], but from the fact that the KMM-Algorithm features some other useful properties that are missing in most other infeasible-interior-point-algorithms: It allows the use of arbitrary starting points and long step sizes that can be di erent in the primal and dual subspaces. Moreover, because of the simple structure of the KMM-Algorithm, it can easily be modi ed to handle inexact search directions. In this paper we give a variant of the KMM-Algorithm that allows the use of search directions that are only calculated to moderate accuracy (inexact search directions) and prove its convergence behavior under the assumption that the iterates are bounded. This is a di erence to the analysis of the exact algorithm in [3], which gives some information about the infeasibility of the given problem if the iterates are unbounded. This (theoretical) information can not be obtained with the analysis of the present paper. The use of inexact search directions is a major di erence to most interior-pointalgorithms, whose convergence is proved under the assumption that the search directions are calculated exactly. Algorithms featuring similar search directions were proposed by Freund, Jarre and Mizuno (see [1], [2] and [8]). After some basic notes and de nitions (x2) we give a motivation for the use of inexact search directions and state our inexact variant of the KMM-Algorithm (x3). In x4 we show that the (polynomial) convergence of the new variant can be proved in almost the same way as the convergence of the original algorithm. After that we give a short analysis of the behavior of the algorithm if unsolvable problems are processed. In x5 we give a method to incorporate the predictor-corrector direction of Mehrotra in the inexact framework of this paper. Finally, we state some numerical results in x6. Throughout the paper we use the following notation: If xk and z k are elements of n , then X k and Z k denote the diagonal matrices X k = diag(xk ) and Z k = diag(z k ). IR
Fachbereich
Mathematik, Universitat Wuppertal, Gaustr. 20, D-42097 Wuppertal, Germany 1
2
JANOS KORZAK
By e we denote the vector e = (1; : : : ; 1)T 2 n , by 0 and I the zero resp. the identity matrix with sizes apparent from the context. As usual, the notation x > 0 (x 0) means that every component of x is greater then zero (nonnegative). For a matrix A 2 mn with rank(A) > 0 we denote by max (A) (min (A)) the largest (smallest positive) singular value of A, and by 2 (A) := max (A)=min (A) we denote the spectral condition number of A. 2. The Problem. In this paper we consider the linear program IR
IR
(PD) :
8 > > > >
> > > :
(x; y; z ) 2 n+m+n Ax = b AT y + z = c ; x0 z0 IR
A is an m n matrix with rank(A) = m, b 2 m and c 2 n . The set PD is called the set of all feasible points. The basic duality theorem states that (x ; y ; z ) is a IR
IR
solution of (PD), i
(x ; y ; z ) 2 f(x; y; z ) 2 n+m+n : x 0; z 0; xT z = 0; Ax = b; AT y + z = cg: IR
For given " > 0, "p > 0 and "d > 0 the infeasible-interior-point-algorithms introduced here try to calculate an element of the set
f(x; y; z ) 2
n+m+n :
x 0; z 0; xT z "; kAx ? bk2 "p ; kAT y + z ? ck2 "dg
IR
and take it as an approximation to a solution of (PD) (a so-called ("; "p ; "d)-solution). In order to ensure convergence towards an ("; "p ; "d )-solution, we force the iterates to lie within a neighborhood of the central path. In this paper we will use the following wide neighborhood proposed by Kojima, Megiddo and Mizuno [3]: Definition 2.1. Let 2 (0; 1), p > 0, d > 0, "p > 0 and "d > 0. The neighborhood N = N ( ; p ; d; "p ; "d ) is de ned by
N = f(x; y; z ) 2
n+m+n :
IR
x > 0; z > 0; xi zi xT z=n (i = 1; : : : ; n); xT z p kAx ? bk2 _ kAx ? bk2 "p ; xT z d kAT y + z ? ck2 _ kAT y + z ? ck2 "d g:
The following trivial Lemma gives a connection between N ( ; p ; d; "p ; "d) and an ("; "p ; "d)-solution: Lemma 2.2. If (x; y; z ) 2 N ( ; p ; d; "p ; "d ) and xT z minf"; "p p ; "d d g, then (x; y; z ) is an ("; "p ; "d )-solution of (PD). 3. An Inexact Infeasible-Interior-Point-Algorithm. To give a motivation for the use of inexact search directions, we review how the KMM-Algorithm tries to calculate an ("; "p ; "d)-solution. Let's assume (just for this motivation) that the set of strictly feasible points PD := f(x; y; z ) 2 PD : x > 0; z > 0g is not empty. Thus for all > 0, 0
F (x; y; z ) = @
Ax ? b AT y + z ? c Xz ? e
1 A
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
3
has a unique zero (x ; y ; z ) with x > 0 and z > 0, and lim !0(x ; y ; z ) is a solution of (PD). The KMM-Algorithm therefore calculates a decreasing sequence k with lim k = 0 and determines in each iteration k an approximation (xk+1 ; yk+1 ; z k+1 ) k!1 of the central point (xk ; yk ; zk ) via a damped Newton-step: For the given vector (xk ; yk ; z k ) 2 N the linear system 0 1 0 1 0 1 A 0 0 xk b ? Axk @ 0 AT I A @ yk A = @ c ? AT yk ? z k A (3.1) k Z 0 Xk z k k e ? X k z k is solved and kp 2 (0; 1] and kd 2 (0; 1] are chosen such that the new iterate
(xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) is an element of N (and satis es an additional descent condition). This procedure together with the fact that the search direction cannot be calculated exactly in practice, leads us to the consideration that it should be sucient to calculate an approximation of (xk ; yk ; z k ) and then proceed as stated before. In this paper we use search directions that can be calculated without any knowledge of (xk ; yk ; z k ) from (3.1): We accept (xk ; yk ; z k ) as an inexact (Newton) search direction, if 0 1 0 1 0 1 0 k 1 A 0 0 xk b ? Axk r (3.2) @ 0 AT I A @ yk A = @ c ? AT yk ? z k A + @ sk A ; Zk 0 Xk z k k e ? X k z k tk where the \residual components" satisfy 8 k (1 ? 1 )kAxk ? bk2; < kr k2 k ks k2 (1 ? 2 )kAT yk + z k ? ck2 ; (3.3)
ktk k1 3 (x )n z ; and 1 2 (0; 1], 2 2 (0; 1] and 3 2 [0; 1) are some appropriately chosen constants. :
k T k
Note that (3.3) can be interpreted as a condition of the relative exactness on each of the three components individually. Therefore, as was pointed out by one of the referees, we can not guarantee that an iterate calculated via some iterative solver applied to system (3.1) will eventually satisfy condition (3.3), even if the iteration is known to converge. However, (3.3) will be satis ed if we perform an iteration on an appropriately reduced system. This will be explained in detail at the beginning of x6, where we report on our numerical experiments. As the new \inexact" algorithm follows the central path in a less rigorous way than the \exact" KMM-Algorithm, we expect an increase in the number of iterations. But the use of inexact search directions can nevertheless result in a decrease of the total processing time, because inexact search directions can sometimes be calculated very eciently (see x6). We are now ready to state our inexact variant of the KMM-Algorithm:
Algorithm 1
1. Choose " > 0, "p > 0, "d > 0, 2 (0; 1), p > 0, d > 0, ! 1 and (x0 ; y0 ; z 0) 2 N ( ; p ; d; "p ; "d) with k(x0 ; z 0 )k1 !. Set k = 0 and choose 0 < 1 < 2 < 3 < 1, 1 2 (0; 1], 2 2 (0; 1] and 3 2 [0; 1) with (a) 1 := (1 ? ) 1 ? (1 + ) 3 > 0; (b) 2 := 1 + 1 ? 3 ? 1 > 0; (c) 3 := 1 + 2 ? 3 ? 1 > 0; (d) 4 := 2 ? 1 ? 3 > 0:
4
JANOS KORZAK
2. If (xk ; yk ; z k ) is an ("; "p ; "d)-solution, stop. 3. Set k = 1 (xk )T z k =n and calculate a search direction (xk ; yk ; z k ) that satis es (3.2) and (3.3). 4. Calculate p;k = supf 2 : xk + xk 0g; d;k = supf 2 : z k + z k 0g; ;k = minf p;k ; d;k g: If (xk ; yk ; z k ) + ;k (xk ; yk ; z k ) is an ("; "p ; "d)-solution, stop. 5. Let k be the maximum of all ~ 2 [0; 1] restricted to: All 2 [0; ~ ] are satisfying (xk ; yk ; z k ) + (xk ; yk ; z k ) 2 N ( ; p ; d ; "p ; "d); (xk + xk )T (z k + z k ) (1 ? (1 ? 2 ))(xk )T z k : 6. Choose kp 2 [0; 1] and kd 2 [0; 1] in such a way that the new iterate (xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) satis es (xk+1 ; yk+1 ; z k+1 ) 2 N ( ; p ; d ; "p ; "d); (xk+1 )T z k+1 (1 ? k (1 ? 3 ))(xk )T z k : 7. If k(xk+1 ; z k+1 )k1 !, stop. Otherwise set k = k + 1 and go to step 2. IR
IR
Remark 3.1.
1. For 1 = 1, 2 = 1 and 3 = 0 Algorithm 1 reduces to the \exact" KMMAlgorithm. 2. jj(x0 ; z 0 )jj1 is generally quite large, hence ! should be chosen to be large too. 3. The conditions (a){(d) in step 1 are met if 1 = 0:1, 2 = 0:99995, 1 = 0:95, 2 = 0:95, 3 = 0:049 and < 51=149. 4. Every vector (xk ; yk ; z k ) which sati es the conditions (3.2) and (3.3) satis es (xk ; z k ) 6 0, hence ;k 2 (0; 1). 5. Since (xk ; yk ; z k ) + k (xk ; yk ; z k ) 2 N and (xk + k xk )T (z k + k z k ) (1 ? k (1 ? 3 ))(xk )T z k ; it is always possible to choose kp = kd = k . 6. Usually it is not necessary to calculate k . The conditions in step 6 are met for a given vector (xk+1 ; yk+1 ; z k+1 ) 2 N , if (xk+1 )T z k+1 (1 ? ;k (1 ? 3 ))(xk )T z k ; because we have (by Theorem 4.2 below) (1 ? ;k (1 ? 3 ))(xk )T z k (1 ? k (1 ? 3 ))(xk )T z k : 4. Convergence Results. We now modify the ideas of Kojima, Megiddo and Mizuno [3] to prove the global convergence of Algorithm 1. We rst state some basic properties of Algorithm 1, then give a lower bound for k (Lemma 4.3), and nally prove global convergence (Theorem 4.4). We start by de ning fik ( ) = (xki + xki )(zik + zik ) ? (xk + xk )T (z k + z k )=n (i = 1; : : : ; n); gpk ( ) = (xk + xk )T (z k + z k ) ? p kA(xk + xk ) ? bk2; gdk ( ) = (xk + xk )T (z k + z k ) ? d kAT (yk + yk ) + z + z k ? ck2 ; hk ( ) = (1 ? (1 ? 2 ))(xk )T z k ? (xk + xk )T (z k + z k )
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
5
and _ k as the maximum of all ~ 2 [0; 1] for which 8 > >
> : dk h ( ) 0 hold for all 2 [0; ~ ].
The following lemma can be easily proved with the help of step 3 and step 6 of Algorithm 1, its proof is therefore omitted. Lemma 4.1. At step 2 of iteration k 1. (xk ; yk ; z k ) 2 N ( ; p ; d ; "p ; "d), 2. k(xk ; z k )k1 !, 3. (xk )T z k (1 ? k?1 (1 ? 3 ))(xk?1 )T z k?1 (if k > 0). At step 4 of iteration k we have for 2 [0; 1] 4. kA(xk + xk ) ? bk2 (1 ? 1 )kAxk ? bk2, 5. kAT (yk + yk ) + z k + z k ? ck2 (1 ? 2 )kAT yk + z k ? ck2, 6. (a) (xk + xk )T (z k + z k ) (1+ ( 1 + 3 ? 1))(xk )T z k + 2 (xk )T z k , (b) (xk + xk )T (z k + z k ) (1+ ( 1 ? 3 ? 1))(xk )T z k + 2 (xk )T z k , 7. (xki + xki )(zik + zik ) (1 ? )xki zik + ( 1 ? 3 )(xk )T z k =n + 2 xki zik for i = 1; : : : ; n. Theorem 4.2. If Algorithm 1 does not terminate at step 4 of iteration k, then k = _ k < ;k . Proof. Suppose that _ k ;k . The conditions (4.1) therefore hold for ;k . Because of the de nition of ;k , there exists an index i with (xki + ;k xki )(zik + ;k zik ) = 0, hence
fik ( ;k ) = ? (xk + ;k xk )T (z k + ;k z k )=n 0; or, equivalently, (xk + ;k xk )T (z k + ;k z k ) 0. Using (xk + ;k xk ) 0 and (z k + ;k z k ) 0 we have (xk + ;k xk )T (z k + ;k z k ) = 0. If gpk ( ;k ) < 0, then kA(xk + ;k xk ) ? bk2 "p , otherwise gpk ( ;k ) = ? pkA(xk + ;k xk ) ? bk2 0; or, equivalently, kA(xk + ;k xk ) ? bk2 = 0. In the same way we can show that kAT (yk + ;k yk )+ z k + ;k z k ? ck2 "d, hence (xk ; yk ; z k )+ ;k (xk ; yk ; z k ) is an ("; "p ; "d)-solution of (PD). This is a contradiction, because Algorithm 1 would have stopped in step 4 of iteration k. So we have shown _ k < ;k . Using the de nitions of _ k and k , one also sees that k = _ k < ;k holds. Lemma 4.3. Let > 0, jxki zik ? (xk )T z k =nj for i = 1; : : : ; n and j(xk )T z k j . At step 5 of iteration k we have
kT k kT k k ~k := min 1; 1 (xn) z ; minf2 ; 3 ; 4 g (x ) z :
Proof. Using part 7, 6(a) and 1 of Lemma 4.1 and the de nition of 1 in step 1, we have for 2 [0; ~ k ] and i = 1; : : : ; n
fik ( ) = (xki + xki )(zik + zik ) ? (xk + xk )T (z k + z k )=n
6
JANOS KORZAK
(1 ? )xki zik + ( 1 ? 3 )(xk )T z k =n + 2 xki zik ? n (1 + ( 1 + 3 ? 1))(xk )T z k ? n 2 (xk )T z k = [xki zik ? (xk )T z k =n] 2 + (1 ? )[xki zik ? (xk )T z k =n] k T k + ( 1 ? 3 ? 1 ? 3 ) (x n) z ? 2 +
kT k ((1 ? ) 1 ? (1 + ) 3 ) (x n) z
kT k = 1 (xn) z ? k )T z k ( x 1 k ? ~ n 0: If gpk (0) = (xk )T z k ? p jjAxk ? bjj2 0, we have by part 4 and 6(b) of Lemma 4.1 for all 2 [0; ~k ] gpk ( ) = (xk + xk )T (z k + z k ) ? p kA(xk + xk ) ? bk2 (1 + ( 1 ? 3 ? 1))(xk )T z k + 2 (xk )T z k ? p (1 ? 1 )kAxk ? bk2 ? 2 + (1 + (2 ? 1 ))(xk )T z k ? p (1 ? 1 )kAxk ? bk2 = ? 2 + [2 (xk )T z k ] + (1 ? 1 )gpk (0) [2 (xk )T z k ? ~k ] 0: If gpk (0) < 0, we have kAxk ? bk2 "p since (xk ; yk ; z k ) 2 N . Using part 4 of Lemma 4.1, we therefore obtain for 2 [0; ~k ] kA(xk + xk ) ? bk2 (1 ? 1 )kAxk ? bk2 (1 ? 1 )"p "p : Similiarly we can prove that for 2 [0; ~k ] 1. gdk ( ) 0 or kAT (yk + yk ) + z k + z k ? ck2 "d , 2. hk ( ) 0, hence _ k ~k . So the lemma follows from applying Theorem 4.2.
Theorem 4.4. Algorithm 1 terminates after a nite number of iterations. Proof. Suppose that Algorithm 1 does not terminate. From (xk ; yk ; z k ) 2
Lemma 2.2 and part 2 of Lemma 4.1 it follows for all k 0 (4.2)
N,
(xk )T z k " := minf"; "p p ; "d d g; k(xk ; z k )k1 ! and (xk ; z k ) 2 M := f(x; z ) 2 n+n : " n! xi ; zi ! for i = 1; : : : ; ng: IR
Using part 4 and 5 of Lemma 4.1, (xk ; z k ) 2 M , ! 1 and step 3 of Algorithm 1, we have
kb ? Axk + rk k2 kAxk ? bk2 + krk k2 kAxk ? bk2 + (1 ? 1 )kAxk ? bk2
7
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
= (2 ? 1 )kAxk ? bk2 kY ?1
(2 ? 1 )
i=0
!
(1 ? ip 1 )
kAx0 ? bk2
(2 ? 1 )kAx0 ? bk2 kc ? AT yk ? z k + sk k2 (2 ? 2 )kAT y0 + z 0 ? ck2; kk e ? X k z k + tk k2 ( 1 + 1 + 3 )pn!2 : Hence the search direction (xk ; yk ; z k ) is the solution of a linear system 0
(4.3)
@
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xk yk z k
1
0
A
=@
h i j
1 A
where
khk2 (2 ? 1 )kAx0 ? bk2; kik2 (2 ? 2 )kAT yp0 + z 0 ? ck2 ; kj k2 ( 1 + 1 + 3 ) n!2 : We note that by (xk ; z k ) 2 M and the fact that M is a compact set, there exists a
compact set K which contains the matrix of equation (4.3) for all k and every matrix B 2 K is regular. Since the inverses of all B 2 K are bounded, (xk ; yk ; z k ) is bounded for all k 0. It follows that there exists a xed > 0 with
jxki zik ? (xk )T z k =nj and j(xk )T z k j for i = 1; : : : ; n and all k 0. Using Lemma 4.3 and (4.2) we have k ~k
" " 1 := min 1; ; minf2; 3 ; 4 g 2 (0; 1];
n
hence applying part 3 of Lemma 4.1 yields (xk )T z k (1 ? (1 ? 3 ))k (x0 )T z 0 : This implies klim (xk )T z k = 0, which contradicts (4.2). !1 Remark 4.5. If ", "p , "d, , p , d , !, 1 , 2 , 3 , 4 , 1 , 2 , 3 , 1 , 2 and 3 are chosen independently of the problem, then it is possible to prove (see [4]) that the quantity in the proof of Theorem 4.4 satis es
0 ? bk2 = O 22 (AAT ) 1 + (2 ? 2 )kAT y0 + z 0 ? ck2 + (2 ? 1 )kAx max (A)
2
!
n6:5 :
Because of this relation it is easy to show that Algorithm 1 terminates after at most
0 O 22 (AAT ) 1 + (2 ? 2 )kAT y0 + z 0 ? ck2 + (2 ? 1 )kAx(A)? bk2 max
2
n7:5 L
!
8
JANOS KORZAK n
l
x) z iterations, if L is de ned as L = max 0; ln minf(";" p p ;"d d g ticular that the KMM-Algorithm terminates after at most
O
22 (AAT )
0
1 + kAT y0 + z 0 ? ck2 + kAx
T 0
mo
0 ? bk2 2
max (A)
: We have in par-
n7:5 L
!
iterations. We nish this section with some notes on the stopping criterion given in step 7 of Algorithm 1: If the algorithm terminates in step 7 of iteration j , the search direction is the exact Newton direction in iterations s to j and ! is chosen large enough, than there exists no strictly feasible point in a certain subregion of IRn (see Theorem 4.3 in [3]). Therefore a termination in step 7 can be viewed as a hint that PD = ;, or, equivalently, (see Robinson [10]) that (PD) is not solution/functional-stable. In this sense we can say that the KMM-algorithm either calculates an ("; "p ; "d)-solution or detects, that (PD) is probably unstable/unsolvable. It does not seem possible to establish a variant of this statement for inexact search directions. Therefore, if one wants to detect instability or unsolvability of a problem, we propose to restart Algorithm 1 with the current iterate as the starting point, 1 = 1, 2 = 1, 3 = 0, and a larger ! if the algorithm terminates in step 7. (One can also use an exact algorithm, e.g. the algorithm of Mizuno and Jarre [9], that guarantees the calculation of an ("; "p ; "d )-solution, if ! is large enough and (PD) is solvable.) Although this is a theoretical disadvantage, this is not important in practice: First, exact search directions can not be calculated in practice. Secondly, even for simple unsolvable problems the exact KMM-Algorithm behaves in the following way (see [4], the search directions are calculated with high accuracy and treated as exact): After a few iterations ( 10) the norm of the search direction becomes very large and the algorithm is forced to use very small step sizes to stay in N . At this point, k(xk ; z k )k1 is usually quite small ( 1020) compared with ! ( 1040 ), and k(xk ; z k )k1 is increased by only approx. 100 in each iteration. Since this means that the stopping criterion in step 7 will not be met in a reasonable time, the exact KMM-Algorithm is unable to detect the instability or the unsolvability of a problem in practice. Since Algorithm 1 is a variant of the KMM-Algorithm, it is not surprising that Algorithm 1 behaves in a similar way when being applied to unsolvable problems. It is therefore natural not to restart Algorithm 1, but to terminate in step 7 with the statement that (PD) is probably unstable or unsolvable. Moreover, because of the behavior of the norm of the search direction, it seems reasonable to use the following stopping criterion: Stop, if k(xk ; z k )k1 > !. If this stopping criterion is used by Algorithm 1, it is easy to prove (see [4]) that under the assumptions of Remark 4.5 Algorithm 1 terminates after at most O(n2 L) iterations. 5. Inexact Predictor-Corrector-Methods. We now give a variant of Algorithm 1 that allows the use of a whole class of inexact search directions. This class includes an inexact variant of the predictor-corrector search direction of Mehrotra [7], which is of one of the most ecient search directions in practice (see e.g. Lustig [5], [6]). The convergence of the given variant is ensured in a simple way: If the current search direction does not allow for suciently large step sizes, we use an inexact Newton search direction instead. After we state our algorithm, we therefore give only one remark that states the convergence of Algorithm 2.
9
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
Algorithm 2
The steps 2, 4 and 6 are identical to steps 2, 4 and 6 of Algorithm 1. 1. Choose all quantities as in step 1 of Algorithm 1 and c 2 (0; 1]. Set N 0 = 0. 3. (a) If N k = 1, calculate (xk ; yk ; z k ) as in step 3 of Algorithm 1. (b) If N k = 0, calculate a search direction (xk ; yk ; z k ) which satis es 0 @
A
0
0
1 0
0 AT I Zk 0 Xk
A @
xk yk z k
1
0
A
=@
b ? Axk c ? AT y k ? z k jk
1
0
A+@
rk sk tk
1 A
where
jjrk jj2 (1 ? 1 )jjAxk ? bjj2 ; jjsk jj2 n(1 ? 2 )jjAT ynk + z k ? cjj2 ; and j k 2 and tk 2 are arbitrary: If (xk ; z k ) 0, set N k = 1 and go to step 3. IR
IR
5. (a) This step is identical to step 5 of Algorithm 1. (b) If N k = 0 and k < c, set N k = 1 and go to step 3. 7. If k(xk+1 ; z k+1 )k1 !, stop. Otherwise set k = k + 1 and N k = 0, and go to step 2. Remark 5.1. If c is chosen independently of the problem, then Remark 4.5 can be applied to Algorithm 2. Remark 5.2. We can modify Algorithm 2 in the following way: If N k = 0, we attempt to nd a vector (xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) which satis es
(xk+1 ; yk+1 ; z k+1 ) 2 N ( ; p ; d; "p ; "d ); (xk+1 )T z k+1 (1 ? c(1 ? 3 ))(xk )T z k ; and reject the search direction only if no appropriate vector can be found through some trials (e.g. kp = 4 p;k and kd = 4 d;k , or kp = kd = 4 ;k with 4 2 (0; 1)). An inexact variant of the search direction of Mehrotra can now be incorporated in the following way: A. Calculate a predictor direction (xka ; yak ; zak ) which satis es 0 @
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xka yak zak
1
0
A
=@
b ? Axk c ? AT y k ? z k ?X k z k
1
0
A+@
where (a) jjrak jj2 (1 ? 1 )jjAxk ? bjj2 , (b) jjska jj2 (1 ? 2 )jjAT yk + z k ? cjj2 , k T k (c) jjtka jj1 3 (x )n z . k T k B. If (xk )T z k < 1, xka 0, or zak 0, de ne k = (x()n)z where (n) =
n2 ; n 5000 : n1:5 ; n > 5000
rak ska tka
1 A
10
JANOS KORZAK
Otherwise set
k ~p = min ? (xxik ) : i 2 f1; : : : ; ng and (xka )i < 0 ; ai k ~d = min ? (zzik ) : i 2 f1; : : : ; ng and (zak )i < 0 ; a i
~p = 0:99995~ p, ~d = 0:99995~ d and k =
(xk + ~p xka )T (z k + ~d zak ) (xk )T z k
!2
(xk + ~p xka )T (z k + ~d zak )
!
n
:
C. Calculate a search direction (xk ; yk ; z k ) which satis es 0 @
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xk yk z k
1
0
A
=@
b ? Axk c ? AT y k ? z k k e ? X k z k ? Xak zak
1
0
A+@
rk sk tk
1 A
where (a) jjrk jj2 (1 ? 1 )jjAxk ? bjj2 , (b) jjsk jj2 (1 ? 2 )jjAT yk + z k ? cjj2 , k T k (c) jjtk jj1 3 (x )n z . This approach has the following drawback: One of the main reasons for the eciency of the exact search direction of Mehrotra is the fact that it can be calculated with the help of only one matrix factorization. But if we determine \inexact" directions (xka ; yak ; zak ) and (xk ; yk ; z k ) with the help of Krylov subspace methods, we have to do the Krylov iteration for two linear systems. The calculation of a search direction with Krylov subspace methods can therefore be more time-consuming than the calculation via direct methods. We make some more notes on this topic in the following section. 6. Numerical Results. We now give some numerical results that were obtained with the algorithms of this paper. Inexact search directions can be calculated in several ways (see [1], [2] and [4]), but in this paper we give results for only two methods. We note that the exact Newton search direction can be calculated via (with Dk = X k (Z k )?1 and (see (4.3)) appropriate h, i and j ) (6.1)
8 < :
yk = (ADk AT )?1 (ADk (?(X k )?1 j + i) + h); z k = i ? AT yk ; xk = ?Dk (z k ? (X k )?1 j ):
For the calculation of an inexact Newton search direction we therefore de ne M = ADk AT and b = ADk (?(X k )?1 j + i) + h, and solve M yk = b with the help of A. the sparse Cholesky decomposition of the (minimum-degree reordered) matrix M, B. the preconditioned CG-algorithm. The CG-algorithm was implemented to use the simple Jacobi-preconditioner and to terminate if the CG-iterate ~yk satis es
kM ~yk ? bk2 (1 ? 1 )kAxk ? bk2:
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
11
yk is then set to ~yk , and for both methods xk and z k are then calculated as stated in (6.1). Note that in the case of Method B the vector (xk ; yk ; z k ) satis es
krk k2 = kAxk ? hk2 = kA((Z k )?1 j ? Dk z k ) ? hk2 = kA(Z k )?1 j ? ADk (i ? AT yk ) ? hk2 = kM ~yk ? bj2 (1 ? 1 )kAxk ? bk2: Furthermore, we have ksk k2 = ktk k1 = 0, hence (xk ; yk ; z k ) satis es (3.3). This
means that (xk ; yk ; z k ) is a valid inexact search direction in the sense of this paper. The search direction which is calculated with the help of Method A is usually treated as an \exact" search direction, although, due to the e ect of round-o , it is sometimes not even an inexact search direction in the sense of this paper. Nevertheless, we always treat the di erences in the results of Method A and Method B as being caused by the use of inexact search directions in Method B. Note that Algorithm 1 reduces to the KMM-Algorithm if the search direction is calculated with Method A. Algorithm 2 always tried to use the inexact predictor-corrector search direction of the proceeding section, which were calculated analogous to inexact Newton search directions. Note that Algorithm 2 reduces to a variant of the algorithm of Mehrotra (see e.g. Lustig [5], [6]) if the search direction is calculated with Method A. Before we state results, we give some details of the implementation: 1. The algorithms were programmed with MATLAB 5.3 on a Sun UltraSparc 60. 2. The starting points are calculated as proposed by Lustig [5]. 3. The parameters in step 1 of the algorithms are chosen as follows (N is enlarged, if this is necessary): "p = 1e-08, "d = 1e-08, = 1e-08, p = 1e-08,
d = 1e-08, ! = 1e+40, c = 1e-10, 1 = 0.1, 2 = 0.99995, 3 = 0.99997, 4 = 0.99995, 1 = 0.95, 2 = 0.95 and 3 = 0.049. " was set to 1e-08, 1e-07 or 5e-03. 4. In each iteration the algorithms try to use the large step sizes kp = 4 p;k and kd = 4 d;k or the step sizes kp = kd = 4 ;k . If the calculated vector does not satisfy the conditions in step 6 of the algorithms, kp = kd = k is used instead. The use of large step sizes is one of the reasons for the eciency of the exact variants of Algorithm 1 and Algorithm 2, because the use of this step sizes results in a low number of iterations. But after a few iterations it also results in a very high condition number of M , because some components of x and z become very small. This forces the CG-algorithm to use a large number of iterations for the determination of yk . We nevertheless use the stated step sizes, because we mainly want to detect di erences in the number of iterations until convergence. 5. We use the (scaled) stopping criterion used by Lustig et al. [5], that is the algorithms terminate if (xk ; yk ; z k ) satis es
(xk )T z k "; kAxk ? bk2 " ; kAT yk + z k ? ck2 " : p d 1 + jbT yk j 1 + kxk k2 1 + kz k k2 First we have a look at some results obtained with the small problem israel of the NETLIB-set. More results for 38 problems of the NETLIB-set and the problems of
12
JANOS KORZAK
the NETGEN-set can be found in [4]. The rst table gives the number of iterations and the times (in seconds) that the algorithms needed for the calculation of a (scaled) ("; "p ; "d)-solution (" was set to 1e-08).
Results for problem israel
Method A Method B iterations seconds iterations seconds Algorithm 1 35 3.60 34 31.87 Algorithm 2 23 2.56 23 27.63 We notice that the number of iterations for Method B can be compared with the number of iterations needed by Method A, but the processing time is higher by a factor of 8.9 and 10.8. We already stated that the reason for this is the high number of iterations of the CG-algorithm within each step of Algorithm 1 and Algorithm 2, which in turn is caused by the high demands on the accuracy of the search direction combined with the high condition number of M for large k. The next two pictures show the primal accuracy demand (1 ? 1 )kAxk ? bk2 and the calculation time for Algorithm 1 and Method B: 8
10
6
10
4
10
primal accuracy demand
2
10
0
10
−2
10
−4
10
−6
10
−8
10
−10
10
0
5
10
15
20
25
30
35
25
30
35
iteration
3
calculation time (seconds)
2.5
2
1.5
1
0.5
0
0
5
10
15
20 iteration
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
13
We can see that the time needed for the calculation of the search direction indeed increases for higher k. Because in the nal iterations the calculation time is 30 times as high as the calculation time in the rst iteration, and because of the fact that the calculation of the search direction via Method A needs a constant time (only 0.09 seconds), it makes sense to use Method B only in the rst few iterations and switch to Method A if Method B needs more time than Method A (this method is henceforth called Method C). Using Method C we obtain the following results: Algorithm 1 (Algorithm 2) needs 34 (23) iterations and 3.42 (2.41) seconds, with a total of only 6 (3) search directions calculated via Method B. This result (a small increase in the number of iterations, a small increase or decrease in the total processing time, and very few search directions calculated with Method B) is typical for nearly all problems of the NETLIB-set (see [4]) and can be explained by the fact that the sparse Cholesky decomposition can be calculated very eciently (in MATLAB) for those small problems. Even for the few NETLIB-problems for which Method C calculated at least two search directions with Method B the di erence in the processing time (and the number of iterations) of Method A and Method C is nearly negligible ( 7%). In order to show the performance of our algorithms on large problems, we give results for some slightly modi ed problems of the NETGEN-set (the upper bounds on x were removed). The matrices of these large problems (n = 15:325; 27:887; 52:809; m = 5:000) have a sparsity of less than 0.06%, the calculation of the Cholesky decomposition of the reordered matrix M takes 41{251 seconds. Since, on the other hand, an iteration of the CG-algorithm takes 0.01{0.08 seconds, we expect that the use of Method C will result in a substantial decrease of the needed calculation time. The following tables give the results for " = 1e-07 or 5e-03. (Smaller values for " produced numerical problems when calculating the Cholesky decomposition of M .)
Algorithm 1: Results for NETGEN-problems, " = 1e-07 Method A iterations seconds
101 102 103 105 108 109
38 40 39 37 49 54
7589 8024 7978 7863 12453 12461
Method C iterations seconds total Method B 43 38 2315 45 32 3111 46 38 2647 44 34 2979 56 30 6887 44 31 3397
speed-up
3.28 2.58 3.01 2.64 1.81 3.67
Algorithm 1: Results for NETGEN-problems, " = 5e-03 Method A iterations seconds
104 106 107 110
24 22 32 23
4932 920 7848 943
Method C iterations seconds total Method B 36 36 1163 26 21 433 42 42 1995 26 20 500
speed-up
4.24 2.08 3.93 1.89
14
JANOS KORZAK
Algorithm 2: Results for NETGEN-problems, " = 1e-07 Method A iterations seconds
20 25 23 20 28 26
101 102 103 105 108 109
Method C iterations seconds total Method B 25 21 1494 27 19 1767 25 21 1640 23 20 1492 30 23 2504 31 24 2427
4262 4891 4796 4317 7366 6341
speed-up
2.85 2.77 2.92 2.89 2.94 2.61
Algorithm 2: Results for NETGEN-problems, " = 5e-03 Method A iterations seconds
15 13 17 14
104 106 107 110
Method C iterations seconds total Method B 23 20 1081 21 16 348 22 22 628 16 14 255
3220 597 4480 617
speed-up
2.98 1.72 7.13 2.42
We notice two facts: First, the use of Method C results in an increase in the number of iterations, but because of the high number of search directions that are calculated with Method B a speed-up between 1.72 and 7.13 (usually 3) is reached. Secondly, the shortest running time is reached with Algorithm 2 and Method C. This is somewhat surprising, because in most iterations it is necessary to use the CGalgorithm for the solution of two linear systems. We nally have a look at the time that Algorithm 2 needs to calculate a search direction via Method A and Method B for the problem NETGEN 103: 3
10
Method A
2
calculation time (seconds)
10
1
10
Method B
0
10
−1
10
0
5
10
15
20
25
iteration
The plot shows that in the rst iterations Method B is much faster than Method A. Because the number of iterations of Algorithm 2 increases by only two if Method C is used instead of Method A, the use of Method C results in a huge decrease in processing
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
15
time.
7. Concluding Remarks. In this paper we proved the (polynomial) complexity of a class of inexact infeasible-interior-point-algorithms. This class includes inexact variants of some practically ecient infeasible-interior-point-algorithms, in particular variants of the algorithms of Kojima et. al. and Mehrotra. The theory developed in this paper usually justi es the use of the Cholesky decomposition for determining a search direction, because the calculated search direction, which is aicted with rounding errors, is in most cases an inexact search direction in the sense of this paper. Furthermore we have seen that the use of Krylov subspace methods results in an increase in the number of iterations, and for large problems in a huge decrease of the processing time. To make this kind of calculation more time-ecient, it seems necessary to use a more sophisticated preconditioner and to calculate the searchdirections with the help of stable linear systems if (xk )T z k approaches zero (e.g. one can try to use (M + I ) for some > 0). Another approach is to use smaller step sizes, because the used long step sizes are at least partially the reason for the high condition number of M even for small k. The use of smaller step sizes will result in an increase of the number of iterations of Algorithm 1 and Algorithm 2, but if the step sizes are chosen carefully, this can nevertheless lead to a decrease in the total processing time. Some numerical evidence for this is, that we can achieve a decrease of 22% in processing time for problem israel and Algorithm 1 with Method B, if we set 4 = 0:99 instead of 4 = 0:99995. Acknowledgments. I am grateful to the referees for suggestions and improvements upon the earlier versions of this manuscript. REFERENCES [1] Roland W. Freund, Florian Jarre, A QMR-Based Interior-Point Algorithm for Solving Linear Programs, Numerical Analysis Manuscript 94-19, AT&T Bell Laboratories, December 1994 [2] Roland W. Freund, Florian Jarre, Shinji Mizuno, Convergence of a Class of Inexact Interior-Point Algorithms for Linear Programs, Interior-Point Methods Online, Argonne National Laboratory, http://www.mcs.anl.gov/home/otc/InteriorPoint/, September 1996 [3] Masakazu Kojima, Nimrod Megiddo, Shinji Mizuno, A primal-dual infeasible-interior-point algorithm for linear programming, Mathematical Programming, Vol. 61, pp. 263{280, 1993 [4] Janos Korzak, Primal-Duale pfadfolgende inexakte Auere-Punkte-Verfahren zur Losung linearer Optimierungsaufgaben, Dissertation, Fachbereich Mathematik, Bergische Universitat - Gesamthochschule Wuppertal, Januar 1997 [5] Irvin J. Lustig, Roy E. Marsten, David F. Shanno, On Implementing Mehrotra's PredictorCorrector Interior Point Method for Linear Programming, SIAM Journal on Optimization, Vol. 2, No. 3, pp. 435{449, August 1992 , Computational experience with a globally convergent primal-dual predictor-corrector [6] algorithm for linear programming, Mathematical Programming, Vol. 66, pp. 123{135, 1994 [7] Sanjay Mehrotra, On the Implementation of a Primal-Dual Interior Point Method, SIAM Journal on Optimization, Vol. 2, No. 4, November 1992, pp. 575{601 [8] Shinji Mizuno, Florian Jarre, Global and Polynomial-Time Convergence of an InfeasibleInterior-Point Algorithm Using Inexact Computation, Research Memorandum 605, The Institute of Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 106, Japan, April 1996 [9] , An Infeasible-Interior-Point Algorithm Using Projections onto a Convex Set, Annals of Opererations Research, Vol. 62, pp. 59-80, 1996 [10] Stephan M. Robinson, A Characterization of Stability in Linear Programming, Operations Research, Vol. 25, No. 3, May-June 1977
Abstract. In this paper we present a convergence analysis for some inexact (polynomial) variants of the infeasible-interior-point-algorithm of Kojima, Megiddo and Mizuno. For this analysis we assume that the iterates are bounded. The new variants allow the use of search directions that are calculated from the de ning linear system with only moderate accuracy, e.g. via the use of Krylov subspace methods like CG or QMR. Furthermore, some numerical results for the proposed methods are given. Key words. Linear programming, infeasible-interior-point method, inexact search direction AMS subject classi cations. 90C05, 65K05, 90C06
1. Introduction. Considering the fact that the primal-dual algorithm of Kojima, Megiddo and Mizuno [3] (henceforth called the KMM-Algorithm) does not use a predictor-corrector approach, it is surprising that it is an ecient algorithm for solving linear programming problems in practice. Of course, this eciency does not only follow from the use of the Newton search direction, which is clearly inferior to the predictor-corrector direction proposed by Mehrotra [7], but from the fact that the KMM-Algorithm features some other useful properties that are missing in most other infeasible-interior-point-algorithms: It allows the use of arbitrary starting points and long step sizes that can be di erent in the primal and dual subspaces. Moreover, because of the simple structure of the KMM-Algorithm, it can easily be modi ed to handle inexact search directions. In this paper we give a variant of the KMM-Algorithm that allows the use of search directions that are only calculated to moderate accuracy (inexact search directions) and prove its convergence behavior under the assumption that the iterates are bounded. This is a di erence to the analysis of the exact algorithm in [3], which gives some information about the infeasibility of the given problem if the iterates are unbounded. This (theoretical) information can not be obtained with the analysis of the present paper. The use of inexact search directions is a major di erence to most interior-pointalgorithms, whose convergence is proved under the assumption that the search directions are calculated exactly. Algorithms featuring similar search directions were proposed by Freund, Jarre and Mizuno (see [1], [2] and [8]). After some basic notes and de nitions (x2) we give a motivation for the use of inexact search directions and state our inexact variant of the KMM-Algorithm (x3). In x4 we show that the (polynomial) convergence of the new variant can be proved in almost the same way as the convergence of the original algorithm. After that we give a short analysis of the behavior of the algorithm if unsolvable problems are processed. In x5 we give a method to incorporate the predictor-corrector direction of Mehrotra in the inexact framework of this paper. Finally, we state some numerical results in x6. Throughout the paper we use the following notation: If xk and z k are elements of n , then X k and Z k denote the diagonal matrices X k = diag(xk ) and Z k = diag(z k ). IR
Fachbereich
Mathematik, Universitat Wuppertal, Gaustr. 20, D-42097 Wuppertal, Germany 1
2
JANOS KORZAK
By e we denote the vector e = (1; : : : ; 1)T 2 n , by 0 and I the zero resp. the identity matrix with sizes apparent from the context. As usual, the notation x > 0 (x 0) means that every component of x is greater then zero (nonnegative). For a matrix A 2 mn with rank(A) > 0 we denote by max (A) (min (A)) the largest (smallest positive) singular value of A, and by 2 (A) := max (A)=min (A) we denote the spectral condition number of A. 2. The Problem. In this paper we consider the linear program IR
IR
(PD) :
8 > > > >
> > > :
(x; y; z ) 2 n+m+n Ax = b AT y + z = c ; x0 z0 IR
A is an m n matrix with rank(A) = m, b 2 m and c 2 n . The set PD is called the set of all feasible points. The basic duality theorem states that (x ; y ; z ) is a IR
IR
solution of (PD), i
(x ; y ; z ) 2 f(x; y; z ) 2 n+m+n : x 0; z 0; xT z = 0; Ax = b; AT y + z = cg: IR
For given " > 0, "p > 0 and "d > 0 the infeasible-interior-point-algorithms introduced here try to calculate an element of the set
f(x; y; z ) 2
n+m+n :
x 0; z 0; xT z "; kAx ? bk2 "p ; kAT y + z ? ck2 "dg
IR
and take it as an approximation to a solution of (PD) (a so-called ("; "p ; "d)-solution). In order to ensure convergence towards an ("; "p ; "d )-solution, we force the iterates to lie within a neighborhood of the central path. In this paper we will use the following wide neighborhood proposed by Kojima, Megiddo and Mizuno [3]: Definition 2.1. Let 2 (0; 1), p > 0, d > 0, "p > 0 and "d > 0. The neighborhood N = N ( ; p ; d; "p ; "d ) is de ned by
N = f(x; y; z ) 2
n+m+n :
IR
x > 0; z > 0; xi zi xT z=n (i = 1; : : : ; n); xT z p kAx ? bk2 _ kAx ? bk2 "p ; xT z d kAT y + z ? ck2 _ kAT y + z ? ck2 "d g:
The following trivial Lemma gives a connection between N ( ; p ; d; "p ; "d) and an ("; "p ; "d)-solution: Lemma 2.2. If (x; y; z ) 2 N ( ; p ; d; "p ; "d ) and xT z minf"; "p p ; "d d g, then (x; y; z ) is an ("; "p ; "d )-solution of (PD). 3. An Inexact Infeasible-Interior-Point-Algorithm. To give a motivation for the use of inexact search directions, we review how the KMM-Algorithm tries to calculate an ("; "p ; "d)-solution. Let's assume (just for this motivation) that the set of strictly feasible points PD := f(x; y; z ) 2 PD : x > 0; z > 0g is not empty. Thus for all > 0, 0
F (x; y; z ) = @
Ax ? b AT y + z ? c Xz ? e
1 A
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
3
has a unique zero (x ; y ; z ) with x > 0 and z > 0, and lim !0(x ; y ; z ) is a solution of (PD). The KMM-Algorithm therefore calculates a decreasing sequence k with lim k = 0 and determines in each iteration k an approximation (xk+1 ; yk+1 ; z k+1 ) k!1 of the central point (xk ; yk ; zk ) via a damped Newton-step: For the given vector (xk ; yk ; z k ) 2 N the linear system 0 1 0 1 0 1 A 0 0 xk b ? Axk @ 0 AT I A @ yk A = @ c ? AT yk ? z k A (3.1) k Z 0 Xk z k k e ? X k z k is solved and kp 2 (0; 1] and kd 2 (0; 1] are chosen such that the new iterate
(xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) is an element of N (and satis es an additional descent condition). This procedure together with the fact that the search direction cannot be calculated exactly in practice, leads us to the consideration that it should be sucient to calculate an approximation of (xk ; yk ; z k ) and then proceed as stated before. In this paper we use search directions that can be calculated without any knowledge of (xk ; yk ; z k ) from (3.1): We accept (xk ; yk ; z k ) as an inexact (Newton) search direction, if 0 1 0 1 0 1 0 k 1 A 0 0 xk b ? Axk r (3.2) @ 0 AT I A @ yk A = @ c ? AT yk ? z k A + @ sk A ; Zk 0 Xk z k k e ? X k z k tk where the \residual components" satisfy 8 k (1 ? 1 )kAxk ? bk2; < kr k2 k ks k2 (1 ? 2 )kAT yk + z k ? ck2 ; (3.3)
ktk k1 3 (x )n z ; and 1 2 (0; 1], 2 2 (0; 1] and 3 2 [0; 1) are some appropriately chosen constants. :
k T k
Note that (3.3) can be interpreted as a condition of the relative exactness on each of the three components individually. Therefore, as was pointed out by one of the referees, we can not guarantee that an iterate calculated via some iterative solver applied to system (3.1) will eventually satisfy condition (3.3), even if the iteration is known to converge. However, (3.3) will be satis ed if we perform an iteration on an appropriately reduced system. This will be explained in detail at the beginning of x6, where we report on our numerical experiments. As the new \inexact" algorithm follows the central path in a less rigorous way than the \exact" KMM-Algorithm, we expect an increase in the number of iterations. But the use of inexact search directions can nevertheless result in a decrease of the total processing time, because inexact search directions can sometimes be calculated very eciently (see x6). We are now ready to state our inexact variant of the KMM-Algorithm:
Algorithm 1
1. Choose " > 0, "p > 0, "d > 0, 2 (0; 1), p > 0, d > 0, ! 1 and (x0 ; y0 ; z 0) 2 N ( ; p ; d; "p ; "d) with k(x0 ; z 0 )k1 !. Set k = 0 and choose 0 < 1 < 2 < 3 < 1, 1 2 (0; 1], 2 2 (0; 1] and 3 2 [0; 1) with (a) 1 := (1 ? ) 1 ? (1 + ) 3 > 0; (b) 2 := 1 + 1 ? 3 ? 1 > 0; (c) 3 := 1 + 2 ? 3 ? 1 > 0; (d) 4 := 2 ? 1 ? 3 > 0:
4
JANOS KORZAK
2. If (xk ; yk ; z k ) is an ("; "p ; "d)-solution, stop. 3. Set k = 1 (xk )T z k =n and calculate a search direction (xk ; yk ; z k ) that satis es (3.2) and (3.3). 4. Calculate p;k = supf 2 : xk + xk 0g; d;k = supf 2 : z k + z k 0g; ;k = minf p;k ; d;k g: If (xk ; yk ; z k ) + ;k (xk ; yk ; z k ) is an ("; "p ; "d)-solution, stop. 5. Let k be the maximum of all ~ 2 [0; 1] restricted to: All 2 [0; ~ ] are satisfying (xk ; yk ; z k ) + (xk ; yk ; z k ) 2 N ( ; p ; d ; "p ; "d); (xk + xk )T (z k + z k ) (1 ? (1 ? 2 ))(xk )T z k : 6. Choose kp 2 [0; 1] and kd 2 [0; 1] in such a way that the new iterate (xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) satis es (xk+1 ; yk+1 ; z k+1 ) 2 N ( ; p ; d ; "p ; "d); (xk+1 )T z k+1 (1 ? k (1 ? 3 ))(xk )T z k : 7. If k(xk+1 ; z k+1 )k1 !, stop. Otherwise set k = k + 1 and go to step 2. IR
IR
Remark 3.1.
1. For 1 = 1, 2 = 1 and 3 = 0 Algorithm 1 reduces to the \exact" KMMAlgorithm. 2. jj(x0 ; z 0 )jj1 is generally quite large, hence ! should be chosen to be large too. 3. The conditions (a){(d) in step 1 are met if 1 = 0:1, 2 = 0:99995, 1 = 0:95, 2 = 0:95, 3 = 0:049 and < 51=149. 4. Every vector (xk ; yk ; z k ) which sati es the conditions (3.2) and (3.3) satis es (xk ; z k ) 6 0, hence ;k 2 (0; 1). 5. Since (xk ; yk ; z k ) + k (xk ; yk ; z k ) 2 N and (xk + k xk )T (z k + k z k ) (1 ? k (1 ? 3 ))(xk )T z k ; it is always possible to choose kp = kd = k . 6. Usually it is not necessary to calculate k . The conditions in step 6 are met for a given vector (xk+1 ; yk+1 ; z k+1 ) 2 N , if (xk+1 )T z k+1 (1 ? ;k (1 ? 3 ))(xk )T z k ; because we have (by Theorem 4.2 below) (1 ? ;k (1 ? 3 ))(xk )T z k (1 ? k (1 ? 3 ))(xk )T z k : 4. Convergence Results. We now modify the ideas of Kojima, Megiddo and Mizuno [3] to prove the global convergence of Algorithm 1. We rst state some basic properties of Algorithm 1, then give a lower bound for k (Lemma 4.3), and nally prove global convergence (Theorem 4.4). We start by de ning fik ( ) = (xki + xki )(zik + zik ) ? (xk + xk )T (z k + z k )=n (i = 1; : : : ; n); gpk ( ) = (xk + xk )T (z k + z k ) ? p kA(xk + xk ) ? bk2; gdk ( ) = (xk + xk )T (z k + z k ) ? d kAT (yk + yk ) + z + z k ? ck2 ; hk ( ) = (1 ? (1 ? 2 ))(xk )T z k ? (xk + xk )T (z k + z k )
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
5
and _ k as the maximum of all ~ 2 [0; 1] for which 8 > >
> : dk h ( ) 0 hold for all 2 [0; ~ ].
The following lemma can be easily proved with the help of step 3 and step 6 of Algorithm 1, its proof is therefore omitted. Lemma 4.1. At step 2 of iteration k 1. (xk ; yk ; z k ) 2 N ( ; p ; d ; "p ; "d), 2. k(xk ; z k )k1 !, 3. (xk )T z k (1 ? k?1 (1 ? 3 ))(xk?1 )T z k?1 (if k > 0). At step 4 of iteration k we have for 2 [0; 1] 4. kA(xk + xk ) ? bk2 (1 ? 1 )kAxk ? bk2, 5. kAT (yk + yk ) + z k + z k ? ck2 (1 ? 2 )kAT yk + z k ? ck2, 6. (a) (xk + xk )T (z k + z k ) (1+ ( 1 + 3 ? 1))(xk )T z k + 2 (xk )T z k , (b) (xk + xk )T (z k + z k ) (1+ ( 1 ? 3 ? 1))(xk )T z k + 2 (xk )T z k , 7. (xki + xki )(zik + zik ) (1 ? )xki zik + ( 1 ? 3 )(xk )T z k =n + 2 xki zik for i = 1; : : : ; n. Theorem 4.2. If Algorithm 1 does not terminate at step 4 of iteration k, then k = _ k < ;k . Proof. Suppose that _ k ;k . The conditions (4.1) therefore hold for ;k . Because of the de nition of ;k , there exists an index i with (xki + ;k xki )(zik + ;k zik ) = 0, hence
fik ( ;k ) = ? (xk + ;k xk )T (z k + ;k z k )=n 0; or, equivalently, (xk + ;k xk )T (z k + ;k z k ) 0. Using (xk + ;k xk ) 0 and (z k + ;k z k ) 0 we have (xk + ;k xk )T (z k + ;k z k ) = 0. If gpk ( ;k ) < 0, then kA(xk + ;k xk ) ? bk2 "p , otherwise gpk ( ;k ) = ? pkA(xk + ;k xk ) ? bk2 0; or, equivalently, kA(xk + ;k xk ) ? bk2 = 0. In the same way we can show that kAT (yk + ;k yk )+ z k + ;k z k ? ck2 "d, hence (xk ; yk ; z k )+ ;k (xk ; yk ; z k ) is an ("; "p ; "d)-solution of (PD). This is a contradiction, because Algorithm 1 would have stopped in step 4 of iteration k. So we have shown _ k < ;k . Using the de nitions of _ k and k , one also sees that k = _ k < ;k holds. Lemma 4.3. Let > 0, jxki zik ? (xk )T z k =nj for i = 1; : : : ; n and j(xk )T z k j . At step 5 of iteration k we have
kT k kT k k ~k := min 1; 1 (xn) z ; minf2 ; 3 ; 4 g (x ) z :
Proof. Using part 7, 6(a) and 1 of Lemma 4.1 and the de nition of 1 in step 1, we have for 2 [0; ~ k ] and i = 1; : : : ; n
fik ( ) = (xki + xki )(zik + zik ) ? (xk + xk )T (z k + z k )=n
6
JANOS KORZAK
(1 ? )xki zik + ( 1 ? 3 )(xk )T z k =n + 2 xki zik ? n (1 + ( 1 + 3 ? 1))(xk )T z k ? n 2 (xk )T z k = [xki zik ? (xk )T z k =n] 2 + (1 ? )[xki zik ? (xk )T z k =n] k T k + ( 1 ? 3 ? 1 ? 3 ) (x n) z ? 2 +
kT k ((1 ? ) 1 ? (1 + ) 3 ) (x n) z
kT k = 1 (xn) z ? k )T z k ( x 1 k ? ~ n 0: If gpk (0) = (xk )T z k ? p jjAxk ? bjj2 0, we have by part 4 and 6(b) of Lemma 4.1 for all 2 [0; ~k ] gpk ( ) = (xk + xk )T (z k + z k ) ? p kA(xk + xk ) ? bk2 (1 + ( 1 ? 3 ? 1))(xk )T z k + 2 (xk )T z k ? p (1 ? 1 )kAxk ? bk2 ? 2 + (1 + (2 ? 1 ))(xk )T z k ? p (1 ? 1 )kAxk ? bk2 = ? 2 + [2 (xk )T z k ] + (1 ? 1 )gpk (0) [2 (xk )T z k ? ~k ] 0: If gpk (0) < 0, we have kAxk ? bk2 "p since (xk ; yk ; z k ) 2 N . Using part 4 of Lemma 4.1, we therefore obtain for 2 [0; ~k ] kA(xk + xk ) ? bk2 (1 ? 1 )kAxk ? bk2 (1 ? 1 )"p "p : Similiarly we can prove that for 2 [0; ~k ] 1. gdk ( ) 0 or kAT (yk + yk ) + z k + z k ? ck2 "d , 2. hk ( ) 0, hence _ k ~k . So the lemma follows from applying Theorem 4.2.
Theorem 4.4. Algorithm 1 terminates after a nite number of iterations. Proof. Suppose that Algorithm 1 does not terminate. From (xk ; yk ; z k ) 2
Lemma 2.2 and part 2 of Lemma 4.1 it follows for all k 0 (4.2)
N,
(xk )T z k " := minf"; "p p ; "d d g; k(xk ; z k )k1 ! and (xk ; z k ) 2 M := f(x; z ) 2 n+n : " n! xi ; zi ! for i = 1; : : : ; ng: IR
Using part 4 and 5 of Lemma 4.1, (xk ; z k ) 2 M , ! 1 and step 3 of Algorithm 1, we have
kb ? Axk + rk k2 kAxk ? bk2 + krk k2 kAxk ? bk2 + (1 ? 1 )kAxk ? bk2
7
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
= (2 ? 1 )kAxk ? bk2 kY ?1
(2 ? 1 )
i=0
!
(1 ? ip 1 )
kAx0 ? bk2
(2 ? 1 )kAx0 ? bk2 kc ? AT yk ? z k + sk k2 (2 ? 2 )kAT y0 + z 0 ? ck2; kk e ? X k z k + tk k2 ( 1 + 1 + 3 )pn!2 : Hence the search direction (xk ; yk ; z k ) is the solution of a linear system 0
(4.3)
@
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xk yk z k
1
0
A
=@
h i j
1 A
where
khk2 (2 ? 1 )kAx0 ? bk2; kik2 (2 ? 2 )kAT yp0 + z 0 ? ck2 ; kj k2 ( 1 + 1 + 3 ) n!2 : We note that by (xk ; z k ) 2 M and the fact that M is a compact set, there exists a
compact set K which contains the matrix of equation (4.3) for all k and every matrix B 2 K is regular. Since the inverses of all B 2 K are bounded, (xk ; yk ; z k ) is bounded for all k 0. It follows that there exists a xed > 0 with
jxki zik ? (xk )T z k =nj and j(xk )T z k j for i = 1; : : : ; n and all k 0. Using Lemma 4.3 and (4.2) we have k ~k
" " 1 := min 1; ; minf2; 3 ; 4 g 2 (0; 1];
n
hence applying part 3 of Lemma 4.1 yields (xk )T z k (1 ? (1 ? 3 ))k (x0 )T z 0 : This implies klim (xk )T z k = 0, which contradicts (4.2). !1 Remark 4.5. If ", "p , "d, , p , d , !, 1 , 2 , 3 , 4 , 1 , 2 , 3 , 1 , 2 and 3 are chosen independently of the problem, then it is possible to prove (see [4]) that the quantity in the proof of Theorem 4.4 satis es
0 ? bk2 = O 22 (AAT ) 1 + (2 ? 2 )kAT y0 + z 0 ? ck2 + (2 ? 1 )kAx max (A)
2
!
n6:5 :
Because of this relation it is easy to show that Algorithm 1 terminates after at most
0 O 22 (AAT ) 1 + (2 ? 2 )kAT y0 + z 0 ? ck2 + (2 ? 1 )kAx(A)? bk2 max
2
n7:5 L
!
8
JANOS KORZAK n
l
x) z iterations, if L is de ned as L = max 0; ln minf(";" p p ;"d d g ticular that the KMM-Algorithm terminates after at most
O
22 (AAT )
0
1 + kAT y0 + z 0 ? ck2 + kAx
T 0
mo
0 ? bk2 2
max (A)
: We have in par-
n7:5 L
!
iterations. We nish this section with some notes on the stopping criterion given in step 7 of Algorithm 1: If the algorithm terminates in step 7 of iteration j , the search direction is the exact Newton direction in iterations s to j and ! is chosen large enough, than there exists no strictly feasible point in a certain subregion of IRn (see Theorem 4.3 in [3]). Therefore a termination in step 7 can be viewed as a hint that PD = ;, or, equivalently, (see Robinson [10]) that (PD) is not solution/functional-stable. In this sense we can say that the KMM-algorithm either calculates an ("; "p ; "d)-solution or detects, that (PD) is probably unstable/unsolvable. It does not seem possible to establish a variant of this statement for inexact search directions. Therefore, if one wants to detect instability or unsolvability of a problem, we propose to restart Algorithm 1 with the current iterate as the starting point, 1 = 1, 2 = 1, 3 = 0, and a larger ! if the algorithm terminates in step 7. (One can also use an exact algorithm, e.g. the algorithm of Mizuno and Jarre [9], that guarantees the calculation of an ("; "p ; "d )-solution, if ! is large enough and (PD) is solvable.) Although this is a theoretical disadvantage, this is not important in practice: First, exact search directions can not be calculated in practice. Secondly, even for simple unsolvable problems the exact KMM-Algorithm behaves in the following way (see [4], the search directions are calculated with high accuracy and treated as exact): After a few iterations ( 10) the norm of the search direction becomes very large and the algorithm is forced to use very small step sizes to stay in N . At this point, k(xk ; z k )k1 is usually quite small ( 1020) compared with ! ( 1040 ), and k(xk ; z k )k1 is increased by only approx. 100 in each iteration. Since this means that the stopping criterion in step 7 will not be met in a reasonable time, the exact KMM-Algorithm is unable to detect the instability or the unsolvability of a problem in practice. Since Algorithm 1 is a variant of the KMM-Algorithm, it is not surprising that Algorithm 1 behaves in a similar way when being applied to unsolvable problems. It is therefore natural not to restart Algorithm 1, but to terminate in step 7 with the statement that (PD) is probably unstable or unsolvable. Moreover, because of the behavior of the norm of the search direction, it seems reasonable to use the following stopping criterion: Stop, if k(xk ; z k )k1 > !. If this stopping criterion is used by Algorithm 1, it is easy to prove (see [4]) that under the assumptions of Remark 4.5 Algorithm 1 terminates after at most O(n2 L) iterations. 5. Inexact Predictor-Corrector-Methods. We now give a variant of Algorithm 1 that allows the use of a whole class of inexact search directions. This class includes an inexact variant of the predictor-corrector search direction of Mehrotra [7], which is of one of the most ecient search directions in practice (see e.g. Lustig [5], [6]). The convergence of the given variant is ensured in a simple way: If the current search direction does not allow for suciently large step sizes, we use an inexact Newton search direction instead. After we state our algorithm, we therefore give only one remark that states the convergence of Algorithm 2.
9
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
Algorithm 2
The steps 2, 4 and 6 are identical to steps 2, 4 and 6 of Algorithm 1. 1. Choose all quantities as in step 1 of Algorithm 1 and c 2 (0; 1]. Set N 0 = 0. 3. (a) If N k = 1, calculate (xk ; yk ; z k ) as in step 3 of Algorithm 1. (b) If N k = 0, calculate a search direction (xk ; yk ; z k ) which satis es 0 @
A
0
0
1 0
0 AT I Zk 0 Xk
A @
xk yk z k
1
0
A
=@
b ? Axk c ? AT y k ? z k jk
1
0
A+@
rk sk tk
1 A
where
jjrk jj2 (1 ? 1 )jjAxk ? bjj2 ; jjsk jj2 n(1 ? 2 )jjAT ynk + z k ? cjj2 ; and j k 2 and tk 2 are arbitrary: If (xk ; z k ) 0, set N k = 1 and go to step 3. IR
IR
5. (a) This step is identical to step 5 of Algorithm 1. (b) If N k = 0 and k < c, set N k = 1 and go to step 3. 7. If k(xk+1 ; z k+1 )k1 !, stop. Otherwise set k = k + 1 and N k = 0, and go to step 2. Remark 5.1. If c is chosen independently of the problem, then Remark 4.5 can be applied to Algorithm 2. Remark 5.2. We can modify Algorithm 2 in the following way: If N k = 0, we attempt to nd a vector (xk+1 ; yk+1 ; z k+1 ) = (xk + kp xk ; yk + kd yk ; z k + kd z k ) which satis es
(xk+1 ; yk+1 ; z k+1 ) 2 N ( ; p ; d; "p ; "d ); (xk+1 )T z k+1 (1 ? c(1 ? 3 ))(xk )T z k ; and reject the search direction only if no appropriate vector can be found through some trials (e.g. kp = 4 p;k and kd = 4 d;k , or kp = kd = 4 ;k with 4 2 (0; 1)). An inexact variant of the search direction of Mehrotra can now be incorporated in the following way: A. Calculate a predictor direction (xka ; yak ; zak ) which satis es 0 @
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xka yak zak
1
0
A
=@
b ? Axk c ? AT y k ? z k ?X k z k
1
0
A+@
where (a) jjrak jj2 (1 ? 1 )jjAxk ? bjj2 , (b) jjska jj2 (1 ? 2 )jjAT yk + z k ? cjj2 , k T k (c) jjtka jj1 3 (x )n z . k T k B. If (xk )T z k < 1, xka 0, or zak 0, de ne k = (x()n)z where (n) =
n2 ; n 5000 : n1:5 ; n > 5000
rak ska tka
1 A
10
JANOS KORZAK
Otherwise set
k ~p = min ? (xxik ) : i 2 f1; : : : ; ng and (xka )i < 0 ; ai k ~d = min ? (zzik ) : i 2 f1; : : : ; ng and (zak )i < 0 ; a i
~p = 0:99995~ p, ~d = 0:99995~ d and k =
(xk + ~p xka )T (z k + ~d zak ) (xk )T z k
!2
(xk + ~p xka )T (z k + ~d zak )
!
n
:
C. Calculate a search direction (xk ; yk ; z k ) which satis es 0 @
A
0
0
0 AT I Zk 0 Xk
1 0
A @
xk yk z k
1
0
A
=@
b ? Axk c ? AT y k ? z k k e ? X k z k ? Xak zak
1
0
A+@
rk sk tk
1 A
where (a) jjrk jj2 (1 ? 1 )jjAxk ? bjj2 , (b) jjsk jj2 (1 ? 2 )jjAT yk + z k ? cjj2 , k T k (c) jjtk jj1 3 (x )n z . This approach has the following drawback: One of the main reasons for the eciency of the exact search direction of Mehrotra is the fact that it can be calculated with the help of only one matrix factorization. But if we determine \inexact" directions (xka ; yak ; zak ) and (xk ; yk ; z k ) with the help of Krylov subspace methods, we have to do the Krylov iteration for two linear systems. The calculation of a search direction with Krylov subspace methods can therefore be more time-consuming than the calculation via direct methods. We make some more notes on this topic in the following section. 6. Numerical Results. We now give some numerical results that were obtained with the algorithms of this paper. Inexact search directions can be calculated in several ways (see [1], [2] and [4]), but in this paper we give results for only two methods. We note that the exact Newton search direction can be calculated via (with Dk = X k (Z k )?1 and (see (4.3)) appropriate h, i and j ) (6.1)
8 < :
yk = (ADk AT )?1 (ADk (?(X k )?1 j + i) + h); z k = i ? AT yk ; xk = ?Dk (z k ? (X k )?1 j ):
For the calculation of an inexact Newton search direction we therefore de ne M = ADk AT and b = ADk (?(X k )?1 j + i) + h, and solve M yk = b with the help of A. the sparse Cholesky decomposition of the (minimum-degree reordered) matrix M, B. the preconditioned CG-algorithm. The CG-algorithm was implemented to use the simple Jacobi-preconditioner and to terminate if the CG-iterate ~yk satis es
kM ~yk ? bk2 (1 ? 1 )kAxk ? bk2:
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
11
yk is then set to ~yk , and for both methods xk and z k are then calculated as stated in (6.1). Note that in the case of Method B the vector (xk ; yk ; z k ) satis es
krk k2 = kAxk ? hk2 = kA((Z k )?1 j ? Dk z k ) ? hk2 = kA(Z k )?1 j ? ADk (i ? AT yk ) ? hk2 = kM ~yk ? bj2 (1 ? 1 )kAxk ? bk2: Furthermore, we have ksk k2 = ktk k1 = 0, hence (xk ; yk ; z k ) satis es (3.3). This
means that (xk ; yk ; z k ) is a valid inexact search direction in the sense of this paper. The search direction which is calculated with the help of Method A is usually treated as an \exact" search direction, although, due to the e ect of round-o , it is sometimes not even an inexact search direction in the sense of this paper. Nevertheless, we always treat the di erences in the results of Method A and Method B as being caused by the use of inexact search directions in Method B. Note that Algorithm 1 reduces to the KMM-Algorithm if the search direction is calculated with Method A. Algorithm 2 always tried to use the inexact predictor-corrector search direction of the proceeding section, which were calculated analogous to inexact Newton search directions. Note that Algorithm 2 reduces to a variant of the algorithm of Mehrotra (see e.g. Lustig [5], [6]) if the search direction is calculated with Method A. Before we state results, we give some details of the implementation: 1. The algorithms were programmed with MATLAB 5.3 on a Sun UltraSparc 60. 2. The starting points are calculated as proposed by Lustig [5]. 3. The parameters in step 1 of the algorithms are chosen as follows (N is enlarged, if this is necessary): "p = 1e-08, "d = 1e-08, = 1e-08, p = 1e-08,
d = 1e-08, ! = 1e+40, c = 1e-10, 1 = 0.1, 2 = 0.99995, 3 = 0.99997, 4 = 0.99995, 1 = 0.95, 2 = 0.95 and 3 = 0.049. " was set to 1e-08, 1e-07 or 5e-03. 4. In each iteration the algorithms try to use the large step sizes kp = 4 p;k and kd = 4 d;k or the step sizes kp = kd = 4 ;k . If the calculated vector does not satisfy the conditions in step 6 of the algorithms, kp = kd = k is used instead. The use of large step sizes is one of the reasons for the eciency of the exact variants of Algorithm 1 and Algorithm 2, because the use of this step sizes results in a low number of iterations. But after a few iterations it also results in a very high condition number of M , because some components of x and z become very small. This forces the CG-algorithm to use a large number of iterations for the determination of yk . We nevertheless use the stated step sizes, because we mainly want to detect di erences in the number of iterations until convergence. 5. We use the (scaled) stopping criterion used by Lustig et al. [5], that is the algorithms terminate if (xk ; yk ; z k ) satis es
(xk )T z k "; kAxk ? bk2 " ; kAT yk + z k ? ck2 " : p d 1 + jbT yk j 1 + kxk k2 1 + kz k k2 First we have a look at some results obtained with the small problem israel of the NETLIB-set. More results for 38 problems of the NETLIB-set and the problems of
12
JANOS KORZAK
the NETGEN-set can be found in [4]. The rst table gives the number of iterations and the times (in seconds) that the algorithms needed for the calculation of a (scaled) ("; "p ; "d)-solution (" was set to 1e-08).
Results for problem israel
Method A Method B iterations seconds iterations seconds Algorithm 1 35 3.60 34 31.87 Algorithm 2 23 2.56 23 27.63 We notice that the number of iterations for Method B can be compared with the number of iterations needed by Method A, but the processing time is higher by a factor of 8.9 and 10.8. We already stated that the reason for this is the high number of iterations of the CG-algorithm within each step of Algorithm 1 and Algorithm 2, which in turn is caused by the high demands on the accuracy of the search direction combined with the high condition number of M for large k. The next two pictures show the primal accuracy demand (1 ? 1 )kAxk ? bk2 and the calculation time for Algorithm 1 and Method B: 8
10
6
10
4
10
primal accuracy demand
2
10
0
10
−2
10
−4
10
−6
10
−8
10
−10
10
0
5
10
15
20
25
30
35
25
30
35
iteration
3
calculation time (seconds)
2.5
2
1.5
1
0.5
0
0
5
10
15
20 iteration
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
13
We can see that the time needed for the calculation of the search direction indeed increases for higher k. Because in the nal iterations the calculation time is 30 times as high as the calculation time in the rst iteration, and because of the fact that the calculation of the search direction via Method A needs a constant time (only 0.09 seconds), it makes sense to use Method B only in the rst few iterations and switch to Method A if Method B needs more time than Method A (this method is henceforth called Method C). Using Method C we obtain the following results: Algorithm 1 (Algorithm 2) needs 34 (23) iterations and 3.42 (2.41) seconds, with a total of only 6 (3) search directions calculated via Method B. This result (a small increase in the number of iterations, a small increase or decrease in the total processing time, and very few search directions calculated with Method B) is typical for nearly all problems of the NETLIB-set (see [4]) and can be explained by the fact that the sparse Cholesky decomposition can be calculated very eciently (in MATLAB) for those small problems. Even for the few NETLIB-problems for which Method C calculated at least two search directions with Method B the di erence in the processing time (and the number of iterations) of Method A and Method C is nearly negligible ( 7%). In order to show the performance of our algorithms on large problems, we give results for some slightly modi ed problems of the NETGEN-set (the upper bounds on x were removed). The matrices of these large problems (n = 15:325; 27:887; 52:809; m = 5:000) have a sparsity of less than 0.06%, the calculation of the Cholesky decomposition of the reordered matrix M takes 41{251 seconds. Since, on the other hand, an iteration of the CG-algorithm takes 0.01{0.08 seconds, we expect that the use of Method C will result in a substantial decrease of the needed calculation time. The following tables give the results for " = 1e-07 or 5e-03. (Smaller values for " produced numerical problems when calculating the Cholesky decomposition of M .)
Algorithm 1: Results for NETGEN-problems, " = 1e-07 Method A iterations seconds
101 102 103 105 108 109
38 40 39 37 49 54
7589 8024 7978 7863 12453 12461
Method C iterations seconds total Method B 43 38 2315 45 32 3111 46 38 2647 44 34 2979 56 30 6887 44 31 3397
speed-up
3.28 2.58 3.01 2.64 1.81 3.67
Algorithm 1: Results for NETGEN-problems, " = 5e-03 Method A iterations seconds
104 106 107 110
24 22 32 23
4932 920 7848 943
Method C iterations seconds total Method B 36 36 1163 26 21 433 42 42 1995 26 20 500
speed-up
4.24 2.08 3.93 1.89
14
JANOS KORZAK
Algorithm 2: Results for NETGEN-problems, " = 1e-07 Method A iterations seconds
20 25 23 20 28 26
101 102 103 105 108 109
Method C iterations seconds total Method B 25 21 1494 27 19 1767 25 21 1640 23 20 1492 30 23 2504 31 24 2427
4262 4891 4796 4317 7366 6341
speed-up
2.85 2.77 2.92 2.89 2.94 2.61
Algorithm 2: Results for NETGEN-problems, " = 5e-03 Method A iterations seconds
15 13 17 14
104 106 107 110
Method C iterations seconds total Method B 23 20 1081 21 16 348 22 22 628 16 14 255
3220 597 4480 617
speed-up
2.98 1.72 7.13 2.42
We notice two facts: First, the use of Method C results in an increase in the number of iterations, but because of the high number of search directions that are calculated with Method B a speed-up between 1.72 and 7.13 (usually 3) is reached. Secondly, the shortest running time is reached with Algorithm 2 and Method C. This is somewhat surprising, because in most iterations it is necessary to use the CGalgorithm for the solution of two linear systems. We nally have a look at the time that Algorithm 2 needs to calculate a search direction via Method A and Method B for the problem NETGEN 103: 3
10
Method A
2
calculation time (seconds)
10
1
10
Method B
0
10
−1
10
0
5
10
15
20
25
iteration
The plot shows that in the rst iterations Method B is much faster than Method A. Because the number of iterations of Algorithm 2 increases by only two if Method C is used instead of Method A, the use of Method C results in a huge decrease in processing
Convergence Analysis of Inexact Infeasible-Interior-Point-Algorithms
15
time.
7. Concluding Remarks. In this paper we proved the (polynomial) complexity of a class of inexact infeasible-interior-point-algorithms. This class includes inexact variants of some practically ecient infeasible-interior-point-algorithms, in particular variants of the algorithms of Kojima et. al. and Mehrotra. The theory developed in this paper usually justi es the use of the Cholesky decomposition for determining a search direction, because the calculated search direction, which is aicted with rounding errors, is in most cases an inexact search direction in the sense of this paper. Furthermore we have seen that the use of Krylov subspace methods results in an increase in the number of iterations, and for large problems in a huge decrease of the processing time. To make this kind of calculation more time-ecient, it seems necessary to use a more sophisticated preconditioner and to calculate the searchdirections with the help of stable linear systems if (xk )T z k approaches zero (e.g. one can try to use (M + I ) for some > 0). Another approach is to use smaller step sizes, because the used long step sizes are at least partially the reason for the high condition number of M even for small k. The use of smaller step sizes will result in an increase of the number of iterations of Algorithm 1 and Algorithm 2, but if the step sizes are chosen carefully, this can nevertheless lead to a decrease in the total processing time. Some numerical evidence for this is, that we can achieve a decrease of 22% in processing time for problem israel and Algorithm 1 with Method B, if we set 4 = 0:99 instead of 4 = 0:99995. Acknowledgments. I am grateful to the referees for suggestions and improvements upon the earlier versions of this manuscript. REFERENCES [1] Roland W. Freund, Florian Jarre, A QMR-Based Interior-Point Algorithm for Solving Linear Programs, Numerical Analysis Manuscript 94-19, AT&T Bell Laboratories, December 1994 [2] Roland W. Freund, Florian Jarre, Shinji Mizuno, Convergence of a Class of Inexact Interior-Point Algorithms for Linear Programs, Interior-Point Methods Online, Argonne National Laboratory, http://www.mcs.anl.gov/home/otc/InteriorPoint/, September 1996 [3] Masakazu Kojima, Nimrod Megiddo, Shinji Mizuno, A primal-dual infeasible-interior-point algorithm for linear programming, Mathematical Programming, Vol. 61, pp. 263{280, 1993 [4] Janos Korzak, Primal-Duale pfadfolgende inexakte Auere-Punkte-Verfahren zur Losung linearer Optimierungsaufgaben, Dissertation, Fachbereich Mathematik, Bergische Universitat - Gesamthochschule Wuppertal, Januar 1997 [5] Irvin J. Lustig, Roy E. Marsten, David F. Shanno, On Implementing Mehrotra's PredictorCorrector Interior Point Method for Linear Programming, SIAM Journal on Optimization, Vol. 2, No. 3, pp. 435{449, August 1992 , Computational experience with a globally convergent primal-dual predictor-corrector [6] algorithm for linear programming, Mathematical Programming, Vol. 66, pp. 123{135, 1994 [7] Sanjay Mehrotra, On the Implementation of a Primal-Dual Interior Point Method, SIAM Journal on Optimization, Vol. 2, No. 4, November 1992, pp. 575{601 [8] Shinji Mizuno, Florian Jarre, Global and Polynomial-Time Convergence of an InfeasibleInterior-Point Algorithm Using Inexact Computation, Research Memorandum 605, The Institute of Statistical Mathematics, 4-6-7 Minami-Azabu, Minato-ku, Tokyo 106, Japan, April 1996 [9] , An Infeasible-Interior-Point Algorithm Using Projections onto a Convex Set, Annals of Opererations Research, Vol. 62, pp. 59-80, 1996 [10] Stephan M. Robinson, A Characterization of Stability in Linear Programming, Operations Research, Vol. 25, No. 3, May-June 1977
