corrector method for solving linear programs from ... - Semantic Scholar

Report 2 Downloads 70 Views
Mathematical Programming67 (1994) 383-406

A quadratically convergent predictor-corrector method for solving linear programs from infeasible starting points a. Florian A. Potra Department of Mathematics, Universityof lowa, Iowa City, IA 52240, USA

Received 31 July 1992; revised manuscriptreceived 23 May 1994

Abstract A predictor-corrector method for solving linear programs from infeasible starting points is analyzed. The method is quadratically convergent and can be combined with Ye's finite termination scheme under very general assumptions. If the starting points are large enough then the algorithm has O(nL) iteration complexity. If the ratio between feasibility and optimality at the starting points is small enough then the algorithm has O ( n ~ ) iteration complexity. For feasible starting points the algorithm reduces to the Mizuno-Todd-Ye predictor-corrector method. AMS Subject Classification: 90C05 Keywords: Linear programming;Pfimal-dual; Polynomialcomplexity;Infeasible;Interior-point;Exterior-point

1. Introduction The vast majority of theoretical results on interior-point methods for linear programming assume that the starting point satisfies exactly the equality constraints and is strictly positive, i.e. it lies in the interior of the region defined by the inequality constraints. All subsequent points generated by the interior-point algorithm will have the same properties. However, in practice it may be very difficult to obtain feasible starting points. Numerical experiments have shown that it is possible to obtain good practical performance by using starting points that lie in the interior of the region defined by the inequality constraints but do not satisfy « This work was supported by an interdisciplinary research grant from the Institute for Advanced Studies of the University of lowa. 0025-5610 © 1994--The Mathematical ProgrammingSociety,Inc. All rights reserved SSD10025 -5610 ( 94 ) 00038-U

384

F.A. Potra / Mathematical Programming 67 (1994) 383-406

the equality constraints (cf. [6] ). The points generated by the algorithm will remain in the interior of the region defined by the inequality constraints but will never satisfy exactly the equality constraints, although the measure of "feasibility" as weil as "optimality" will improve at each step. This property is reflected in the name "infeasible-interior-point algorithm" that has been recently suggested for such methods. The first theoretical result on infeasible-interior-point algorithms was obtained by Kojima, Meggido and Mizuno [5]. They showed that an infeasible-interior-point variant of the primal-dual (feasible) interior-point method studied in [ 11 ] is globally convergent. Their method can also be used to detect regions that do not contain optimal solutions. The first polynomial complexity result was obtained by Zhang [ 19] who proved that with proper initialization an infeasible-interior-point algorithm has O(n2L)-iteration complexity. Shortly after that Mizuno [ 10] proved that the Kojima-Meggido-Mizuno algorithm has also O(n2L)-iteration complexity. In [ 10] and [ 14] two globally convergent predictorcorrector infeasible-interior-point algorithms having O(nL)-iterafion complexity are proposed. The direction of the predictor from [ 10] coincides with the direction of the KojimaMeggido-Mizuno algorithm while the corrector coincides with the Mizuno--Todd-Ye [ 11 ] corrector. The algorithm of [ 14] requires one matrix factorization for the predictor and another one for the corrector, as in the original Mizuno-Todd-Ye algorithm and the predictor-corrector method of [ 10]. However two back-solves are used for the predictor and only one for the corrector. This is due to the fact that the predictor needs the computation of two directions, one for improving "optimality" and the other one for improving "feasibility". The step-lengths for the two directions are different but they are chosen in such a way that both "optimality" and "feasibility" are improved at the same rate. We mention that this is not the case with the predictor-corrector method of [ 10], where "optimality" may be improved of a slightly slower rate than "feasibility". The determination of the steplength in [ 14] requires the solution of an eighth order polynomial equation. We have not been able to prove the quadratic convergence of that algorithm without imposing additional conditions. In the present paper we propose a new infeasible-interior-point predictor-corrector algorithm that requires two matrix factorizations and at most three back-solves per iteration step. If the primal and dual search directions of the predictor are orthogonal, which is always the case with feasible starting points, only two back-solves are needed. In fact for feasible starting points our algorithm reduces to the Mizuno-Todd-Ye predictor-corrector algorithm [ 11 ]. If the primal and dual search directions of the predictor are not orthogonal then in contrast with the algorithm from [ 14] we use only one back-solve for the predictor and two for the corrector. The step-length is determined by solving at most five quadratic equations and therefore is given in explicit form. The algorithm shares with the algorithm from [ 14] the property that feasibility and optimality are improved at the same rate at each iteration, so that the balance between feasibility and optimality existing at the starting point is preserved along the trajectory. Hence the two algorithms are "balanced" in the sense of Freund [2]. We show that our algorithm is globally and quadratically convergent when started from

F.A. Potra / Mathematical Programming 67 (1994) 383-406

385

any positive primal-dual pair that is approximately centered under the sole assumption that an optimal primal-dua_a_l_pair exists. According to Ostrowski [ 13 ] the asymptotic efficiency index is defined as ~/ord, where ord is the order of the iterative procedure and the c is the cost of an iteration. As the main cost associated with interior-point methods is represented by matrix factorizations, and since our algorithm requires two matrix factorizations per step it is natural to take c = 2. Because in our case ord = 2 it follows that the asymptotic efficiency index of out infeasible-interior-point method is equal to X/2. We also show that our algorithm can be combined with Ye's finite termination scheme [ 17]. The computational complexity or our infeasible-interior-point method is expressed in terms of the ratio between "feasibility" and "optimality" at the starting point. If this ratio is smaller than a certain bound then the algorithm has O(~nL)-iteration complexity, the same as the best iteration complexity of any (feasible) interior-point algorithm. For starting points that are "large enough" the algorithm has O (nL)-iteration complexity. We note that the other theoretical results mentioned above assume that the starting points have all their components equal and "large enough". On the other hand, as shown by recent numerical experiments [ 7 ], infeasible-interior-point algorithms are especially efficient when " w a r m started". The results of the present paper could be used to explain this phenomenon either by the lower complexity associated with starting points with favorable " f e a s i b i l i t y " / "optimality" ratio, or by the quadratic convergence results that show that if the primal dual gap is small enough then quadratic convergence prevails. To our knowledge the algorithm presented in this paper is the first infeasible-interior-point algorithm with quadratic convergence. For feasible interior-point algorithms with quadratic convergence see [8] and [ 18]. The present paper supersedes the technical reports [ 15] and [ 16] which are combined here for publication in one comprehensive and self-sufficient paper. The notation used throughout the paper is rather standard: capital letters denote matrices, script capital letters denote sets, lower case letters denote vectors, and Greek letters denote scalars. All vectors are considered to be column vectors. The components of a vector u ~ ~ n will be denoted by [u]» i = 1 . . . . . n. The relation u > 0 is equivalent to [ u ] i > 0, i = 1 . . . . . n, while u >/0 means [ u ] i >~0, i = 1 ..... n. If u ~ R n, w ~ R '~ then ( u, w ) denotes the column vector formed by the components of u and w, i.e. (u, w) ~~n+m, [(U, W)]i = [U]i for l 0, g > 0. Hence (2, g) ~ ~ a . With z + defined by (A10) we can write

(42)

392

F.A. Potra / Mathematical Programming 67 (1994) 383--406

x +s + = 2 g + (Yü+20) + ü6 =~fg+ ( £ü +.~O) + ö~( #ü + üO) + üO =/2e+ Ö2( Üü-~ ~Ü) -~ ~Ô .

(43)

By using the fact that ffrü = ürü = ü~O = 0 we get

/z+ = # = ( i - ö)~.

(44)

Therefore, according to (12), (14), (33), (42)-(44), Lemmas 2.2 and 2.3 we have

II#ü + ü011< II#üll + Ilüell = II(äü) (d- '~ II+ II(dü) (ä-'0)II 1

~< ~

(Iläall IIg - 'õll + Iläüll Ilg - '~11)

1

~< ~(lldüll2+llä-'all

2)

'/2(ll,~üll

2+

IIg-'Ollb'/2

,

and

IIx + s + - ~ + e II < õ" II Õü + ü,~ll + II aOll Ò, < ~ ( Iläüll ~ + lid- 'Sll ") ';"( Ildüll ~ + lid- 'Oll ~) ,,2 1

+~--~ (lldüll ~ + IId-'Oll 2)

~ Il(Xs-) -' ~q,{õ2 uv lll~/2e- ecu + Il#~ell 2)_,__ ,

0, independent of n, such that

(84)

õk >1 1/ (~n ~') •

In what follows we will prove that there is a constant 31, which may depend on n such that

ök >~l - y t z k ,

k>~O.

(85)

The above relation, together with (30) implies Q-quadratic convergence, in the sense that tzk+l~Y~~,

IIpk+ll]0, i = 1, 2 . . . . . n} .

(87)

It is weil known that there is a unique partition B t J ù 4 z = { 1, 2 . . . . . n},

~'NAz=0,

(88)

such that for any (x,s,y) ~ ~ c wehave ( [ x ] i > 0 , [ s ] i = 0 , V i ~ ~ ) and ( [ x ] i = 0 , I s ] i > 0 , Vi ~ J ) . This means that the "basic" and "non-basic" variables are invariant for any strictly-complementary solution. Let us denote by B and N t h e corresponding partition of the columns of A. Also for any vector v ~ R n we denote by v» the vector of components [ v] » i ~ :~, and by VN the vector of components [ v ] » i ~ Az. Then the optimal face for the primal is

g2p = {x~ R": Bx» =b, xB >10, X N =0}

(89)

and the optimal face for the dual is g2d = {(y, s) ~ ~,n+ù: SB = C B - - B l y = 0 , sN =cN--NTy~O} •

(90)

Define ~p = min m a x [ x ] j , j~s~ x~g2p

~a = min max [s]1, j ~ ~ (s,y)~g2a

~=min{~p, ~a}.

(91)

Obviously, ~> 0. We will show that there is a number ~> 0 such that the sequence generated by Algorithm 2.1 satisfies

x~>~~e,

s~>~~e,

k>~l.

(92)

F.A. Potra / Mathematical Programming 67 (1994) 383-406

400

Lemma 4.1. Let (x*, s*, y * ) ~ J « and denote ~* =min{[ [x*]i: i ~ ~ } ,

«* =min{ [s*]i: i~JY} ,

(93)

(1-a)~* ~*=min{~*, ~J'},

~=

n((2_Õo)/Õo+C, ) ,

(94)

where «* = ( (x ° ) ~'s* + (s o ) ~ x * ) / ( (x ° ) ~s °) .

(95)

Then (92) is satisfied. Proof. We drop the index k and use the notation from the proof of Lemma 3.1. According to (55), (30) and (31) we have

s~'x * +x~'s * ( 1 - o t ) [ x * ] i /> ( 1 - a ) ~ * n÷

~>~,

Vi~~,



which proves the first inequality in (92). The second inequality can be proved in a similar manner. [] Before stating our next lemma let us introduce some useful notation that will be used in the proof, as weil as later on in the paper. Let B be any set of rows of B having maximal rank (/~ = B if the rows of B are independent). Let ~r, Ä, and 6 denote the corresponding rows of N and A, and components of b. By introducing the rectangular matrix P whose rows are formed by the unit vectors of R '~ corresponding to/~ we can write

F.A. Potra/ MathematicalProgramming67 (1994) 383-406 B=PB ,

A=PA ,

N=PN,

b=Pb.

401 (98)

Now, by generalizing some of the techniques developed in [ 18] we can prove the following result: L e m m a 4.2. There is a constant ~-~ (0, 1), independent of k, such that if u k, v k are the vectors produced in sub-step (A3) of Algorithm 2.1, then

Ilukll ~< ~-~k,

Ilvkll ~< ~-~»

k>~0.

(99)

Proof. For simplicity we drop the index k, and use the notation of Algorithm 2.1 as weil as (16) and (17). From Theorem 2.5 we have (x, s) ~ ù4~«so that by using (92) we can write ( l + a ) t z 1/2

IIdB Il~ ~ 1-- Y6~k,

403

F.A. Potra / Mathematical Programming 67 (1994) 383-406

while from (A7) we deduce that O3>~ ~

/3

1

>~l-- ~ ß ô k •

By applying Lemma 4.2 we deduce that (85) holds with T= max(Ts, ~/6, ~ßB) "r2. Then according to Theorem 2.5 we deduce that (86) is satisfied. By taking .r/= maùx{31, To, ,/q } we obtain the following quadratic convergence result. Theorem 4.3. I f Y * is nonempty, then there is a constant ~ independent of k such that

tzk+~ 0.

[]

In all the theorems p r o v e n so far we have a s s u m e d that Y * is nonempty. It turns out that A l g o r i t h m 2.1 can be modified in such a way that it can detect whether 9 - * contains points of n o r m less than a quantity chosen in advance.

Theorem 5.4. Suppose that the following instruction " ( A 2 . 5 ) if (x °) rsk + (s °) rxk>~ n(tz« + tZo) + ( 1 - tzk/I.to) (Pd IIx°ll minate"

+ t~,, IIs°ll )

then ter-

is inserted in between instructions ( A 2 ) and ( A3 ) o f Algorithm 2.1. Then the new algorithm terminates either at ( A 2 ) with z ~ g - « or at ( A 2 . 5 ) , and in the latter case there is no Z*=(x*,s*,y

*) ~ ~ - * suchthat

IIx*ll